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