PRO dustem_fit_spectro_example,model=model $ ,sed_file=sed_file $ ,Nitermax=Nitermax $ ,fits_save_and_restore=fits_save_and_restore $ ,ucfactor=ucfactor $ ,wait=wait $ ,noobj=noobj $ ,help=help $ ,verbose=verbose ;+ ; NAME: ; dustem_fit_spectro_example ; ; PURPOSE: This routine is an example of how to fit an observational ; spectrometer data in emission (StokesI) with DustEM and ; DustEMWrap. ; ; See the DustEMWrap User Guide for more information. ; ; CATEGORY: ; DustEMWrap, Distributed, High-Level, User Example ; ; CALLING SEQUENCE: ; dustem_fit_spectro_example[,model=][sed_file=][ext_file=][,Nitermax=][,fits_save_and_restore=][,/help,/wait,/verbose,/noobj] ; ; 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_2.xcat' is used. ; Nitermax = maximum number of fit iterations. Default is 5. ; ucfactor = Multiplicative factor for uncertainties on ; spectrometer data. Default is 1. ; 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 ; noobj = if set, runs with no object graphics ; ; 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: ; ; EXAMPLES ; dustem_fit_spectro_example ; dustem_fit_spectro_example,Nitermax=10,fits_save_and_restore='/tmp/mysavefile.fits' ; dustem_fit_spectro_example,model='DBP90' ; ; MODIFICATION HISTORY: ; Written by AH Oct 2022 ; 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_spectro_example' goto,the_end END IF keyword_set(model) THEN BEGIN use_model=strupcase(model) ENDIF ELSE BEGIN use_model='WD01_RV5P5B' 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. ;default is no polarization in models 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_ucfactor=1. if keyword_set(ucfactor) then use_ucfactor=ucfactor 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 ;=== Refer to the DustEM and DustEMWrap userguides for an explanation ; of the different grain types ;=== 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. ;=== Example is provided for WD01_RV5P5B ;=== If you choose a different dust model, you will need to adjust the ;=== parameter information accordingly. pd = [ '(*!dustem_params).G0', $ ; multiplicative factor of the radiation field '(*!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 = [2., 5.4e-4, 5.4e-4, 1.8e-4, 1.65e-3, 7.e-3] ; starting guesses Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(1.e-15,Npar) ;=== Fixed parameters fpd=[] ;=== Initial parameter values for fixed parameters fiv=[] 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 IF keyword_set(noobj) THEN !dustem_noobj=1 !EXCEPT=2 ; for debugging ;=== READ EXAMPLE EMISSION DATA TO FIT dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/' file=dir+'example_SED_2.xcat' if keyword_set(sed_file) then file=sed_file sed=read_xcat(file,/silent) ;=== TWEAK THE UNCERTAINTIES ON THE SPECTROMETER DATA TO ADJUST THEIR ;=== RELATIVE WEIGHT IN THE FIT fidx=where(sed.filter eq 'SPECTRUM',fct) if fct gt 0 then sed[fidx].sigmaii = use_ucfactor*sed[fidx].sigmaii ;== SET THE OBSERVATIONAL STRUCTURE dustem_set_data, sed,sed ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT ;== AND ACTIVATE ANY PLUGINS dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims,polarization=use_polarization 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-10 ;fit tolerance ;=== INFORMATION TO MAKE THE PLOT ;=== _x/_extinction means extinction plots, _m/_emissions means emission plots xr_m = [1.,5e5] yr_m = [5e-8,1.00e6] tit_m='Spectral Energy Distribution' ;=== RUN THE FIT t1=systime(0,/sec) res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $ ,/xlog,/ylog,xr_m=xr_m,yr_m=yr_m,xr_x=xr_x,yr_x=yr_x,xtit=xtit,ytit=ytit,tit_m=tit_m $ ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $ ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot) 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 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 dustem_read_fits_table,filename=fits_save_and_restore,dustem_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,tit_m=tit_m+' (From Saved FITS file)' ENDIF ELSE BEGIN dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (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 ENDIF ENDIF message,'Finished dustem_fit_spectro_example',/info the_end: END