FUNCTION dustem_compute_stokext,p_dim,st,$
QEXT_spec=QEXT_spec,$
UEXT_spec=UEXT_spec,$
PSIEXT_spec=PSIEXT_spec,$
dustem_psi_ext=dustem_psi_ext,$
out_st=out_st,$
_extra=extra

;+
; NAME:
;    dustem_compute_stokext
; PURPOSE:
;    Provides Stokes Q and Stokes U "extinction" optical depths for a given dustem spectrum and wrapper SED 
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    out=dustem_compute_stokext(p_dim[,st=][,_extra=][,/help])
; INPUTS:
;    p_dim      = parameter values
; OPTIONAL INPUT PARAMETERS:
;    stp        = Dustem structure
; OUTPUTS:
;
;
;
;
; OPTIONAL OUTPUT PARAMETERS:
;    stp        = Dustem output structure
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.

IF keyword_set(help) THEN BEGIN
    doc_library,'dustem_compute_stokext'
    out=0.
    goto,the_end
ENDIF

IF not keyword_set(st) THEN BEGIN

    dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st
    ;st=dustem_run(p_dim)

ENDIF

out=0. ; no specific output is required.

If keyword_set(result) THEN result=st

degtorad = !pi/180

EXT_spec = st.ext.ext_tot * (*!dustem_HCD)/1.e21     ;This is Total extinction (for a default value of 10^20 NH/cm^2)
POLEXT_spec = st.polext.ext_tot * (*!dustem_HCD)/1.e21   ;This is Polarized extinction

;Hard-coded test to set to zero potential negative values in the polarized extinction arrays:
;so far this has been seen in absorption arrays  
;This test needs to include all grain species that polarize

IF !run_pol THEN BEGIN
    
    ;THIS FALSE CORRECT  IT USING POLEXT_spec
    ;POLEXT_spec.abs_grain[2,*]
    ind_ngtv_plxt = where (POLEXT_spec LT 0, ct_ngtv_plxt)
    
    IF ct_ngtv_plxt NE 0 THEN (POLEXT_spec)[ind_ngtv_plxt] = 0.
        
        
ENDIF

Nwaves=(size(EXT_SPEC))[1]

frac=POLEXT_spec/EXT_spec
tes=where(finite(frac) eq 0)
frac(tes)=0.

scopes=tag_names((*!dustem_plugin))
IF scopes[0] NE 'NONE' THEN BEGIN

    FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN   
        
        IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'REPLACE_POLEXT') THEN BEGIN
            
            QEXT_spec = (*(*!dustem_plugin).(i).spec)[*,1]
            UEXT_spec = (*(*!dustem_plugin).(i).spec)[*,2]

            
        ENDIF 
    
    ENDFOR
ENDIF    

    
IF ~isa(QEXT_spec) && ~isa(UEXT_spec) THEN BEGIN

    polar_ippsi2iqu,EXT_spec,QEXT_spec,UEXT_spec,frac,replicate(!dustem_psi_ext,Nwaves);PSI_spec 

ENDIF    
        

IF scopes[0] NE 'NONE' THEN BEGIN        
    FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN
        
        IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'ADD_POLEXT') THEN BEGIN
            
            QEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,1]
            UEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,2]

        ENDIF

    ENDFOR
ENDIF

;-----------------------------------------

;Hard-coded block to avoid cases where Q=-0.0 and U=0.
testnzero = where(QEXT_spec EQ -0., ct_nzero)
if ct_nzero ne 0 then QEXT_spec[testnzero] = 0

PSIEXT_spec = 0.5*atan(UEXT_spec,QEXT_spec)/degtorad

;INITIALIZING THE STOKES SEDS

dustem_qext = (*(*!dustem_data).qext).values * 0.
dustem_uext = (*(*!dustem_data).uext).values * 0.

;Using a list because dustem_qext and dustem_uext can have different dimensions
out=list(fltarr(n_elements(dustem_qext)),fltarr(n_elements(dustem_uext)))

if not isarray(EXT_spec) THEN stop 


;COLOR CORRECTION NOT PERFORMED YET. ONLY SPECTRUM DATA POINTS ARE CONSIDERED FOR NOW
; 
; ;Performing color corrections on the dustem spectra and generating the SEDs at the given filters 
; 
; IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
;     message,'DOING color correction calculations',/info
; ENDIF ELSE BEGIN
;           message,'SKIPPING color correction calculations',/info
;       ENDELSE
; 

; ind_qsed=where((*(*!dustem_data).qsed).filt_names NE 'SPECTRUM',count_qsed)
; ind_used=where((*(*!dustem_data).used).filt_names NE 'SPECTRUM',count_used)

;   
;   
; IF count_qsed NE 0 THEN BEGIN 

;   filter_names=((*(*!dustem_data).qsed).filt_names)(ind_qsed)    

;   if isa(dustem_qsed) then begin 
;     sqsed=dustem_cc(st.polsed.wav,Q_spec,filter_names,cc=cc)
;     dustem_qsed(ind_qsed)=sqsed
;   endif
; ENDIF
; 
; IF count_used NE 0 THEN BEGIN

; filter_names=((*(*!dustem_data).used).filt_names)(ind_used)

; if isa(dustem_used) then begin
;   sused=dustem_cc(st.polsed.wav,U_spec,filter_names,cc=cc)
;   dustem_used(ind_used)=sused
; endif
; ENDIF


;For spectrum data points, interpolate in log-log
;Linear interpolation leads to wrong values, in particular where few
;wavelengths points exist in the model (long wavelengths).

ind_qspec=where((*(*!dustem_data).qext).filt_names EQ 'SPECTRUM',count_qspec)
ind_uspec=where((*(*!dustem_data).uext).filt_names EQ 'SPECTRUM',count_uspec)

IF count_qspec NE 0 THEN BEGIN

  if isa(dustem_qext) then dustem_qext(ind_qspec)=interpol(QEXT_spec,st.polext.wav,(((*(*!dustem_data).qext).wav)(ind_qspec)))  

ENDIF

IF count_uspec NE 0 THEN BEGIN

  if isa(dustem_uext) then dustem_uext(ind_uspec)=interpol(UEXT_spec,st.polext.wav,(((*(*!dustem_data).uext).wav)(ind_uspec)))

ENDIF

out_st=st

;generating interpolates for Psi_ext
dustem_psi_ext = 0.5*atan(dustem_uext,dustem_qext)/degtorad

out[0] = dustem_qext
out[1] = dustem_uext

;clean pointers
heap_gc

the_end:

RETURN,out

END