FUNCTION dustem_compute_stokext,p_dim,$
                                st=st,$
                                QEXT_spec=QEXT_spec,$
                                UEXT_spec=UEXT_spec,$
                                PSIEXT_spec=PSIEXT_spec,$
                                dustem_psi_ext=dustem_psi_ext,$
                                help=help
                               

;+
; NAME:
;    dustem_compute_stokext
;
; PURPOSE:
;    Computes Stokes Parameters for extinction in units of optical depth for a given Dustem spectrum
;
; CATEGORY:
;    Dustem
;
; CALLING SEQUENCE:
;    out=dustem_compute_stokext(p_dim[,st=][,QEXT_spec=QEXT_spec][,UEXT_spec=UEXT_spec][,PSIEXT_spec=PSIEXT_spec][,dustem_psi_ext=dustem_psi_ext][,/help])
;
; INPUTS:
;    p_dim      = parameter values
;
; OPTIONAL INPUT PARAMETERS:
;    st        = Dustem structure
;
; OUTPUTS:
;    out        = List Containing the Stokes parameters extinction for spectrum points in !dustem_data 
;    out[0]     = Computed Stokes Q extinction for spectrum points in !dustem_data
;    out[1]     = Computed Stokes U extinction for spectrum points in !dustem_data 
;
; OPTIONAL OUTPUT PARAMETERS:
;    st        = Dustem output structure
;    QEXT_spec        = Dustem Stokes Q extinction 
;    UEXT_spec        = Dustem Stokes U extinction
;    PSIEXT_spec      = Dustem Polarization angle output in extinction
;    dustem_psi_ext   = Computed Polarization angle ouput for spectrum points in !dustem_data
;
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
;
; PROCEDURE:
;    None
;
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard 2022
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-


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
    ;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.

plugin_tags=(*!dustem_plugin).name
ind=where(plugin_tags NE 'NONE',Nplugins)
IF Nplugins NE 0 THEN BEGIN
    FOR i=0L,Nplugins-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 Nplugins NE 0 THEN BEGIN
    FOR i=0L,Nplugins-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