PRO dustem_run_readme,model=model,help=help,plot_it=plot_it,png=png,postcript=postcript

;+
; NAME:
;    dustem_run_readme
; PURPOSE:
;    This is an example of how to run dustem SEDs with the dustem wrapper.
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_run_readme,model=model,help=help,/postcript,/help,/plot_it
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    model = Selects one of the dust mixture used by dustem
;           'COMPIEGNE_ETAL2010' from Compiegne et al 2010 (default)
;           'DBP90' from Desert et al 1990
;           'DL01' from Draine & Li 2001
;           'DL07' from Draine & Li 2007
;    postcript = if set, plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem fortran code must be installed
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    dustem_run_readme,model='COMPIEGNE_ETAL2010'
; MODIFICATION HISTORY:
;    Written by J.P. Bernard April 1st 2011
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

;This Readme shows how to run dustem from the wrapper
;Obviously, the dustem package must have been installed succesfully (see
;dustem_cvs_readme.txt for install instructions).

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_run_readme'
  goto,the_end
ENDIF

loadct,13

IF keyword_set(model) THEN BEGIN
  use_model=strupcase(model)
ENDIF ELSE BEGIN
  use_model='COMPIEGNE_ETAL2010'    ;Default is last dustem model
ENDELSE

;== This just initializes 
dustem_init,model=use_model
!dustem_verbose=1
!dustem_show_plot=1

IF keyword_set(png) THEN BEGIN
  dir_png=!dustem_dat+'/Figures/'
  force_mkdir,dir_png
ENDIF
IF keyword_set(postcript) THEN BEGIN
  dir_ps=!dustem_dat+'/Figures/'
  force_mkdir,dir_ps
ENDIF

IF not keyword_set(dir) THEN BEGIN
;  CASE getenv('DUSTEM_WHICH') OF
  CASE !dustem_which OF
    'DESERT':dir='DESERT_POST_ISO_MORE/'
    'COMPIEGNE': dir='MC_DAT/'
    'VERSTRAETE': dir='d_3.5/'
    'VERSTRAETE': dir='LV_DAT/' 
    ;'WEB3p8': dir_in='dustem3.8_web/' ; this is the original 
    'WEB3p8': dir_in=!dustem_soft_dir ; I changed this 
    ELSE: dir='les_DAT/'
  ENDCASE
ENDIF
;== This is the input directory selected for Dustem
;dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/'+dir
;dir_in=!dustem_soft_dir+'/src/'+dir
;== The following structure contains the inputs to the model
st_model=dustem_read_all(dir_in)
;help,st_model,/str

;=== You can here modify the model inputs
;=== e.g.
;stop
;st_model.G0=2.
;===== This below is to test spinning. 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 writes inputs to the Fortran
dustem_write_all,st_model,!dustem_dat
;== This runs the Fortran. The ouput structure 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 plots the resulting emission SED
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='DUSTEMWrapper 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(postcript) THEN BEGIN
  file_ps=dir_ps+'Last_dustem_run_readme_nuinuem.ps'
  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 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='DUSTEMWrapper extinction '+use_model
dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
IF keyword_set(postcript) THEN BEGIN
  file_ps=dir_ps+'Last_dustem_run_readme_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
  file_png=dir_png+'Last_dustem_run_readme_nuinuem.png'
  xtit=textoidl('\lambda (\mum)')
  ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
  tit='DUSTEMWrapper 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+'Last_dustem_run_readme_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='DUSTEMWrapper 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