PRO dustem_fit_sed_ext_pol_example,model=model $ ,sed_file=sed_file $ ,ext_file=ext_file $ ,Nitermax=Nitermax $ ,fits_save_and_restore=fits_save_and_restore $ ,wait=wait $ ,help=help ;+ ; NAME: ; dustem_fit_sed_ext_pol_example ; ; PURPOSE: ; This routine is an example of how to fit an observational emission + extinction SED ; (StokesI,Q,U) 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. ; ; See the DustEMWrap User Guide for more information. ; ; CATEGORY: ; DustEMWrap, Distributed, High-Level, User Example ; ; CALLING SEQUENCE: ; dustem_fit_ext_pol_example[,model=][sed_file=][ext_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_2.xcat' is used. ; ext_file = string naming the path to text file in .xcat format that ; describes the observational extinction SED. If not set, the file ; 'Data/EXAMPLE_OBSDATA/Mathis1990_DISM.xcat' is used. ; 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_sed_ext_pol_example ; dustem_fit_sed_ext_pol_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits' ; dustem_fit_sed_ext_pol_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_sed_ext_pol_example' goto,the_end END IF keyword_set(model) THEN BEGIN use_model=strupcase(model) ENDIF ELSE BEGIN use_model='G17_MODELA' ;Default is last dustem model ENDELSE use_polarization=0. ;default is no polarization in models ;IC: NOT NEEDED SINCE PLOTTING HAPPENS INTERNALLY ;Leaving as is because maybe it's used for something else? use_window=2 ; default graphics window number to use for plotting the results use_verbose=0 if keyword_set(verbose) then use_verbose=1 ;=== 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 CASE use_model OF 'DBP90':BEGIN 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_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0, 4.3e-4, 4.7e-4,6.4e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'DL01':BEGIN 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', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'WD01_RV5p5B':BEGIN ;; ***COMMENT AH*** ;; we need to implement this, or remove message, 'WD01 model not yet implemented in DustEMWrap',/info ;; ***END COMMENT AH*** END 'DL07':BEGIN 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', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001];,10,1.,10.,1.] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'MC10':BEGIN pd = [ $ ;'(*!dustem_params).gas.G0', $ ;G0 ; 'dustem_plugin_continuum_2', $ ;intensity of NIR continuum ; 'dustem_plugin_synchrotron_2',$ '(*!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 ] ;initial parameter values for parameters to be fitted rv = [ $ ; 1.0, $ ; 0.002, $ ;intensity of NIR continuum ; 0.01,$ 7.8e-4, $ ;mass fraction of PAH0 7.8e-4, $ ;mass fraction of PAH1 1.65e-4, $ ;mass fraction of amCBEx 1.45e-3, $ ;mass fraction of amCBEx 7.8e-3 $ ;mass fraction of aSilx ] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) ;=== Fixed parameters ;parameter description of parameters to be set to a non-default value ; fpd=[ $ ; '(*!dustem_params).gas.G0' , $ ;multiplicative factor to total ISRF ; 'dustem_plugin_continuum_2' $ ;intensity of NIR continuum ; ] ; ;initial parameter values for fixed parameters ; fiv=[ $ ; 1. , $ ;multiplicative factor to total ISRF ; 3.e-3$ ;intensity of NIR continuum ; ] ;; ***END COMMENT AH*** END 'J13':BEGIN pd = [ $ '(*!dustem_params).G0', $ ;G0 '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 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', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0, 7.8e-4, 7.8e-4,1.65e-4,1.45e-3,7.8e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) END 'G17_MODELA':BEGIN 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',$ ;aSil 'dustem_plugin_modify_dust_polx_2' ] ;'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [5.4e-4, 5.4e-4,1.8e-4,50] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) use_polarization=0 ;=== Fixed parameters ;parameter description of parameters to be set to a non-default value fpd=[ $ '(*!dustem_params).gas.G0' , $ ;multiplicative factor to total ISRF 'dustem_plugin_continuum_2' $ ;intensity of NIR continuum ] ;initial parameter values for fixed parameters fiv=[ $ 1. , $ ;multiplicative factor to total ISRF 3.e-3$ ;intensity of NIR continuum ] END 'G17_MODELB':BEGIN 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', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) use_polarization=1 END 'G17_MODELC':BEGIN 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', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) use_polarization=1 END 'G17_MODELD':BEGIN 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', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum rv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001] Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) llims=replicate(0.,Npar) use_polarization=1 END 'ELSE':BEGIN message,'model '+model+' unknown',/continue message,'Known models are MC10,DBP90,DL01,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue stop END ENDCASE ;== INITIALISE DUSTEM ;; ***COMMENT AH*** ;; do we need the use_polarization key word activated here? --> ;; len(kwords) must match Npar_dust (i.e. without plugins), otherwise ;; fails in dustem_read_grain -- should set this within the above case statement? ;; IC: yes good idea. Maybe by default set 'kwords' to the original dust keywords and let the user modify kwords(i) where i is the index of the dust species. dustem_init,model=use_model,/pol;,kwords=['logn','plaw','plaw'];,polarization=use_polarization ;; ***END COMMENT AH*** !dustem_nocatch=1 !dustem_verbose=use_verbose ;!dustem_show_plot=1 ;!dustem_noobj=1 !EXCEPT=2 ;so I can locate the plotting error... ;=== READ EXAMPLE DATA: EXTINCTION ;NB:HERE WE ARE NOT FITTING THE DATA. WE'RE JUST USING THE EXISTING FILTERS. dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/' file=dir+'Mathis1990_DISM.xcat' ;We should add the serkowski xcat's content into the MATHIS one right? if keyword_set(ext_file) then file=ext_file ; not really needed for the tests ext=read_xcat(file,/silent) dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/' file=dir+'example_SED_2.xcat' if keyword_set(sed_file) then file=sed_file ; not really needed for the tests sed=read_xcat(file,/silent) ;=== GENERATE EXAMPLE DATA: EXTINCTION ;filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3'] ;Nfilt=n_elements(filters) ;ext=dustem_initialize_ext(Nfilt) ;ext.filter=replicate('SPECTRUM',Nfilt) ;Using these filters to set the wavelengths ;ext.wave=dustem_filter2wav(filters) ;ext.instru=replicate('EXTINCTION', Nfilt);don't really know what to put here ;dustem_filter2instru(filters) ;=== initializing IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro ext.EXT_I=1. for i=4l,n_tags(ext)-1 do begin ext.(i) = ext.EXT_I/1E13 endfor ;=== initializing IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro sed.StokesI=1. for i=4l,n_tags(sed)-1 do begin sed.(i) = sed.StokesI/1E10 endfor ; ind=where(spec.sigextII EQ 0.,count) ; IF count NE 0 THEN spec[ind].sigextII=(0.2*spec(ind).EXT_I)^2 ; ind=where(spec.instru EQ 'FIRAS',count) ; I do not know what this is for (WHY FIRAS)... maybe an old remnant of the older procedures ; IF count NE 0 THEN spec[ind].sigextII=(0.2*spec(ind).EXT_I)^2 ;###Before merging any xcat files we'll just displace the initial parameter array. ;###Creating fake extinction data: (already replaced iv by rv) ;###Creating observational strcture using xcat file (for simplicity) ;###Initializing the structure tags so that dustem_set_data sets said tags in !dustem_data ; for i=4l,n_tags(spec) -1 do begin ; ; ;initial loop ; ;actually only data pertaining to Q and U need to be set. ; ;but dustem_set_data will take care of the tags of the !dustem_data structure pertaining to other polarization quantities ; ;The content of the xcat file shouldn't not/isn't used after this in the main procedure ; if (tag_names(spec))[i] ne 'SIGEXTII' then begin ; spec.(i) = spec.(3)/1E12 ;just to fill the tags... ; ;This gives really high sigma values ; endif ; ; endfor ;###first call to dustem_set_data so that the !dustem_data tags that are used in the 'compute_' procedures are set ;dustem_set_data,m_fit,m_show,ext,ext dustem_set_data,sed,sed,ext,ext iv = rv+[10.00E-4,10.00E-4,10.00E-4,-30] ;###setting the initial parameter vector ;== 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=fpd dustem_init_fixed_params,fpd,fiv ;generation of the fake data (since dustem_set_data takes ONE structure containing all the info) ;!dustem_psi_ext=25 ;using the plugin instead ;stop ext.EXT_I = dustem_compute_ext(rv,st=st) toto = dustem_compute_stokext(rv,st=st,Qext_spec=dustem_qext,Uext_spec=dustem_uext) ext.EXT_Q = toto[0] ext.EXT_U = toto[1] sed.StokesI = dustem_compute_sed(rv,st=st) toto = dustem_compute_stokes(rv,st=st,Q_spec=dustem_qsed,U_spec=dustem_used) sed.StokesQ = toto[0] sed.StokesU = toto[1] ;; ***END COMMENT AH*** ;== SET THE OBSERVATIONAL STRUCTURE dustem_set_data, sed,sed,ext,ext ; dustem_init_fixed_params,fpd,fiv ;== RUN THE FIT tol=1.e-10 use_Nitermax=30 ;maximum number of iterations. IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax xr_x = [0.01,30] yr_x = [5.00E-8,10] xr_m = [1.,5e5] yr_m = [5e-8,1.00e6] ;Because we're in emission+extinction ;two seperate titles have to be provided ;if only one is present, the default titles are kept. tit_m='Spectral Energy Distribution' ; combinations that work: anthing that starts with 'tit' and ends with '_m' or '_emission' tit_x='Dust Optical Depth' ;same for emission but ends with '_x' or '_extinction'. Not case sensitive. ;These two are not needed anymore for any of the plottings.... ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') xtit=textoidl('\lambda (\mum)') ;Set show_plot to 0 to hide plot ;Commented or set to 1 is the same since !dustem_show_plot (existing sysvar) is initialized to 1 in dustem_init ;show_plot = 0 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,tit_x=tit_x $ ,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(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 ;stop 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,tit_m=tit_m+' (From Saved FITS file)',tit_x=tit_x+' (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)',tit_x=tit_x+' (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_sed_ext_pol_example',/info the_end: END