PRO dustem_fit_sed_ext_pol_example,model=model $ ,sed_file=sed_file $ ,ext_file=ext_file $ ,Nitermax=Nitermax $ ,fits_save=fits_save $ ,postscript=postscript $ ,wait=wait $ ,noobj=noobj $ ,help=help $ ,verbose=verbose ;+ ; NAME: ; dustem_fit_sed_ext_pol_example ; ; PURPOSE: This routine is an example of how to fit an observational ; data in emission + extinction SED (StokesI,Q,U) with DustEM and ; DustEMWrap. ; ; See the DustEMWrap User Guide for more information. ; ; CATEGORY: ; DustEMWrap, Distributed, High-Level, User Example ; ; CALLING SEQUENCE: ; dustem_fit_sed_ext_pol_example[,model=][sed_file=][ext_file=][postscript=][,Nitermax=][,fits_save=][,/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. 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_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/Mathis90Fitz99_DISM_NH20.xcat' is used. ; postscript = if set, final plot is saved as postscript file ; Nitermax = maximum number of fit iterations. Default is 5. ; fits_save = 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_sed_ext_pol_example ; dustem_fit_sed_ext_pol_example,Nitermax=1,fits_save='/tmp/mysavefile.fits' ; dustem_fit_sed_ext_pol_example,model='G17_MODELC' ; ; 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_MODELD' ENDELSE exists=dustem_test_model_exists(use_model,/pol) if exists ne 1 then $ message,'Unknown/unpolarized dust model: '+use_model 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 ;=== 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 G17_MODELD ;=== If you choose a different dust model, you will need to adjust the ;=== parameter information accordingly. pd = [ $ '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;amCBE mass fraction '(*!dustem_params).grains(2).mdust_o_mh',$ ;aSil 'dustem_plugin_modify_dust_pol_2', $ 'dustem_plugin_modify_dust_polx_2' ] true_vals = [5.4e-4, 6.4e-4,4.8e-4, 48., 52] ; these are the true values that will be used to generate the data iv = true_vals+[2.00E-4,-2.00E-4,2.00E-4,10, -10] ; this is the vector of starting guesses for the fit Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) & llimed[3]=0 & llimed[4]=0 llims=replicate(1.e-15,Npar) use_polarization=1 ;=== Fixed parameters fpd=[ $ '(*!dustem_params).G0' , $ ;multiplicative factor to total ISRF 'dustem_plugin_continuum_2' $ ;intensity of a BB continuum ] ;=== Initial parameter values for fixed parameters fiv=[ $ 1. , $ ;multiplicative factor to total ISRF 3.e-3$ ;intensity of BB continuum ] 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 OBSERVATIONAL EXTINCTION DATA ;=== IF THE USER DOES NOT PROVIDE AN EXTINCTION CURVE FILE , WE READ MATHIS TO INITIALISE THE OBSERVATIONAL STRUCTURE. dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/' use_ext_file=dir+'Mathis90Fitz99_DISM_NH20.xcat' IF keyword_set(ext_file) THEN use_ext_file=ext_file ext=read_xcat(use_ext_file,/silent) IF not keyword_set(ext_file) THEN BEGIN message, 'You have not provided an extinction curve. I will generate one from the '+use_model+' model',/info ;;=== CREATE ENTRIES FOR ALL STOKES (OTHERWISE NEXT CALL TO ;;=== DUSTEM_SET_DATA WILL COMPLAIN SINCE WE HAVE INITIALISED ;;=== DUSTEM IN POLARIZATION MODE ;;=== WE ASSUME THAT A USER-PROVIDED EXT CURVE WILL SATISFY ALL REQUIREMENTS ext.EXT_I=1. for i=4l,n_tags(ext)-1 do begin ext.(i) = 1.E-15 endfor ENDIF ;=== READ OBSERVATIONAL EMISSION DATA ;=== IF THE USER DOES NOT PROVIDE A SED FILE , WE READ AN EXAMPLE_SED FILE TO INITIALISE THE OBSERVATIONAL STRUCTURE. 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) IF not keyword_set(sed_file) THEN BEGIN message, 'You have not provided an SED. I will generate one from the '+use_model+' model',/info ;;=== CREATE ENTRIES FOR ALL STOKES (OTHERWISE NEXT CALL TO ;;=== DUSTEM_SET_DATA WILL COMPLAIN SINCE WE HAVE INITIALISED ;;=== DUSTEM IN POLARIZATION MODE ;;=== WE ASSUME THAT A USER-PROVIDED SED CURVE WILL SATISFY ALL REQUIREMENTS sed.StokesI=1. for i=4l,n_tags(sed)-1 do begin sed.(i) = 1.E-15 endfor ENDIF ;=== HERE WE SET THE DATA SO THAT THE COMPUTE_ FUNCTIONS BELOW HAVE ;=== THE NECESSARY FIELDS AND INFORMATION dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT ;== AND ACTIVATE ANY PLUGINS (HENCE WE NEED TO DO THIS BEFORE THE COMPUTE_ CALLS) 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 if not keyword_set(ext_file) then begin ;== GENERATE THE EXTINCTION DATA USING THE MODEL AND TRUE VALUES ext.EXT_I = dustem_compute_ext(true_vals,st=st) toto = dustem_compute_stokext(true_vals,st=st,Qext_spec=dustem_qext,Uext_spec=dustem_uext) ext.EXT_Q = toto[0] ext.EXT_U = toto[1] end IF not keyword_set(sed_file) THEN BEGIN ;== GENERATE THE EMISSION DATA USING THE MODEL AND TRUE VALUES sed.StokesI = dustem_compute_sed(true_vals,st=st) toto = dustem_compute_stokes(true_vals,st=st,Q_spec=dustem_qsed,U_spec=dustem_used) sed.StokesQ = toto[0] sed.StokesU = toto[1] end ;=== RESET THE OBSERVATIONAL STRUCTURE IF WE HAVE GENERATED DATA USING ;=== A MODEL if not keyword_set(ext_file) or not keyword_set(sed_file) THEN BEGIN dustem_set_data, m_fit=sed,m_show=sed,x_fit=ext,x_show=ext if keyword_set(wait) then begin message,'Finished generating data from the model '+use_model,/info wait,wait end end ;== INFORMATION TO RUN THE FIT tol=1.e-10 ;fit tolerence ;=== INFORMATION TO MAKE THE PLOT ;=== _x/_extinction means extinction plots, _m/_emissions means emission plots xr_x = [0.01,30] yr_x = [5.00E-8,10] xr_m = [1.,5e5] yr_m = [5e-8,1.00e6] tit_m='Spectral Energy Distribution' tit_x='Dust Optical Depth' ;=== 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,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(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='./' mydevice=!d.name set_plot,'PS' 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,tit_m=tit_m+' (Final fit)',tit_x=tit_x+' (Final fit)' ENDIF ELSE BEGIN dustemwrap_plot,(*(*!dustem_fit).CURRENT_PARAM_VALUES),st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (Final fit)',tit_x=tit_x+' (Final fit)' ENDELSE IF keyword_set(postscript) THEN BEGIN device,/close 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 end 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 ;; 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 the user prefers non-OOP plotting then reset !dustem_noobj=1 here since dustem_init is called in reads_fits_table. ;; IF !dustem_noobj THEN BEGIN ;; dustemwrap_plot_noobj,res,st=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,st=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,/info wait,wait ENDIF ENDIF the_end: message,'Finished dustem_fit_sed_ext_pol_example',/info END