PRO dustem_fit_intensity_example,model=model $ ,sed_file=sed_file $ ,Nitermax=Nitermax $ ,postscript=postscript $ ,fits_save_and_restore=fits_save_and_restore $ ,help=help $ ,wait=wait $ ,verbose=verbose ;+ ; NAME: ; dustem_fit_intensity_example ; ; PURPOSE: ; This routine is an example of how to fit an observational SED ; (StokesI only) with DustEM and DustEMWrap. The objective is to ; illustrate how to use DustEMWrap and not to do science -- the fit ; obtained by running this example is likely to be poor. ; ; For this example, the code uses the SED in the file example_SED_1.xcat, ; which is distributed in the Data/EXAMPLE_OBSDATA/ directory ; ; The example SED has Stokes I photometric data points from ; IRAC, MIPS and IRAS. Examples illustrating running DustEMWrap to ; fit spectral data, polarisation data and extinction data ; are provided in other _example routines in the src/idl/ ; directory. See the DustEMWrap User Guide for more information. ; ; CATEGORY: ; DustEMWrap, Distributed, High-Level, User Example ; ; CALLING SEQUENCE: ; dustem_fit_intensity_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save_and_restore=][,/help,/wait,/verbose] ; ; INPUTS: ; None ; ; OPTIONAL INPUT PARAMETERS: ; None ; ; OUTPUTS: ; None ; ; OPTIONAL OUTPUT PARAMETERS: ; Plots, results structure in binary FITS table format ; ; ACCEPTED KEY-WORDS: ; model = specifies the interstellar dust mixture used by DustEM ; 'MC10' model from Compiegne et al 2010 ; 'DBP90' model from Desert et al 1990 ; 'DL01' model from Draine & Li 2001 ; 'WD01_RV5p5B' model from Weingartner & Draine 2002 with Rv=5.5 ; 'DL07' model from Draine & Li 2007 ; 'J13' model from Jones et al 2013, as updated in ; Koehler et al 2014 ; 'G17_ModelA' model A from Guillet et al (2018). Includes ; polarisation. See Tables 2 and 3 of that paper for details. ; 'G17_ModelB' model B from Guillet et al (2018) ; 'G17_ModelC' model C from Guillet et al (2018) ; 'G17_ModelD' model A from Guillet et al (2018) ; sed_file = string naming the path to text file in .xcat format that ; describes the observational SED. If not set, the file ; 'Data/EXAMPLE_OBSDATA/example_SED_1.xcat' is used. ; postscript = if set, final plot is saved as postscript in the ; current working directory ; Nitermax = maximum number of fit iterations. Default is 5. ; fits_save_and_restore = if set, save the fit results in a binary ; FITS file. The code then restore this file and plots ; the results using the saved results information. ; help = if set, print this help ; wait = if set, wait this many seconds between each step of ; the code (for illustration purposes) ; verbose = if set, subroutines will run in verbose mode ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; None ; ; RESTRICTIONS: ; The DustEM fortran code must be installed ; The DustEMWrap IDL code must be installed ; ; PROCEDURES AND SUBROUTINES USED: ; *** COMMENT AH --> is this really NONE? **** ; ; EXAMPLES ; dustem_fit_intensity_example ; dustem_fit_intensity_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits' ; dustem_fit_intensity_example,model='DBP90' ; ; MODIFICATION HISTORY: ; Written by JPB Apr-2011 ; 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_fit_intensity_example' goto,the_end END IF keyword_set(model) THEN BEGIN use_model=strupcase(model) ENDIF ELSE BEGIN use_model='DBP90' ;Example with default keywords uses the DBP90 model ENDELSE known_mdls=['MC10','DBP90','DL01','WD01_RV5P5B','DL07','J13','G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD'] test_model = where(known_mdls eq use_model,ct) if ct eq 0 then begin message,'ISM dust model '+use_model+' unknown',/continue message,'Known models are MC10,DBP90,DL01,WD01_RV5P5B,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue stop end use_polarization=0 ; initialize Dustemwrap in no polarization mode use_window=2 ; default graphics window number to use for plotting the results use_verbose=0 if keyword_set(verbose) then use_verbose=1 use_Nitermax=5 ; maximum number of iterations for the fit IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax dustem_define_la_common ;=== Set the (model-dependent) parameters that you want to fit (pd), ;=== their initial values (iv) ;=== and whether they are bounded (ulimed,llimed,llims,ulims). ;=== Fixed parameters (fpd) and their values (fiv) are also set here. ;=== Refer to the DustEM and DustEMWrap User guides for an explanation ;=== of the physical meaning of dust model and plug-in parameters, and ;=== how to specify them. ;=== Examples are provided for some of the dust models. ;=== To try them, uncomment the model that you want to try and re-run ;=== AN EXAMPLE FOR DBP90 ;=== Here we fit the dust abundances of the DBP90 model, the ;=== intensity of the dust-heating radiation field, as well as several plug-ins: ;=== (i) continuum due to a blackbody ;=== (ii) free-free emission ;=== (iii) synchrotron emission ;=== The free parameters are all lower-bounded at zero. ;=== use_model='DBP90' ; you should specify this above! pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;VSG mass fraction '(*!dustem_params).grains(2).mdust_o_mh'];,$ ;BG mass fraction ; 'dustem_plugin_synchrotron_1', $ ;Spectral index of CREs ; 'dustem_plugin_synchrotron_2'] ;Synchrotron amplitude at 10 mm ; initial values vector for run without synchrotron plugin ; [G0,PAH0,VSG,BG] iv = [1.6, 2.2e-4, 5.7e-4, 3.4e-3] ; initial values vector for run including the synchrotron plugin ; [G0,PAH0,VSG,BG,alpha_CR,Amp_syn] ;iv = [1.6, 2.2e-4, 5.7e-4, 3.4e-3, 2.7,0.01] ; initial values vector for run including the synchrotron plugin and ; fixing the G0 and BB continuum parameters ;iv = [2.2e-4, 5.7e-4, 3.4e-3, 2.7,0.01] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) fpd=[] & fiv=[] ;Example using fixed parameters (ISRF and a NIR continuum) ; Uncomment the following lines to include these fixed parameters ;fpd = ['(*!dustem_params).G0', $ ;G0 ; 'dustem_plugin_continuum_1', $ ;Temperature of a BB ; 'dustem_plugin_continuum_2'] ;Peak amplitude of a BB ; fixed values vector ; [G0,T_BB,Amp_BB] ;fiv=[0.5,750,1.e-2] ;; ;=== AN EXAMPLE FOR DL01 (or DL07 or WD01_RV5P5B) ;; ;=== Here we fit the dust abundances of the model. The ;; ;=== intensity of the dust-heating radiation field is fixed to 1.2*G0. ;; ;=== The free parameters are all lower-bounded at zero. ;; use_model='DL01' ; you should specify this above! ;; pd = [ $ ;; '(*!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', $ ;Gra ;; '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra ;; '(*!dustem_params).grains(4).mdust_o_mh'] ;aSil ;; iv = [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3] ;; Npar=n_elements(pd) ;; ulimed=replicate(0,Npar) ;; llimed=replicate(1,Npar) ;; llims=replicate(0.,Npar) ;; fpd=['(*!dustem_params).G0'] ;; fiv=[1.2] ;; ;=== AN EXAMPLE FOR MC10 ;; ;=== Here we fit the dust abundances of the MC10 model, the ;; ;=== intensity of the dust-heating radiation field as well as a plug-in: ;; ;=== (i) continuum due to a blackbody ;; ;=== The intensity of the dust-heating radiation field is fixed to ;; ;=== 1.5*G0 and the tmperature of the blackbody is fixed to 1200K ;; ;=== The free parameters in the fit are lower-bounded at zero. ;; use_model='MC10' ; you should specify this above! ;; pd = [ $ ;; '(*!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_plugin_continuum_2' $ ;Intensity at peak of the continuum ;; ] ;; iv = [ 7.8e-4, 7.8e-4, 1.65e-4, 1.45e-3, 6.7e-3, 0.003] ;; Npar=n_elements(pd) ;; ulimed=replicate(0,Npar) ;; llimed=replicate(1,Npar) ;; llims=replicate(0.,Npar) ;; fpd=[ $ ;; '(*!dustem_params).G0' $ ; ISRF intensity (multiplicative factor x G0) ;; ,'dustem_plugin_continuum_1' $ ;temperature of blackbody the produces the continuum ;; ] ;; fiv=[1.5, 1200.] ;; ;=== AN EXAMPLE FOR J13 ;; ;=== Here we fit the dust abundances of the MC10 model, the ;; ;=== intensity of the dust-heating radiation field as well as two plug-ins: ;; ;=== (i) free-free emission ;; ;=== (ii)continuum due to a blackbody ;; ;=== The temperature of the blackbody is fixed to 1000K ;; ;=== The free parameters in the fit are all lower-bounded at zero. ;; use_model='J13' ; you should specify this ! ;; pd = [ $ ;; '(*!dustem_params).G0' $ ;G0 ;; ,'(*!dustem_params).grains(0).mdust_o_mh'$ ;CM20 -- power law size distribution ;; ,'(*!dustem_params).grains(1).mdust_o_mh'$ ;CM20 -- logN size distribution ;; ,'(*!dustem_params).grains(2).mdust_o_mh' $ ;aPyM5 ;; ,'(*!dustem_params).grains(3).mdust_o_mh' $ ;aOlM5 ;; ,'dustem_plugin_freefree_1' $ ;ionized gas temperature ;; ,'dustem_plugin_freefree_2' $ ;free-free amplitude ;; ,'dustem_plugin_continuum_2'] ;Intensity at peak of the BB continuum ;; iv = [1.2,1.7e-3, 6.3e-4, 2.55e-3, 2.55e-3, 7500., 0.4, 0.001] ;; Npar=n_elements(pd) ;; ulimed=replicate(0,Npar) ;; llimed=replicate(1,Npar) ;; llims=replicate(0.,Npar) ;; fpd=[ $ ;; 'dustem_plugin_continuum_2'] ; Temperature of the BB ;; fiv=[1000.] ;;Nfix=n_elements(fpd) ;;if n_elements(fiv) ne Nfix then begin ;; message,'Number of fixed parameters (fpd) does not equal number of initial values of fixed parameters (fiv)',/info ;; stop ;;end if keyword_set(wait) then begin message,'Finished setting dust model and plug-in parameters: '+use_model,/info wait,wait end ;== INITIALISE DUSTEM dustem_init,model=use_model,polarization=use_polarization !dustem_nocatch=1 !dustem_verbose=use_verbose !dustem_show_plot=1 ;=== READ EXAMPLE SED DATA dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/' file=dir+'example_SED_1.xcat' IF keyword_set(sed_file) THEN file=sed_file sed=read_xcat(file,/silent) if keyword_set(wait) then begin message,'Finished reading SED data: '+file,/info wait,wait end ;;=== ADJUST THE UNCERTAINTIES FROM WITHIN THE CODE (FOR ILLUSTRATION) ind=where(sed.sigmaII LT (0.2*sed.StokesI)^2,count) IF count NE 0 THEN sed[ind].sigmaII=(0.2*sed[ind].StokesI)^2 ;== SET THE OBSERVATIONAL STRUCTURE ;== sed is passed twice -- the first occurrence is the SED that you ;== wish to fit, the second occurrence is the SED that you wish to visualise. dustem_set_data,sed,sed ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE ;== ADJUSTED DURING THE FIT ;dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims ;== INITIALIZE ANY PLUGINS ;dustem_init_plugins,pd,fpd ;== INITIALIZE ANY FIXED PARAMETERS FOR PLUGINS ;if Nfix gt 0 then dustem_init_fixed_params,fpd,fiv ;== INITIALIZE ANY FREE AND FIXED PARAMETERS, PLUGINS dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims if keyword_set(wait) then begin message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info wait,wait end ;=== INFORMATION TO RUN THE FIT tol=1.e-16 ;fit tolerence ;=== INFORMATION TO MAKE THE PLOT yr=[1.00e-4,1.00E2] ; y-axis limits xr=[1.00E0,6.00e4] ; x-axis limits tit='FIT INTENSITY EXAMPLE' ; plot title ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title xtit=textoidl('\lambda (\mum)') ; x-axis title ;=== RUN THE FIT t1=systime(0,/sec) res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $ ,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit $ ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $ ,errors=errors,chi2=chi2,rchi2=rchi2) t2=systime(0,/sec) if keyword_set(wait) then begin message,'Finished running DustEMWrap, using Niters: '+strtrim(string(use_Nitermax),2),/info message,'Time taken [sec]: '+sigfig(t2-t1,2,/sci),/info wait,wait end ;IC: - ;The final plot is automatically plotted ;and a (Final run) string is automatically added ;to modify this go to dustemwrap_plot and search for the string '(FINAL_RUN)' ;This is true for OOP plotting. ;However, even thought the last iteration for noobj plotting is plotted, the '(FINAL_RUN)' string isn't applied ;to modify this, got to dustemwrap_plot_noobj and search for the '(RUNNING)' string. ;an extra condition using !dustem_end (signaling end of fit) will have to be used. ;=== MAKE THE FINAL PLOT IF keyword_set(postscript) THEN BEGIN dir_ps='./' set_plot,'PS' ps_file=dir_ps+postscript device,filename=ps_file,/color ENDIF IF !dustem_noobj THEN BEGIN dustemwrap_plot_noobj,*(*!dustem_fit).CURRENT_PARAM_VALUES,dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)' ENDIF ELSE BEGIN dustemwrap_plot,*(*!dustem_fit).CURRENT_PARAM_VALUES,dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)' ENDELSE IF keyword_set(postscript) THEN BEGIN set_plot,'X' device,/close message,'Wrote '+ps_file,/info ENDIF if keyword_set(wait) then begin message,'Made the plot of the final results',/info wait,wait end ;IC: Last iteration is automatically plotted. ;I think modifications will have to be applied to dustemwrap_plot... ;idea: create a system variable that signals if the fits file has been successfully saved and use it in the dustemwrap_plot(obj/noobj) routine(s) ;to add the '(From saved FITS file)' string to the title ; IF keyword_set(fits_save_and_restore) THEN BEGIN message,'Writing out results structure: '+fits_save_and_restore,/info dustem_write_fits_table,filename=fits_save_and_restore,help=help ;=== At this point, you could erase all dustem system variables, or exit idl... all the ;=== information needed to recover the results and remake the plots has been saved in the FITS table dustem_read_fits_table,filename=fits_save_and_restore,dustem_spectra_st=dustem_spectra_st ;==== plot result taken from the saved fits table res=*(*!dustem_fit).CURRENT_PARAM_VALUES IF !dustem_noobj THEN BEGIN dustemwrap_plot_noobj,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)' ENDIF ELSE BEGIN dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)' ENDELSE if keyword_set(wait) then begin message,'Saved the results as FITS in the file: '+fits_save_and_restore+', and made a plot using the data in this file',/info wait,wait end ENDIF message,'Finished dustem_fit_intensity_example',/info the_end: END