dustem_fit_intensity_example.pro 16.4 KB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409
PRO dustem_fit_intensity_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

;This routine is an example of how to fit an observational SED (Stokes
;I only) 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!)
;  
;For this example, the code uses the SED in the file example_SED_1.xcat,
;which is distributed in the Data/EXAMPLE_OBSDATA/ directory 
;
;The example SED has Stokes I photometric data points from SPITZER
;IRAC and 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.
  
;+
; NAME:
;    dustem_fit_intensity_example  
; PURPOSE:
;
;  This routine is an example of how to fit an observational SED
;  (Stokes I only) using 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!)
;  
;  For this example, the code uses the SED in the file example_SED_1.xcat,
;  which is distributed in the Data/EXAMPLE_OBSDATA/ directory 
;
;  The example SED has Stokes I photometric data points from SPITZER
;  IRAC and 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.
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
; CALLING SEQUENCE:
;    dustem_fit_intensity_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save_and_restore=][,/help]
; 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 (default)
;           '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/SEDs/sample_SED.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, 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_intensity_example
;    dustem_fit_intensity_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits'
;    dustem_fit_intensity_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_intensity_example'
  goto,the_end
END

IF keyword_set(model) THEN BEGIN
  use_model=strupcase(model)
ENDIF ELSE BEGIN
  use_model='MC10'    ;Default is the Compeigne et al (2011) model
ENDELSE

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

;=== 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 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',$  ;VSG mass fraction
             '(*!dustem_params).grains(2).mdust_o_mh', $ ;BG mass fraction
             'dustem_plugin_continuum_2']                ;Intensity of NIR continuum
        iv =   [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)
        fpd=[] & fiv=[]
     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
        iv =   [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)
        fpd=[] & fiv=[]
      END
      'WD01_RV5p5B':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
        iv =   [1.0,5.4e-4, 5.4e-4,1.8e-4,1.64e-3,6.99e-3,0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        fpd=[] & fiv=[]
       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
        iv =   [1.0,4.97e-4, 4.97e-4,1.66e-4,2.21e-3,7.64e-3,0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        fpd=[] & fiv=[]
      END
      'MC10':BEGIN
; parameters of the model to fit
         pd = [ $
              '(*!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 values of the model parameters that are to be varied during the fit
         iv =   [ 7.8e-4, 7.8e-4, 1.65e-4, 1.45e-3, 6.7e-3]
        Npar=n_elements(pd)
; flags for whether the parameters are upper/lower bounded
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
; lower bound value
        llims=replicate(0.,Npar)
; parameters of the model that are fixed
        fpd=[ $
          '(*!dustem_params).gas.G0' ,  $   ; ISRF intensity (multiplicative factor x G0)
          'dustem_plugin_continuum_2' $     ;intensity of NIR continuum
            ]
; values of the fixed model parameters
        fiv=[2.4, 3.e-3]
      END
     'J13':BEGIN
        pd = [ $
             '(*!dustem_params).G0', $      ;G0
             '(*!dustem_params).grains(0).mdust_o_mh',$     ;CM20 -- power law size distribution
             '(*!dustem_params).grains(1).mdust_o_mh',$     ;CM20 -- logN size distribution
             '(*!dustem_params).grains(2).mdust_o_mh', $    ;aPyM5
             '(*!dustem_params).grains(3).mdust_o_mh', $    ;aOlM5
             'dustem_plugin_continuum_2']                ;Intensity of NIR continuum
        iv =   [1.0, 1.7e-3, 6.3e-4,2.55e-3,2.55e-3,0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        fpd=[] & fiv=[]
     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',$     ;amCB_0.3333x
             '(*!dustem_params).grains(2).mdust_o_mh', $    ;aSil2001_0.3333_p20B
             'dustem_plugin_continuum_2']                ;Intensity of NIR continuum
        iv =   [1.0, 7.094483e-4, 1.319896e-3, 5.588265e-3, 0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        use_polarization=1
        fpd=[] & fiv=[]
      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',$     ;amCB_0.3333x
             '(*!dustem_params).grains(2).mdust_o_mh', $    ;aSil_0.3333_p20B
             'dustem_plugin_continuum_2']                ;Intensity of NIR continuum
        iv =   [1.0, 6.480655e-4, 1.327672e-3, 5.267745e-3, 0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        use_polarization=1
        fpd=[] & fiv=[]
     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',$     ;amCBE_0.3333x
             '(*!dustem_params).grains(2).mdust_o_mh', $    ;aSil2001_0.3333_p20B
             'dustem_plugin_continuum_2']                ;Intensity of NIR continuum
        iv =   [1.0, 6.979060e-4, 1.555919e-3, 4.366849e-3, 0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        use_polarization=1
        fpd=[] & fiv=[]
      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',$     ;amCBE_0.3333x
             '(*!dustem_params).grains(2).mdust_o_mh', $    ;aSil2001BE6pctG_0.4x
             'dustem_plugin_continuum_2']                ;Intensity of NIR continuum
        iv =   [1.0, 7.235580e-4, 1.222725e-3, 5.324418e-3, 0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        use_polarization=1
        fpd=[] & fiv=[]
     END
      'ELSE':BEGIN
         message,'model '+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
ENDCASE
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

if keyword_set(wait) then begin
   message,'Finished setting dust model: '+model,/info
   wait,wait
end


;== INITIALISE DUSTEM
dustem_init,model=use_model,polarization=use_polarization
!dustem_nocatch=1
!dustem_verbose=use_verbose
!dustem_show_plot=1


;=== READ EXAMPLE SED DATA
dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
file=dir+'example_SED_1.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 / plot. These are often, but not
;== always!, the same SED.
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

;== 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

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
use_Nitermax=5 ;maximum number of iterations. 
IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax

;=== INFORMATION TO MAKE THE PLOT
yr=[1.00e-4,1.00E2] ; y-axis limits
xr=[1.00E0,6.00e4] ; x-axis limits
tit='SED -- Example 1' ; 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)
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=tit+' (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 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=tit+' (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_fit_intensity_example',/info

the_end:

END