FUNCTION dustem_compute_polsed,p_dim,$ st=st,$ P_spec=P_spec,$ SP_spec=SP_spec,$ dustem_polfrac=dustem_polfrac,$ dustem_sed=dustem_sed,$ help=help ;+ ; NAME: ; dustem_compute_polsed ; ; PURPOSE: ; Computes a Polarized SED for a given Dustem spectrum ; ; CATEGORY: ; DustEMWrap, Mid-Level, Distributed ; ; CALLING SEQUENCE: ; sed=dustem_compute_polsed(p_dim[,st=][,P_spec=][,SP_spec=][,dustem_polfrac=][,dustem_sed=][,/help]) ; ; INPUTS: ; p_dim = parameter values ; ; OPTIONAL INPUT PARAMETERS: ; st = Dustem structure ; dustem_sed = Computed emission for spectrum points in !dustem_data ; ; OUTPUTS: ; sed = computed SED for filters in !dustem_data ; dustem_polsed = Computed polarized emission for spectrum points in !dustem_data ; P_spec = Dustem polarized emission output ; SP_spec = Dustem polarized emission fraction ; dustem_polfrac = Computed polarization emission fraction for spectrum points in !dustem_data ; ; OPTIONAL OUTPUT PARAMETERS: ; st = Dustem structure ; dustem_sed = Computed emission 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 ; 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_polsed' dustem_polsed=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 If keyword_set(result) THEN result=sti fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.e20/1.e7 P_spec = st.polsed.em_tot*fact ;This is Polarized intensity P stI = st.sed.em_tot*fact ;This is Total intensity I out=0. Nwaves=(size(P_spec))[1] frac=P_spec/stI tes=where(finite(frac) eq 0) frac(tes)=0. ;This block was added to generate the polarization fraction spectrum with the plugins contributions. ;smallp (polfrac) is naturally not an additive quantity. ;NB: if plugins that replace the total emission given by the fortran executable exist, ;the need to add a block with a 'REPLPACE_SED' scope as for the 'REPLACE_POLSED' below, may arise. ;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 plugin_tags=(*!dustem_plugin).name ind=where(plugin_tags NE 'NONE',Nplugins) IF Nplugins NE 0 THEN BEGIN ;IF ptr_valid(!dustem_plugin) THEN BEGIN for i=0L,Nplugins-1 do begin if total(strsplit(((*!dustem_plugin)[i].scope),'+',/extract) eq 'ADD_SED') then stI+=(*(*!dustem_plugin)[i].spec)[*,0] endfor ENDIF IF Nplugins NE 0 THEN BEGIN FOR i=0L,Nplugins-1 DO BEGIN IF total(strsplit(((*!dustem_plugin)[i].scope),'+',/extract) EQ 'REPLACE_POLSED') THEN BEGIN Q_spec=(*(*!dustem_plugin).(i).spec)[*,1] U_spec=(*(*!dustem_plugin).(i).spec)[*,2] ;make sure sti hasn't been replaced. But this is not happening in the near future ENDIF ENDFOR ENDIF IF ~isa(Q_spec) && ~isa(U_spec) THEN BEGIN polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(!dustem_psi,Nwaves) ;the user needs to set this pointer somewhere SP_spec = frac ENDIF IF Nplugins NE 0 THEN BEGIN FOR i=0L,Nplugins-1 DO BEGIN IF total(strsplit(((*!dustem_plugin)[i].scope),'+',/extract) EQ 'ADD_POLSED') THEN BEGIN Q_spec+=(*(*!dustem_plugin)[i].spec)[*,1] U_spec+=(*(*!dustem_plugin)[i].spec)[*,2] ENDIF ENDFOR ENDIF P_spec=sqrt(Q_spec^2+U_spec^2) SP_spec=P_spec/Sti dustem_polsed = (*(*!dustem_data).qsed).values * 0. if not isarray(P_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 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).qsed).filt_names EQ 'SPECTRUM',count_spec) IF count_spec NE 0 THEN dustem_polsed(ind_spec)=interpol(P_spec,st.polsed.wav,(((*(*!dustem_data).qsed).wav)(ind_spec))) ;GENERATING THE INTERPOLATES FOR POLFRAC (dustem_polfrac) if !run_lin then begin ;stop if not keyword_set(dustem_sed) then dustem_sed = dustem_compute_sed(p_dim,st=st,SED_spec=SED_spec) If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin dustem_polsed_x = dustem_polsed dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified nwaves = n_elements(dustem_polsed) for i=0L,nwaves-1 do begin j=where((*(*!dustem_data).sed).wav EQ (*(*!dustem_data).polfrac).wav(i),testwav) if testwav ne 0 then dustem_sed_x(i) = dustem_sed(j(0)) endfor endif ELSE begin dustem_polsed_x = dustem_sed dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified nwaves = n_elements(dustem_sed) for i=0L,nwaves-1 do begin j=where((*(*!dustem_data).sed).wav EQ (*(*!dustem_data).polfrac).wav(i),testwav) if testwav ne 0 then dustem_polsed_x(i) = dustem_polsed(j(0)) endfor ENDELSE endif ELSE BEGIN dustem_polsed_x = dustem_polsed dustem_sed_x = dustem_sed ENDELSE dustem_polfrac = dustem_polsed_x/dustem_sed_x ENDIF ;clean pointers heap_gc the_end: RETURN,dustem_polsed END