PRO dustem_fit_intensity_mbb_example,model=model $ ,sed_file=sed_file $ ,Nitermax=Nitermax $ ,postscript=postscript $ ,fits_save=fits_save $ ,noobj=noobj $ ,help=help $ ,wait=wait $ ,verbose=verbose ;+ ; NAME: ; dustem_fit_intensity_mbb_example ; ; PURPOSE: ; This routine is an example of how to fit an observational SED ; (StokesI only) with DustEM and DustEMWrap. ; ; For this example, the code uses the SED in the file example_SED_1.xcat, ; which is distributed in the Data/EXAMPLE_OBSDATA/ directory ; ; In addition, this example uses the plugin dustem_plugin_mbbdy.pro to ; fit the FIR part of the spectrum using a modified black-body rather ; than a physical ISM dust model for the big grain population. ; ; 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_mbb_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save=][,/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. See userguide or dustem_test_model_exists.pro ; for more details about available models in current release. ; 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 ; Nitermax = maximum number of fit iterations. Default is 5. ; fits_save = if set, save the fit results in a binary ; FITS file. ; 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_intensity_mbb_example ; dustem_fit_intensity_mbb_example,Nitermax=1,fits_save='/tmp/mysavefile.fits' ; dustem_fit_intensity_mbb_example,model='DBP90' ; ; MODIFICATION HISTORY: ; Written by AH Apr-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_intensity_mbb_example' goto,the_end END IF keyword_set(model) THEN BEGIN use_model=strupcase(model) ENDIF ELSE BEGIN use_model='MC10' ;Default is the MC10 model ENDELSE exists=dustem_test_model_exists(use_model) if exists ne 1 then $ message,'Unknown dust model' 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 MC10 ;=== Here we fit the ;=== intensity of the dust-heating radiation field as well as several plug-ins: ;=== (i) a modified black body (initial values correspond to a ; blackbody peaking in the far-infrared) ;=== The free parameters are all lower-bounded at zero. ;=== We invoke the MC10 dust model, and fit the abundances of the PAH components. The BG ;=== components are fixed to ~zero since we want to fit this part of the ;=== observed SED by a MBB. ;use_model='MC10' ; you should specify this when you run the routine! 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_plugin_mbbdy_1', $ ;Amplitude of mBB 'dustem_plugin_mbbdy_2', $ ;T of mBB 'dustem_plugin_mbbdy_3'] ;Emissivity index of mBB iv = [1.8, 7.8e-4, 7.8e-4, 4.e-5, 22., 1.8] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(1.e-15,Npar) fpd = [ $ '(*!dustem_params).grains(2).mdust_o_mh',$ ;amCBEx '(*!dustem_params).grains(3).mdust_o_mh',$ ;amCBEx '(*!dustem_params).grains(4).mdust_o_mh'] ;aSilx fiv = [1.e-12,1.e-12,1.e-12] 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 ;=== READ EXAMPLE SED DATA dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/' file=dir+'example_SED_4.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 (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,m_fit=sed,m_show=sed ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE 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-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='MBB FIT 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,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 ENDIF ;=== MAKE THE FINAL PLOT IF keyword_set(postscript) THEN BEGIN ;dir_ps='./' mydevice=!d.name set_plot,'PS' ;ps_file=dir_ps+postscript ps_file=postscript device,filename=ps_file,/color ENDIF IF !dustem_noobj THEN BEGIN dustemwrap_plot_noobj,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)' ENDIF ELSE BEGIN dustemwrap_plot,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)' ENDELSE IF keyword_set(postscript) THEN BEGIN set_plot,'X' set_plot,mydevice message,'Wrote '+ps_file,/info ENDIF IF keyword_set(wait) THEN BEGIN message,'Made the plot of the final results',/info wait,wait ENDIF IF keyword_set(fits_save) THEN BEGIN message,'Writing out results structure: '+fits_save,/info dustem_write_fits_table,filename=fits_save,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 ;;=== Moved following to the dustem_fitsio_example ;; dustem_read_fits_table,filename=fits_save,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,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)' ;; ENDIF ELSE BEGIN ;; dustemwrap_plot,res,st=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,/info wait,wait ENDIF ENDIF the_end: message,'Finished dustem_fit_intensity_mbb_example',/info END