PRO dustem_fit_polarization_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_polarization_example ; ; PURPOSE: ; This is an example of how to fit observational SEDs (with ; measurements in Stokes IQU) using DustEMWrap. ; It is intended as an example to follow when writing your own ; programs to analyse data with DustEMWrap. ; ; CATEGORY: ; DustEMWrap, Distributed, High-Level, User Example ; ; CALLING SEQUENCE: ; dustem_fit_polarization_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 save structure in binary FITS table format ; ; ACCEPTED KEY-WORDS: ; model = specifies the interstellar dust mixture used by DustEM ; 'G17_ModelA' model A from Guillet et al (2018). 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). Default. ; 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_3.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, saves results in a fits file, ; restore the file and plot restored results ; 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_polarization_example ; dustem_fit_polarization_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits' ; dustem_fit_polarization_example,model='DBP90' ; ; MODIFICATION HISTORY: ; Written by JPB 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_polarization_example' goto,the_end END IF keyword_set(model) THEN BEGIN use_model=strupcase(model) ENDIF ELSE BEGIN use_model='G17_MODELD' ;Default is one of the Guillet et al (2017) models since they treat dust polarization properties ENDELSE known_mdls=['MC10','DBP90','DL01','WD01_RV5P5B','DL07','J13','G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD'] pol_mdls=['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 test_model = where(pol_mdls eq use_model,polct) if polct eq 0 then begin message,'The only models with polarization are G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue stop end 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 ;; ;=== AN EXAMPLE FOR G17_MODELA/B/C/D ;; ;=== Here we fit the dust abundances of the model, the dust ;; ;=== polarization and a plug-in: ;; ;=== (i) synchrotron emission (also polarized) ;; ;=== The free parameters in the fit are all lower-bounded at zero. The ;; ;=== intensity of the dust-heating radiation field is fixed to 1.0*G0. pd = [ $ '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;amorphous carbon mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;amorphous silicate mass fraction 'dustem_plugin_modify_dust_pol_2', $ ;Dust polarization angle 'dustem_plugin_synchrotron_2', $ ;Synchrotron amplitude at 10 mm 'dustem_plugin_synchrotron_3', $ ;Synchrotron polarization fraction 'dustem_plugin_synchrotron_4' $ ;Synchrotron polarization angle ] iv=[7.2e-4,7.2e-4,8.5e-4,-10.,0.015,0.25,35.] Npar=n_elements(pd) ulimed=[0,0,0,0,1,1,1] llimed=[1,1,1,0,1,1,1] ulims=[0,0,0,0,0.02,0.5,60.] llims=[0,0,0,0,0.005,0.05,30.] fpd=[ '(*!dustem_params).G0' , $ ; intensity of the dust-heating radiation field 'dustem_plugin_synchrotron_1'] ;Synchrotron CR spectral index fiv=[1.0,3.0] use_polarization=1 ;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 ;== INITIALISE DUSTEM ;dustem_init,mode=use_model,grain_keywords=['logn-chrg-spin','?','plaw-pol'],polarization=1 dustem_init,model=use_model,pol=use_polarization !dustem_verbose=1 !dustem_show_plot=1 !dustem_nocatch=1 ;=== READ EXAMPLE SED DATA dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/' file=dir+'example_SED_3.xcat' IF keyword_set(sed_file) THEN file=sed_file sed=read_xcat(file,/silent) ;== 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;,tied=tied ;== 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 dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims,pol=use_polarization ;=== RUN fit tol=1.e-16 xtol=1.e-16 Nitermax=5; an initialization IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax xrange=[1.,2.e4] yrange=[1e-4,1.e3] loadct,13 t1=systime(0,/sec) title='FIT POLARIZATION EXAMPLE' res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []',title=title) 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 ;=== 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 dustemwrap_plot,*(*!dustem_fit).CURRENT_PARAM_VALUES,dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=title+' (Final fit)' 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 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 dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=title+' (From saved FITS file)' 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_polarization_example',/info stop the_end: END