FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=_extra ;+ ; NAME: ; dustem_mpfit_run_vg ; PURPOSE: ; function used to fit dustem model with dustem_mpfit_sed ; CATEGORY: ; Dustem ; CALLING SEQUENCE: ; res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help]) ; INPUTS: ; x = necessary for interface but unused ; params = Dustem Model parameters normalized to initial values ; OPTIONAL INPUT PARAMETERS: ; cc_sed = SED color correction factors ; _extra = extra parameters for the plot routine ; OUTPUTS: ; None ; OPTIONAL OUTPUT PARAMETERS: ; None ; ACCEPTED KEY-WORDS: ; help = If set, print this help ; COMMON BLOCKS: ; None ; SIDE EFFECTS: ; SED and model are plotted if !dustem_show_plot ; RESTRICTIONS: ; The dustem idl wrapper must be installed ; PROCEDURE: ; None ; EXAMPLES ; res=dustem_mpfit_run_vg(x,params[,_EXTRA=extra][,cc_sed=cc_sed][,/help]) ; ODIFICATION HISTORY: ; Adapted to extinction and polarization by Vincent Guillet (IAS) ; Written by J.-Ph. Bernard ; see evolution details on the dustem cvs maintained at CESR ; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems. ;- ;note that x is not used. IF keyword_set(help) THEN BEGIN doc_library,'dustem_mpfit_run' sed=0. goto,the_end ENDIF p_dim=p_min*(*(*!dustem_fit).param_init_values) IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac() ;RUN THE MODEL AND GET THE SPECTRUM dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values) st = dustem_run(p_dim) st_model=dustem_read_all(!dustem_soft_dir) chi2_sed = 0 chi2_ext = 0 chi2_polext = 0 chi2_polsed = 0 chi2_polfrac = 0 chi2_qsed =0 chi2_used =0 chi2_qext =0 chi2_uext =0 rchi2_sed = 0 rchi2_ext = 0 rchi2_polext = 0 rchi2_polsed = 0 rchi2_polfrac = 0 rchi2_qsed =0 rchi2_used =0 rchi2_qext =0 rchi2_uext =0 wrchi2_sed = 0 wrchi2_ext = 0 wrchi2_polext = 0 wrchi2_polsed = 0 wrchi2_polfrac = 0 wrchi2_qsed =0 wrich2_used =0 wrchi2_qext =0 wrich2_uext =0 n_sed = 0 n_ext = 0 n_polext = 0 n_polsed = 0 n_polfrac = 0 n_qsed =0 n_used =0 n_qext =0 n_uext =0 win = 0 tagnames=tag_names(!dustem_data) dustem_out = [0] ind_data_tag=dustem_data_tag(count=count_data_tag) ;if exist('./noplot') then !dustem_show_plot = 0 else if !dustem_show_plot eq 0 then !dustem_show_plot = 1 ;stop FOR i_tag=0,count_data_tag-1 DO BEGIN ;stop ;print,i_tag,ind_data_tag(i_tag) ;print,tagnames(ind_data_tag(i_tag)) CASE tagnames(ind_data_tag(i_tag)) OF 'SED' : BEGIN dustem_sed = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron) chi2_sed = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2) n_sed = n_elements((*!dustem_data.sed).values) rchi2_sed = chi2_sed / n_sed wrchi2_sed = rchi2_sed * !fit_rchi2_weight.sed message,"chi2 SED = "+strtrim(chi2_sed,2)+" rchi2 SED = "+strtrim(rchi2_sed,2),/continue;," weighted rchi2 SED=",wrchi2_sed ;print,'wl ',(*!dustem_data.sed).wav ;print,"Delta SED per wl: ",(*!dustem_data.sed).values-dustem_sed ;print,"error per wl: ",(*!dustem_data.sed).sigma ;print,"chi2 per wl: ",((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2 alpha = total(((*!dustem_data.sed).values*dustem_sed)/(*!dustem_data.sed).sigma^2) $ / total(((*!dustem_data.sed).values)^2/(*!dustem_data.sed).sigma^2) beta = total(((*!dustem_data.sed).values*dustem_sed)/(*!dustem_data.sed).sigma^2) $ / total(dustem_sed^2/(*!dustem_data.sed).sigma^2) print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta (*!dustem_fit).chi2=chi2_sed (*!dustem_fit).rchi2=rchi2_sed ; Plot if needed IF !dustem_show_plot NE 0 THEN BEGIN ;title='DUSTEM WRAP('+strtrim(string(win),2)+')' window,win,title='DUSTEM WRAP (SED)' ;cgwindow,WMULTI=[2,2,1],/CURRENT ;I.C ;undeprecated the use of polfrac here ;********** ;I added this block so I can plot polfrac inside the sed window ;I will probably remove it ;********** ;test if polfrac is in tagnames: ;another problem is that sed and pol data don't have the same filters. varvar=where(tagnames eq 'POLFRAC', testfpol) if testfpol ne 0 then begin if !run_lin then begin dustem_polsed = dustem_compute_polsed(p_dim,st) ; this quantity hasn't been computed yet that is why this line is present here. ind_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt) ;locating the filters of polfrac=smallp in SED xcat ind_fspec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',count_fspec) ;locating the spectrum points in polfrac data. Later testing will have to be on the wavelengths ;dustem_sed and dustem_polsed might not have the same length thus this block. ;I'm not so sure about this... 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 ;matching the filter data points IF count_filt ne 0 then begin for i=0L,count_filt-1 do begin j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt) if testfilt ne 0 then dustem_sed_x(ind_filt(i)) = dustem_sed(j(0)) endfor endif ;matching the spectrum data points IF count_fspec ne 0 then begin for i=0L,count_fspec-1 do begin j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec) if testspec ne 0 then dustem_sed_x((ind_fspec(i))) = dustem_sed(j(0)) endfor endif endif ELSE begin ;I think I will have errors here!! The thing is I'm trying to make this work without (*!dustem_data.polsed) which is hard. dustem_polsed_x = dustem_sed dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified ;matching the filter data points IF count_filt ne 0 then begin for i=0L,count_filt-1 do begin j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt) if testfilt ne 0 then dustem_polsed_x(ind_filt(i)) = dustem_polsed(j(0)) endfor endif ;matching the spectrum data points IF count_fspec ne 0 then begin for i=0L,count_fspec-1 do begin j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec) if testspec ne 0 then dustem_polsed_x((ind_fspec(i))) = dustem_polsed(j(0)) endfor endif ENDELSE endif ELSE BEGIN dustem_polsed_x = dustem_polsed dustem_sed_x = dustem_sed ENDELSE dustem_polfrac = dustem_polsed_x/dustem_sed_x ;stop ;stop chi2_polfrac = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2) n_polfrac = n_elements((*!dustem_data.polfrac).values) rchi2_polfrac = chi2_polfrac / n_polfrac wrchi2_polfrac = rchi2_polfrac * !fit_rchi2_weight.polfrac message,"chi2 POLFRAC = "+strtrim(chi2_polfrac,2)+" rchi2 POLFRAC = "+strtrim(rchi2_polfrac,2),/continue ;," weighted rchi2 POLFRAC=",wrchi2_polfrac dustem_out = [ dustem_out, dustem_polfrac ] fpol=dustem_polfrac ; problem because extra structure does not get updated so you actually need the 'extra' procedure you haven't finihed coding. endif endif ;stop dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2,fpol=fpol ;stop ENDIF dustem_out = [ dustem_out, dustem_sed ] END 'EXT' : BEGIN ;ADDING PLUGIN TO SPECTRUM---------------- IF ptr_valid(*!dustem_plugin) THEN BEGIN for i=0L,n_tags(*!dustem_scope)-1 do begin if total(strsplit(*(*!dustem_scope).(i),'+',/extract) eq 'ADD_EXT') then st.ext.ext_tot+=(*(*!dustem_plugin).(i)) endfor endif ;------------------------------------------ dustem_ext = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav) chi2_ext = total(((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2) n_ext = n_elements((*!dustem_data.ext).values) rchi2_ext = chi2_ext / n_ext wrchi2_ext = rchi2_ext * !fit_rchi2_weight.ext print,"chi2 EXT = ",chi2_ext," rchi2 EXT = ",rchi2_ext;," weighted rchi2 EXT=",wrchi2_ext ;print,"chi2 per wl: ",((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2 ; Plot if needed IF !dustem_show_plot NE 0 THEN BEGIN win+=1 ;dustem_plot_ext,st,yr=[1.e-6,1.e0],/ysty,xr=[0.5,40],/xsty,win=win xr=[0.3,40] yr=[1.e-26,1.e-21] dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1 dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2 win+=1 ;dustem_plot_ext,st,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty,win=win xr=[0.2,10] yr=[0,3e-21] dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1 dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2 ENDIF ;dustem_plot_ext,st,/print_Rv ;commented this recently dustem_out = [ dustem_out, dustem_ext ] END 'POLEXT': BEGIN ;---------------the procedure within this block has not been updated for a vert long time IF !dustem_show_plot NE 0 and !run_circ then begin win+=1 dustem_plot_circ,st,win=win,xr=[0.1,5],aligned=st_model.align.grains.aligned,yr=[-0.03,0.03] endif ;-------------------------------------------------------------------------------------- IF !run_lin then begin ;compute_polext needs updating. It's been months. dustem_polext = dustem_compute_polext(p_dim,st,dustem_qext,dustem_uext,Q_ext,U_ext) chi2_polext = total(((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2) n_polext = n_elements((*!dustem_data.polext).values) rchi2_polext = chi2_polext / n_polext wrchi2_polext = rchi2_polext * !fit_rchi2_weight.polext print,"chi2 POLEXT = ",chi2_polext," rchi2 POLEXT = ",rchi2_polext;," weighted rchi2 POLEXT=",wrchi2_polext ;print,"chi2 per wl: ",((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2 ; Fit in (pmax,lmax,K) x=st.polext.wav logx=alog(x) logy=alog(st.polext.ext_tot) ; Whittet+1992 (Table 1) Uband = 0.36 Bband = 0.43 Vband = 0.55 Rband = 0.63 Iband = 0.78 Jband = 1.21 Hband = 1.64 Kband = 2.04 bands = [Uband,Bband,Vband,Rband,Iband,Jband,Hband] nbands = n_elements(bands) jband = indgen(nbands)*0 for k=0,nbands-1 do begin xx=min(abs(st.polext.wav-bands(k)),j) jband(k) = j endfor ;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1) res = poly_fit(logx(jband),logy(jband),2) print print,"Serkowski Fit with free K over ",string(st.polext(jband).wav,format='(100F5.2)') print,"-------------------------" K = -res(2) lmax = exp(res(1)/(2*K)) pmax = exp(res(0)+K*alog(lmax)^2) print,"* lmax = ",lmax," micron" print,"* K = ",K print,"* pmax = ",pmax," for NH=1e21" print ; Plot if needed IF !dustem_show_plot NE 0 THEN BEGIN ; ; Polarization fraction in the UV-IR ; win+=1 ; ; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,50],/xsty,win=win;multi=0,ps=filename ; ; ; Polarisation curve (Serkowski) win+=1 ;xr=[0.3,1e4] xr=[1e-4,3.5] yr=[1e-30,5e-24] dustem_plot_polext,st,p_dim,dustem_polext,aligned=st_model.align.grains.aligned,win=win ; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,/donotclose,multi=1 ; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,multi=2 ; win+=1 ; dustem_plot_fpol,win=win,xr=[10,1000],yr=[0,1.1] ; ENDIF ;dustem_plot_polar,st,/print_ratio,cont=cont dustem_out = [ dustem_out, dustem_polext ] ENDIF END ;USE OF POLFRAC DEPRECATED--------------------------------------- ; 'POLFRAC': BEGIN ; ; ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_POLSED AND DUSTEM_COMPUTE_SED ; ; IF !run_lin then begin ; ; dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st.polsed,cont=cont,freefree=freefree,synchrotron=synchrotron) ; dustem_sed_tmp = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron) ; dustem_sed = dustem_polsed ; ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt) ; for i = 0, count_filt-1 do begin ; j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count) ; if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0)) ; endfor ; dustem_polfrac = dustem_polsed/dustem_sed ; ; ; We normalize at 353GHz ; ;x=min(((*!dustem_data.polfrac).wav-850)^2,j353) ; ;dustem_polfrac = dustem_polfrac/dustem_polfrac(j353)*(*!dustem_data.polfrac).values(j353) ; ; chi2_polfrac = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2) ; n_polfrac = n_elements((*!dustem_data.polfrac).values) ; rchi2_polfrac = chi2_polfrac / n_polfrac ; wrchi2_polfrac = rchi2_polfrac * !fit_rchi2_weight.polfrac ; print,"chi2 POLFRAC = ",chi2_polfrac," rchi2 POLFRAC = ",rchi2_polfrac;," weighted rchi2 POLFRAC=",wrchi2_polfrac ; ;print,"chi2 per wl: ",((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2 ; ; dustem_out = [ dustem_out, dustem_polfrac ] ; ; ; Plot if needed ; IF !dustem_show_plot NE 0 THEN BEGIN ; ; ; Also show prediction for polarized P/I ; win+=1 ; dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,yr=[0,0.25],/ysty,xr=[10,5d3];dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,yr=[0,0.25],/ysty,xr=[10,5d3] ; ; ENDIF ; ; ENDIF ; ; END ;----------------------------------------------------------------- 'POLSED': BEGIN ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED IF !run_lin THEN BEGIN dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron) chi2_polsed = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2) n_polsed = n_elements((*!dustem_data.polsed).values) rchi2_polsed = chi2_polsed / n_polsed wrchi2_polsed = rchi2_polsed * !fit_rchi2_weight.polsed message,"chi2 POLSED = "+strtrim(chi2_polsed,2)+" rchi2 POLSED = "+strtrim(rchi2_polsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2 ;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma dustem_out = [ dustem_out, dustem_polsed ] ; classic command ;----------------------------------------------- ; Plot if needed IF !dustem_show_plot NE 0 THEN BEGIN win+=1 ; xr=[10,8e3] yr=[1e-34,1e-28] dustem_plot_polsed,st,p_dim,dustem_polsed,aligned=st_model.align.grains.aligned,win=win ; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,/donotclose,multi=1 ; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,multi=2 ENDIF ENDIF END 'QSED': BEGIN toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec) chi2_qsed =total(((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2) n_qsed = n_elements((*!dustem_data.qsed).values) rchi2_qsed = chi2_qsed / n_qsed wrchi2_qsed = rchi2_qsed * !fit_rchi2_weight.qsed dustem_out = [ dustem_out, dustem_qsed ] ; classic command message,"chi2 QSED = "+strtrim(chi2_qsed,2)+" rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2 if !dustem_show_plot NE 0 then begin fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40 fact1 = 1.00E-20;(3.e10/(1.e-4*((*!dustem_data.qsed).wav)))/1.d40 ; for data plotting ; fact=1/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.00e17 ; fact1=1.E-20 ; win+=1 window ,win ,title='DUSTEM WRAP (QSED)' titstq='Stokes Q Polarization Intensity' ytitstq=textoidl('Q_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2') xtit=textoidl('\lambda (\mum)') xr=[10,2e4] normpol=Q_spec*fact normpol1=dustem_qsed*fact1 ;normpol1=(*!dustem_data.qsed).values*fact1 yr=[1e-28,1e-17] ;yr=[1e-34,1e-23] indq=where(Q_spec lt 0, countq) ;indq=where((*!dustem_data.qsed).values lt 0, countq) ylog=1 if countq ne 0 then begin ylog=0 yr=[-2e-31,5e-31] endif cgplot,st.polsed.wav,Q_spec*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)' cgoplot,st.polsed.wav,Q_spec*fact,color='black' cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,color='Tomato',psym=16,symsize=1,thick=2 err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1,color=cgColor('Tomato') cgoplot,(*!dustem_data.qsed).wav,dustem_qsed*fact1,psym=6,color='red',symsize=2 cgplot,st.polsed.wav,Q_spec*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1 cgoplot,st.polsed.wav,Q_spec*fact/normpol,color='black' cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,color='Tomato',psym=16,symsize=1,thick=2 err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1/normpol1,color=cgColor('Tomato') endif END 'USED': BEGIN toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec) chi2_used =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2) n_used = n_elements((*!dustem_data.used).values) rchi2_used = chi2_used / n_used wrchi2_used = rchi2_used * !fit_rchi2_weight.used dustem_out = [ dustem_out, dustem_used] ; classic command message,"chi2 USED = "+strtrim(chi2_used,2)+" rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2 if !dustem_show_plot ne 0 then begin fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40 fact1 =1.00E-20;(3.e10/(1.e-4*((*!dustem_data.used).wav)))/1.d40 ; for data plotting win+=1 window ,win ,title='DUSTEM WRAP (USED)' titstu='Stokes U Polarization Intensity' ytitstu=textoidl('U_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2') xtit=textoidl('\lambda (\mum)') xr=[10,6e4] normpol2=U_spec*fact normpol3=dustem_used*fact1 ;yr=[min(normpol2)-min(normpol2)/10,max(normpol2)+max(normpol2)*10] ;conditions here are on computed spectra and they should be on the data itself. ;This even impairs the visualization of the data. Correct this!!!!! yr=[1e-28,1e-17];[1e-34,1e-23] indu=where(U_spec lt 0, countu) ;indu=where((*!dustem_data.used).values lt 0, countu) ylog=1 if countu ne 0 then begin ylog=0 yr=[-2e-31,5e-31] endif cgplot,st.polsed.wav,U_spec*fact,xtit='',ytit=ytitstu,tit=titstu,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black' cgoplot,st.polsed.wav,U_spec*fact,color='black' cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,color='Teal',psym=16,symsize=1,thick=2 err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1,color=cgColor('Teal') cgoplot,(*!dustem_data.used).wav,dustem_used*fact1,psym=6,color='red',symsize=2 cgplot,st.polsed.wav,U_spec*fact/normpol2,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1 cgoplot,st.polsed.wav,U_spec*fact/normpol2,color='black' cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,color='Teal',psym=16,symsize=1,thick=2 err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1/normpol3,color=cgColor('Teal') ;cgwindow, 'plot', endif END 'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION chi2_qext =total(((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2) n_qext = n_elements((*!dustem_data.qext).values) rchi2_qext = chi2_qext / n_qext wrchi2_qext = rchi2_qext * !fit_rchi2_weight.qext dustem_out = [ dustem_out, dustem_qext] ; classic command message,"chi2 QEXT = "+strtrim(chi2_qext,2)+" rchi2 QEXT = "+strtrim(rchi2_qext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2 if !dustem_show_plot ne 0 then begin win+=1 window ,win ,title='DUSTEM WRAP (QEXT)' titqex='Q-polarized Cross-section' ytitqex=textoidl('\sigma_{Q}/N_H (cm^2/H)') xtit=textoidl('\lambda^{-1} (\mum^{-1})') ;xr=[0.3,1e4] xr=[1e-4,3.5] yr=[1e-34,5e-24] indqex=where(Q_ext lt 0, countqex) ylog=1 if countqex ne 0 then begin ylog=0 yr=[-9e-27,5e-27] endif normpol = st.polext.ext_tot * 0. + 1.d21 cgplot,1/st.polext.wav,Q_ext,xtit='',ytit=ytitqex,tit=titqex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black' cgoplot,1/st.polext.wav,Q_ext/normpol,color='black' cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,color='Indian Red',psym=16,symsize=1,thick=2 err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,yrms=3.*((*!dustem_data.qext).sigma)/2/normpol,color=cgColor('Indian Red') cgoplot,1/(*!dustem_data.qext).wav,dustem_qext/normpol,psym=6,color='red',symsize=2 cgplot,1/st.polsed.wav,Q_ext/Q_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1 cgoplot,1/st.polsed.wav,Q_ext/Q_ext,color='black' cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,color='Indian Red',psym=16,symsize=1,thick=2 err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,yrms=3.*((*!dustem_data.qext).sigma)/2.,color=cgColor('Indian Red') endif END ; 'UEXT': BEGIN chi2_uext =total(((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2) n_uext = n_elements((*!dustem_data.uext).values) rchi2_uext = chi2_uext / n_uext wrchi2_uext = rchi2_uext * !fit_rchi2_weight.uext dustem_out = [ dustem_out, dustem_uext] ; classic command message,"chi2 UEXT = "+strtrim(chi2_uext,2)+" rchi2 UEXT = "+strtrim(rchi2_uext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2 if !dustem_show_plot ne 0 then begin win+=1 window ,win ,title='DUSTEM WRAP (UEXT)' tituex='U-polarized Cross-section' ytituex=textoidl('\sigma_{U}/N_H (cm^2/H)') xtit=textoidl('\lambda^{-1} (\mum^{-1})') ;xr=[0.3,1e4] xr=[1e-4,3.5] yr=[1e-34,5e-24] induex=where(U_ext lt 0, countuex) ylog=1 if countuex ne 0 then begin ylog=0 yr=[-9e-27,5e-27] endif normpol = st.polext.ext_tot * 0. + 1.d21 cgplot,1/st.polext.wav,U_ext,xtit='',ytit=ytituex,tit=tituex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black' cgoplot,1/st.polext.wav,U_ext/normpol,color='black' cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,color='Steel Blue',psym=16,symsize=1,thick=2 err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue') cgoplot,1/(*!dustem_data.uext).wav,dustem_uext/normpol,psym=6,color='red',symsize=2 cgplot,1/st.polsed.wav,U_ext/U_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1 cgoplot,1/st.polsed.wav,U_ext/U_ext,color='black' cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,color='Steel Blue',psym=16,symsize=1,thick=2 err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue') endif END ELSE : IF tagnames(ind_data_tag(i_tag)) NE 'POLFRAC' THEN stop ENDCASE ENDFOR DOF = n_elements(p_min) total_chi2 = 0 total_rchi2 = 0 total_n = 0 if !fit_rchi2_weight.ext ne 0 then begin total_n += n_ext total_chi2 += chi2_ext total_rchi2 += rchi2_ext endif if !fit_rchi2_weight.sed ne 0 then begin total_n += n_sed total_chi2 += chi2_sed total_rchi2 += rchi2_sed endif if !run_pol then begin if !fit_rchi2_weight.polext ne 0 then begin total_n += n_polext total_chi2 += chi2_polext total_rchi2 += rchi2_polext endif if !fit_rchi2_weight.polsed ne 0 then begin total_n += n_polsed total_chi2 += chi2_polsed total_rchi2 += rchi2_polsed endif if !fit_rchi2_weight.polfrac ne 0 then begin total_n += n_polfrac total_chi2 += chi2_polfrac total_rchi2 += rchi2_polfrac endif if !fit_rchi2_weight.qsed ne 0 then begin total_n += n_qsed total_chi2 += chi2_qsed total_rchi2 += rchi2_qsed endif if !fit_rchi2_weight.used ne 0 then begin total_n += n_used total_chi2 += chi2_used total_rchi2 += rchi2_used endif if !fit_rchi2_weight.qext ne 0 then begin total_n += n_qext total_chi2 += chi2_qext total_rchi2 += rchi2_qext endif if !fit_rchi2_weight.uext ne 0 then begin total_n += n_uext total_chi2 += chi2_uext total_rchi2 += rchi2_uext endif endif total_wchi2 = chi2_ext * !fit_rchi2_weight.ext $ + chi2_sed * !fit_rchi2_weight.sed if !run_pol then begin total_wchi2 += chi2_polext * !fit_rchi2_weight.polext $ + chi2_polsed * !fit_rchi2_weight.polsed $ + chi2_polfrac * !fit_rchi2_weight.polfrac $ + chi2_qsed * !fit_rchi2_weight.qsed $ + chi2_used * !fit_rchi2_weight.used $ + chi2_qext * !fit_rchi2_weight.qext $ + chi2_uext * !fit_rchi2_weight.uext endif message,'======================================================================',/continue print,"Total weighted chi2 driving mpfitfun = ", total_wchi2 print,"Nb of data points = ", total_n print,"DOF = ", DOF print,"Reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",rchi2_ext,rchi2_sed,rchi2_polext,rchi2_polsed,rchi2_polfrac print,"Weighted reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",wrchi2_ext,wrchi2_sed,wrchi2_polext,wrchi2_polsed,wrchi2_polfrac if !run_pol then print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed,!fit_rchi2_weight.polext,!fit_rchi2_weight.polsed,!fit_rchi2_weight.polfrac $ else print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF) print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF) print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF);total_rchi2/4*total_n/(total_n-DOF) original message,'======================================================================',/continue ;print,"Total reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac)/(n_sed+n_ext+n_polext+n_polsed+n_polfrac-DOF) ; Reduced chi2 is like if there were the same nb of data points for each curve ;print,"Mean reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac) * total_n/(total_n-DOF) ; Plot size distribution win+=1 xr=[0.3,5e3] yr=[0.1,100] ;IF !dustem_show_plot NE 0 THEN dustem_plot_sdist,xr=xr,yr=yr,ps=filename,win=win;,/donotclose ; Suppress the artificial [0] at the beginning dustem_out = dustem_out[1:*] ;IF defined(extra) THEN BEGIN ; IF keyword_set(extra.plot_it) THEN BEGIN ;; message,'params='+strtrim(p_min,2),/info,/continue ;; stop ; print,'dustem_mpfit_run_vg: Current parameters values:',p_dim ; ENDIF ;ENDIF the_end: RETURN,dustem_out END