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