PRO dustem_fit_spin_test_readme,postcript=postcript,mode=mode,help=help ;This Readme describes how to fit SEDs with a modified Black-Body ;It runs on the SED stored into the file sample_SED.xcat ;Remember that the goal here is just to illustrate the method. ;Note that its is not necessary to use the .xcat file format, and data SED can be provided ;manually, but the observation structure must have the structure shown below. ;Obviously, the dustem package must have been installed succesfully ;see DustemWrap documentation for install instructions. ;+ ; NAME: ; dustem_fit_modified_bb_readme ; PURPOSE: ; This is an example of how to fit SEDs with a modified Black-Body ; using the dustem wrapper. ; It is meant to be an example to follow when writing your own ; programs using the dustem IDL wrapper. ; The SED used here is in sample_SED.xcat ; CATEGORY: ; Dustem ; CALLING SEQUENCE: ; dustem_fit_spin_test_readme,postcript=postcript,help=help ; INPUTS: ; None ; OPTIONAL INPUT PARAMETERS: ; None ; OUTPUTS: ; None ; OPTIONAL OUTPUT PARAMETERS: ; None ; ACCEPTED KEY-WORDS: ; postcript = if set plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps ; help = If set print this help ; COMMON BLOCKS: ; None ; SIDE EFFECTS: ; None ; RESTRICTIONS: ; The dustem idl wrapper must be installed ; PROCEDURE: ; None ; EXAMPLES ; dustem_fit_spin_test_readme ; MODIFICATION HISTORY: ; Written by J.P. Bernard June 9th 2011 ; see evolution details on the dustem cvs maintained at CESR ; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems. ;- IF keyword_set(help) THEN BEGIN doc_library,'dustem_fit_spin_test_readme' goto,the_end ENDIF ;=== initialise dustem dustem_init,mode='COMPIEGNE_ETAL2010' !dustem_verbose=1 !dustem_show_plot=1 ;=== make a fake one out for given filters filters=['DIRBE1','DIRBE2','DIRBE3','DIRBE4','IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3'] Nfilt=n_elements(filters) sed=dustem_initialize_sed(Nfilt) sed.filter=filters sed.wave=dustem_filter2wav(filters) sed.instru=dustem_filter2instru(filters) st=dustem_set_data(sed=sed) ;1.000000 ;PAH0_MC10 25 logn-chrg-spin-zm 7.8000E-04 2.2400E+00 3.5000E-08 1.2000E-07 6.4000E-08 1.0000E-01 ;PAH1_MC10 25 logn-chrg-spin-zm 7.8000E-04 2.2400E+00 3.5000E-08 1.2000E-07 6.4000E-08 1.0000E-01 ;amCBEx 15 logn 1.6500E-04 1.8100E+00 6.0000E-08 2.0000E-06 2.0000E-07 3.5000E-01 ;amCBEx 25 plaw-ed 1.4500E-03 1.8100E+00 4.0000E-07 2.0000E-04 -2.8000E+00 1.5000E-05 1.5000E-05 2.0000E+00 ;aSilx 25 plaw-ed 6.7000E-03 3.0000E+00 4.0000E-07 2.0000E-04 -3.4000E+00 2.0000E-05 2.0000E-05 2.0000E+00 goto,freefree_tag ;=== Set which parameters you want to fit pd = [ $ ; '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSilx '(*!dustem_params).gas.G0', $ ;Gas G0 '(*!dustem_params).gas.Tgas', $ ;Gas temperature '(*!dustem_params).gas.nH', $ ;Gas temperature ; 'dustem_create_freefree_1', $ ;Gas temperature for free-free emission ; 'dustem_create_freefree_2', $ ;H-alpha for free-free emission (R) 'dustem_create_continuum_2'] ;Intensity of NIR continuum ;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001] ;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,0.01,0.1] p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1] ;=== set initial value of parameters you want to fit ;iv=p_truth ;exact solution ;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,0.01,0.001] ;shifted iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001] ;shifted ulimed=[0 , 0 ,0,0,0,0,0,0,0] llimed=[1 , 1 ,1,1,1,1,1,1,1] llims =[0. ,0. ,0.,0.,0.,0.,0.,0.,0.] goto,out_of_there freefree_tag: ;=== Set which parameters you want to fit pd = [ $ ; '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSilx '(*!dustem_params).gas.G0', $ ;Gas G0 '(*!dustem_params).gas.Tgas', $ ;Gas temperature '(*!dustem_params).gas.nH', $ ;Gas temperature 'dustem_create_freefree_1', $ ;Gas temperature for free-free emission 'dustem_create_freefree_2', $ ;H-alpha for free-free emission (R) 'dustem_create_continuum_2'] ;Intensity of NIR continuum ;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001] p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,1.e-11,0.1] ;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1] ;=== set initial value of parameters you want to fit ;iv=p_truth ;exact solution iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,1.e-11,0.001] ;shifted ;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001] ;shifted ulimed=[0 , 0 ,0,0,0,0,0,0,0,0,0] llimed=[1 , 1 ,1,1,1,1,1,1,1,1,1] llims =[0. ,0. ,0.,0.,0.,0.,0.,0.,0.,0,0] synchrotron_tag: ;=== Set which parameters you want to fit pd = [ $ ; '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSilx '(*!dustem_params).gas.G0', $ ;Gas G0 '(*!dustem_params).gas.Tgas', $ ;Gas temperature '(*!dustem_params).gas.nH', $ ;Gas temperature 'dustem_create_synchrotron_1', $ ;CR spectral index for synchrotron emission 'dustem_create_synchrotron_2', $ ;Amplitude for synchrotron emission (R) 'dustem_create_continuum_2'] ;Intensity of NIR continuum ;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001] p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,3.,0.01,0.1] ;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1] ;=== set initial value of parameters you want to fit ;iv=p_truth ;exact solution iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,3.,0.01,0.1] ;shifted ;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001] ;shifted ulimed=[0 , 0 ,0,0,0,0,0,0,0,0,0] llimed=[1 , 1 ,1,1,1,1,1,1,1,1,1] llims =[0. ,0. ,0.,0.,0.,0.,0.,0.,0.,0,0] out_of_there: ;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo) dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims (*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin' (*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin' help,*!dustem_fit ;stop sed.spec=dustem_compute_sed(p_truth,sst,_extra=extra) sed.error=(sed.spec*0.1);>0.1 ;sed.error[*]=100. ;=== SET THE OBSERVATION STRUCTURE ;ind=where(sed.wave GE 5.) ind=where(sed.wave GE 0.1) st=dustem_set_data(sed=sed[ind]) ;=== RUN fit ;tol=1.e-20 tol=1.e-10 xtol=1.e-10 Nitermax=5 ;maximum number of iteration. This is the criterium which will stop the fit procedure xrange=[1.,2.e4] yrange=[1e-4,1.e3] loadct,13 ;!y.range=[1e-8,100] ;This is to ajust plot range from outside the routine t1=systime(0,/sec) res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []') t2=systime(0,/sec) print,(res-p_truth)/p_truth*100. ;stop ;=== SAVE FIT RESULTS file_out='/tmp/DUSTEM_spinning_fit_example.sav' dustem_save_system_variables,file_out message,'Saved '+file_out,/continue stop errors=(*(*!dustem_fit).current_param_errors) chi2=(*!dustem_fit).chi2 rchi2=(*!dustem_fit).rchi2 loadct,13 dustem_sed_plot,(*(*!dustem_fit).current_param_values), $ ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty, $ res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog ;stop ;====================================== ;====You can exit IDL here and re-enter ;====================================== file='/tmp/DUSTEM_spinning_fit_example.sav' dustem_restore_system_variables,file window,1 xrange=[1.,2.e4] yrange=[1e-4,1.e3] ;=== recover best fit values res=*(*!dustem_fit).current_param_values chi2=(*!dustem_fit).chi2 rchi2=(*!dustem_fit).rchi2 ;errors=*(*!dustem_fit).current_param_errors errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values) ;=== Plot best fit tit='DUSTEM spinning dust Example (Saved)' ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') xtit=textoidl('\lambda (\mum)') loadct,13 IF keyword_set(postcript) THEN BEGIN set_plot,'PS' ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_spinning_fit.ps' device,filename=ps_file,/color ENDIF dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog ;dustem_sed_plot,*(*!dustem_fit).current_param_values, $ ; ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog IF keyword_set(postcript) THEN BEGIN device,/close set_plot,'X' message,'wrote '+ps_file,/info ENDIF print,'dustem_mpfit_sed executed in ',t2-t1,' sec' the_end: END