FUNCTION dustem_compute_polext,p_dim,$ st=st,$ POLEXT_spec=POLEXT_spec,$ SPEXT_spec=SPEXT_spec,$ dustem_fpolext=dustem_fpolext,$ dustem_ext=dustem_ext,$ help=help ;+ ; NAME: ; dustem_compute_polext ; ; PURPOSE: ; Computes Polarized Extinction in units of optical depth for a given Dustem spectrum ; ; CATEGORY: ; Dustem, Distributed, ; ; CALLING SEQUENCE: ; polext=dustem_compute_polext(p_dim[,st=][,polext_spec=][,spext_spec=][,dustem_fpolext=][,dustem_ext=][,/help]) ; ; INPUTS: ; p_dim = parameter values ; ; OPTIONAL INPUT PARAMETERS: ; st = Dustem output structure ; dustem_ext = Computed polarized extinction for spectrum points in !dustem_data ; ; OUTPUTS: ; polext = computed Polarized Extinction for spectrum points in !dustem_data ; ; OPTIONAL OUTPUT PARAMETERS: ; st = Dustem output structure ; polext_spec = Dustem polarized extiction output ; spext_spec = Dustem polarized exinction fraction ; dustem_fpolext = Computed polarized extinction fraction for spectrum points in !dustem_data ; ; ACCEPTED KEY-WORDS: ; help = If set, print this help ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; None ; ; RESTRICTIONS: ; The DustEMWrap IDL code must be installed ; ; MODIFICATION HISTORY: ; Written by JPB,NF,DP Jan-2007 ; 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_polext' dustem_polext=0. goto,the_end ENDIF IF not keyword_set(st) THEN BEGIN dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st ENDIF POLEXT_spec = st.polext.ext_tot * (*!dustem_HCD)/1.0e21 ; ;Hard-coded test to set to potential zero 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 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 EXT_spec = st.ext.ext_tot * (*!dustem_HCD)/1.0e21 ;This is Total intensity I out=0. Nwaves=(size(POLEXT_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 ;IF ptr_valid(!dustem_plugin) THEN BEGIN for i=0L,n_tags(*!dustem_plugin)-1 do begin if total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) eq 'ADD_EXT') then EXT_spec+=(*(*!dustem_plugin).(i).spec)[*,0] endfor 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 '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) SPEXT_spec = frac 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 POLEXT_spec=sqrt(QEXT_spec^2+UEXT_spec^2) SPEXT_spec = POLEXT_SPEC/EXT_spec dustem_polext = (*(*!dustem_data).qext).values * 0. if not isarray(POLEXT_spec) THEN stop ;I don't understand this test. The only thing it can indicate is a problem with the dust parameters or the fortran executable. But in my opinion these tests would be ;NO COLOR CORRECTION FOR EXTINCTION. FOR NOW ... ; 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_polsed=where((*(*!dustem_data).qsed).filt_names NE 'SPECTRUM',count_polsed) ; ; IF count_polsed NE 0 THEN BEGIN ; filter_names=((*(*!dustem_data).qsed).filt_names)(ind_polsed) ; spolsed=dustem_cc(st.polsed.wav,P_spec,filter_names,cc=cc) ; dustem_polsed[ind_polsed]=spolsed ; 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). ;I believe this to be automatic because the provided dustem spectra are already sampled on a log-log grid. ind_spec=where((*(*!dustem_data).qext).filt_names EQ 'SPECTRUM',count_spec) IF count_spec NE 0 THEN dustem_polext(ind_spec)=interpol(POLEXT_spec,st.polext.wav,(((*(*!dustem_data).qext).wav)(ind_spec))) ;GENERATING THE INTERPOLATES FOR POLFRAC (dustem_polfrac) if !run_lin then begin if not keyword_set(dustem_ext) then dustem_ext = dustem_compute_ext(p_dim,st=st,EXT_spec=EXT_spec) If n_elements(dustem_ext) ne n_elements(dustem_polext) then begin if n_elements(dustem_ext) gt n_elements(dustem_polext) then begin dustem_polext_x = dustem_polext dustem_ext_x = dustem_polext ;meaning dustem_sed needs to be modified nwaves = n_elements(dustem_polext) for i=0L,nwaves-1 do begin j=where((*(*!dustem_data).ext).wav EQ (*(*!dustem_data).fpolext).wav(i),testwav) if testwav ne 0 then dustem_ext_x(i) = dustem_ext(j(0)) endfor endif ELSE begin dustem_polext_x = dustem_ext dustem_ext_x = dustem_ext ;meaning dustem_polsed needs to be modified nwaves = n_elements(dustem_ext) for i=0L,nwaves-1 do begin j=where((*(*!dustem_data).ext).wav EQ (*(*!dustem_data).fpolext).wav(i),testwav) if testwav ne 0 then dustem_polext_x(i) = dustem_polext(j(0)) endfor ENDELSE endif ELSE BEGIN dustem_polext_x = dustem_polext dustem_ext_x = dustem_ext ENDELSE dustem_fpolext = dustem_polext_x/dustem_ext_x ENDIF the_end: return, dustem_polext END