From 906f3c27c86fd44cb347e92baeef51ef0eac525c Mon Sep 17 00:00:00 2001 From: Jean-Philippe Bernard Date: Fri, 7 May 2021 23:00:56 +0200 Subject: [PATCH] First commit --- src/idl/dustem_fit_spin_test_readme.pro | 202 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 202 insertions(+), 0 deletions(-) create mode 100644 src/idl/dustem_fit_spin_test_readme.pro diff --git a/src/idl/dustem_fit_spin_test_readme.pro b/src/idl/dustem_fit_spin_test_readme.pro new file mode 100644 index 0000000..0a030ba --- /dev/null +++ b/src/idl/dustem_fit_spin_test_readme.pro @@ -0,0 +1,202 @@ +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_modified_bb_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_modified_bb_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=['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 + +;=== 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_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] + +;=== 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,90.,500.,600.,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.] + +;== 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-spin' +(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-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.) +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 + +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=[1.,2.e4],/xstyle,yrange=[1e-4,1.e3],/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=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav' +file_out=!dustem_res+'DUSTEM_bb_fit_example.sav' +dustem_save_sed_fit,file_out + +;=== Plot best fit +yr=[1e-3,100] +xr=[50,4000] +tit='DUSTEM Modified BB Example' +ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') +xtit=textoidl('\lambda (\mum)') +;errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values) +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=yr,xr=xr,/ysty,/xsty, $ + res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit' + +;stop + +;====================================== +;====You can exit IDL here and re-enter +;====================================== +window,1 + +;=== RESTORE FIT RESULTS +;=== initialise dustem +dustem_init +file=!dustem_res+'DUSTEM_bb_fit_example.sav' +dustem_restore_sed_fit,file + +;=== recover best fit values +res=*(*!dustem_fit).current_param_values +chi2=(*!dustem_fit).chi2 +rchi2=(*!dustem_fit).rchi2 +errors=*(*!dustem_fit).current_param_errors + +;=== Plot best fit +yr=[1e-3,100] +xr=[50,4000] +tit='DUSTEM Modified BB 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_bb_fit.ps' + device,filename=ps_file,/color +ENDIF +dustem_sed_plot,*(*!dustem_fit).current_param_values, $ + ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit' +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 -- libgit2 0.21.2