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]) ; MODIFICATION 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. ;COMMON MPFIT_ERROR, ERROR_CODE IF keyword_set(help) THEN BEGIN doc_library,'dustem_mpfit_run_vg' sed=0. goto,the_end ENDIF p_dim=p_min*(*(*!dustem_fit).param_init_values) ;Run dustem with as an interface to mpfitfun IF !dustem_idl_continuum THEN cont=dustem_create_continuum() IF !dustem_idl_freefree THEN freefree=dustem_create_freefree() IF !dustem_idl_synchrotron THEN synchrotron=dustem_create_synchrotron() ;CALL THE DUSTEM_CREATE FUNCTIONS ;Actually, res is never used ;THIS line below was added by NF, but removed by JPB on Dec 5 2009 ;dustem_set_params,p_dim ;NF on October 2009 (params have to be updated before calling some functions !) ;dustem_create_functions,p_dim/(*!IDO),cont=cont,res=res dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),cont=cont,freefree=freefree,synchrotron=synchrotron,res=res ;if ISRF has been modified and IONFRAC has not been modified, then call dustem_create_ionfrac ;NF on October 2009 : handling of DUSTEM_CREATE_IONFRAC has been changed ;this call only if IONFRACPAH master-keyword in GRAIN.DAT ;Removed for now (JPB). Use of ionfrac should be fully revisited. IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac() ;RUN THE MODEL AND GET THE SPECTRUM 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 rchi2_sed = 0 rchi2_ext = 0 rchi2_polext = 0 rchi2_polsed = 0 rchi2_polfrac = 0 wrchi2_sed = 0 wrchi2_ext = 0 wrchi2_polext = 0 wrchi2_polsed = 0 wrchi2_polfrac = 0 n_sed = 0 n_ext = 0 n_polext = 0 n_polsed = 0 n_polfrac = 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 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,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_plot_nuinu_em,st,dustem_sed,cont,/print_radiance,p_dim=p_dim ; Plot if needed IF !dustem_show_plot NE 0 THEN BEGIN window,0 dustem_plot_fit_sed,st,dustem_sed,cont,freefree,synchrotron,_extra=extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2 ENDIF dustem_out = [ dustem_out, dustem_sed ] END 'EXT' : BEGIN 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 dustem_out = [ dustem_out, dustem_ext ] END 'POLEXT': BEGIN 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 dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav) 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) ;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1) ; 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,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,4],/xsty,win=win ; Polarisation curve (Serkowski) win+=1 xr=[0.1,50] yr=[5e-25,5e-23] dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,/donotclose,multi=1 dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,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 'POLFRAC': BEGIN IF !run_lin then begin dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont,freefree=freefree,synchrotron=synchrotron) dustem_sed_tmp = 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,cont=cont,freefree=freefree,synchrotron=synchrotron,yr=[0,0.25],/ysty,xr=[10,5d3] ENDIF ENDIF END 'POLSED': BEGIN IF !run_lin THEN BEGIN dustem_polsed = 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 ] ; Plot if needed IF !dustem_show_plot NE 0 THEN BEGIN win+=1 xr=[10,5e3] yr=[1e-34,1e-28] dustem_plot_polar,st,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,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 ELSE : 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 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 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) 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