PRO dustem_run_example,model=model $
                      ,png=png $
                      ,postscript=postscript $
                      ,help=help

;+
; NAME:
;    dustem_run_example
; PURPOSE:
;    This is an example of how to run the DustEM fortran code using DustEMWrap
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
; CALLING SEQUENCE:
;    dustem_run_example,model=model,postscript=postscript,png=png,help=help
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; 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)
;    postscript = if set, final plot is saved as postscript in the
;                 current working directory
;    png = if set, final plot is saved as png  in the
;          current working directory
;    help      = if set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    dustem_run_example,model='DBP90',png='test_dbp90.png'
; MODIFICATION HISTORY:
;    Written 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_run_example'
  goto,the_end
ENDIF

loadct,13

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

;== INITIALISE DUSTEM
dustem_init,model=use_model
!dustem_verbose=1
!dustem_show_plot=1

CASE !dustem_which OF
   'DESERT': begin
      dir='DESERT_POST_ISO_MORE/'
      message,'DESERT version disabled for this routine',/info
   end
   'COMPIEGNE': begin
      dir='MC_DAT/'
      message,'COMPIEGNE version disabled for this routine',/info
   end
   'VERSTRAETE': begin
      dir='LV_DAT/'
      message,'VERSTRAETE version disabled for this routine',/info
   end
   'RELEASE': dir_in=!dustem_soft_dir  
   ELSE: dir_in=!dustem_soft_dir  
   ;dir='les_DAT/'
ENDCASE

;== Fill the following structure contains with default inputs to the model
st_model=dustem_read_all(dir_in)
;help,st_model,/str

;=== You can modify the model inputs by changing the values in the
;=== st_model structure

;=== e.g. modify the default value of the ISRF intensity used by DustEM
;st_model.G0=2.

;===== e.g. modify keywords related to spinning dust. Does not work with DBP90 !
;st_model.grains[0].type_keywords='logn-chrg-spin'
;st_model.gas.NH=1.e2
;st_model.gas.Tgas=1.e5

;== The following line saves the modified inputs for use by the Fortran
dustem_write_all,st_model,!dustem_dat

;== The following line runs the Fortran. The ouput structure st contains the results
st=dustem_run()

;== Optionaly, this creates an SED to be overploted
;stop
;filters=['IRAS1','IRAS4']
;Nfilt=n_elements(filters)
;sed=dustem_initialize_sed(Nfilt)
;sed.filter=filters
;sed.wave=dustem_filter2wav(filters)
;sed.instru=dustem_filter2instru(filters)
;stop
;stt=dustem_set_data(sed=sed)

;== The following lines plot the resulting SED (emission)
win=0L
window,win & win=win+2
xtit=textoidl('\lambda (\mum)')
;=== same plot as in Compiegne et al. 2010
ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
tit='DustEMWrap Emitted Intensity '+use_model
yr=[1e-11,6.e-7]
xr=[1,5e3]
dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit
IF keyword_set(postscript) THEN BEGIN
   dir_ps='./'
   file_ps=dir_ps+postscript
   previous_device=!d.name
   set_plot,'PS'
   device,filename=file_ps
   dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit
   device,/close
   set_plot,previous_device
   message,'Wrote '+file_ps,/continue
ENDIF

;== The following plots the resulting SED (extinction)
window,win & win=win+1
yr=[0,2.5]  ;range of 1/lambda
xr=[1,10]   ;range of sigma
xtit=textoidl('1/\lambda (\mum^{-1})')
ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)')
tit='DustEMWrap extinction '+use_model
dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
IF keyword_set(postscript) THEN BEGIN
  file_ps=dir_ps+file_basename(file_ps,'.ps')+'_extinction.ps'
  previous_device=!d.name
  set_plot,'PS'
  device,filename=file_ps
  dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
  device,/close
  set_plot,previous_device
  message,'Wrote '+file_ps,/continue
ENDIF

IF keyword_set(png) THEN BEGIN
  win=10L
  window,win ;& win=win+2
  dir_png='./'
  file_png=dir_png+png
  xtit=textoidl('\lambda (\mum)')
  ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
  tit='DustEMWrap Emitted Intensity '+use_model
  yr=[1e-11,6.e-7]
  xr=[1,5e3]
  dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit
  write_png,file_png,tvrd(/true)
  message,'Wrote '+file_png,/info
  ;=====
  file_png=dir_png+file_basename(file_png,'.png')+'_extinction.png'

  yr=[0,2.5]  ;range of 1/lambda
  xr=[1,10]   ;range of sigma
  xtit=textoidl('1/\lambda (\mum^{-1})')
  ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)')
  tit='DustEMWrap extinction '+use_model
  dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
  write_png,file_png ,tvrd(/true),/verbose
  message,'Wrote '+file_png,/info
ENDIF

the_end:

END