PRO dustem_run_example, model $
                        ,show=show $
                        ,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,show=show,postscript=postscript,help=help
;
; INPUTS:
;        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 D from Guillet et al (2018)
;
; OPTIONAL INPUT PARAMETERS:
;       show = vector specifying the type of plot(s) to show. Possible
;       values (for all dust models) are
;       ["emis", "extuv", "extir", "alb", "sdist","rfield"]
;       Additionally, for the G17_MODEL?s with polarisation, you can specify
;        [ "polext", "polsed", "align"]
;       "emis" -- Stokes I SED in emission predicted by the model
;       "extuv" -- Extinction in the UV predicted by the model
;       "extir" -- Extinction in the IR predicted by the model
;       "alb" -- Albedo of the grain types
;       "sdist" -- Size distribution of the grain types
;       "rfield" -- Dust-heating radiation field
;       "polsed" -- Polarised intensity SED predicted by the model
;       "polext" -- Polarised extinction predicted by the model
;       "align" -- Grain alignment fraction predicted by the model
;       "all" -- All plots appropriate to the model  
;
; OUTPUTS:
;    Plots
;
; OPTIONAL OUTPUT PARAMETERS:
;    None
;
; ACCEPTED KEY-WORDS:
;    postscript = on/off if set, plots are saved as postscript in the
;                 current working directory (and not shown on screen)
;    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,'DBP90',show="all",/post
;
; 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

IF keyword_set(show) THEN BEGIN
  use_show=show
ENDIF ELSE BEGIN
  use_show=['all']    ;Default is to show everything
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 strupcase(model),ct)
if ct eq 0 then begin
   message,'ISM dust 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
test_pol_model = where(pol_mdls eq model,polct)

;== INITIALISE DUSTEM
dustem_init,model=model
!dustem_verbose=1
!dustem_show_plot=1
!dustem_which='RELEASE'
dir_in=!dustem_soft_dir  

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


;=== You could modify the model inputs here by directly editing the values in the
;=== st_model structure. To inspect contents, type IDL> help,st_model,/str
;=== e.g.
;=== st_model.G0=2.0 ; change the default value of the ISRF intensity used by DustEM
;=== st_model.isrf.isrf=sqrt(st_model.isrf.isrf) ; change the shape of ISRF
;=== st_model.grains[0].mdust_O_mh=5.e-4 ; change the abundance of grain[0]
;=== st_model.grains[1].type_keywords='plaw' ; change the size distribution description of grain[1]

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

known_plots=["emis", "extuv", "extir", "alb", "sdist"]
if polct gt 0 then known_plots=["emis", "extuv", "extir", "alb", "sdist","polext", "polsed", "align"]
if n_elements(use_show) eq 1 then $
   if strupcase(use_show) eq "ALL" then use_show=[known_plots,"rfield"]

;== The following lines generate plots that show the outputs of the fortran.
;=== For this release (V2.0), we use the old method (dustem_show_fortran) to
;=== access model quantities.
;=== This method will be deprecated in a future release and plots will be
;=== constructed directly from the st and st_model structures.
match,known_plots,use_show,a,b
if total(a,/nan) ne -1 and keyword_set(postscript) then $
   dustem_show_fortran,model=model,st=st,show=known_plots[a],hard=1,tit="DustEM_run_example "+model+" :"
if total(a,/nan) ne -1 and not keyword_set(postscript) then $
   dustem_show_fortran,model=model,st=st,show=known_plots[a],tit="DustEM_run_example "+model+" :"

;== The following lines generate a plot of the ISRF using the
;== st_model structure. For illustration purposes.
use_win=12
rfield_plot=["rfield"]
match,rfield_plot,use_show,a,b
if total(a,/nan) ne -1 then begin
      xtit=textoidl('\lambda (\mum)')
      ytit=textoidl('log ISRF [4 \pi I_\nu (erg/cm^2/s/Hz)]')
      tit='DustEM_run_example ISRF'
      yr=[min(st_model.isrf.isrf),5.*max(st_model.isrf.isrf)]
      xr=[min(st_model.isrf.lambisrf),max(st_model.isrf.lambisrf)]
   IF not keyword_set(postscript) THEN BEGIN
      window,use_win,xs=600,ys=400,tit='DUSTEM_SHOW_FORTRAN: ISRF '+strtrim(use_win,2) & use_win=use_win+2
      cgplot,st_model.isrf.lambisrf,st_model.isrf.isrf $
             ,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog $
             ,title=tit,xtit=xtit,ytit=ytit,color=cgcolor('black'),/nodata
      oplot,st_model.isrf.lambisrf,st_model.isrf.isrf,col=cgcolor('black')
   END ELSE BEGIN
      file_ps="isrf.ps"
      previous_device=!d.name
       !y.thick = 5
       !x.thick = 5
       !p.thick = 5
       !p.charsize = 1.3
       !p.charthick = 5
       !x.ticklen = 0.04
      set_plot,'PS'
      device,filename=file_ps,/portrait,/color
      cgplot,st_model.isrf.lambisrf,st_model.isrf.isrf $
             ,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog $
             ,title=tit,xtit=xtit,ytit=ytit,color=cgcolor('black'),/nodata
      oplot,st_model.isrf.lambisrf,st_model.isrf.isrf ,col=cgcolor('black')
      device,/close
      !y.thick = 0
      !x.thick = 0
      !p.thick = 0
      !p.charsize = 1
      !p.charthick = 0
      !x.ticklen = 0.02
      set_plot,previous_device
   END
end

goto,the_end     

;;== The following lines plot the resulting SED (emission)
;;== using st/st_model structure and a dedicated plotting function 
;;== Prototype for future releases.
;; win=0L
;; window,win & win=win+2
;; xtit=textoidl('\lambda (\mum)')
;; ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
;; tit='DustEMWrap Emission: '+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,post=postscript

;;== The following plots the resulting SED (extinction)
;;== using st/st_model structure and a dedicated plotting function 
;;== Prototype for future releases.
;; 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: '+model
;; dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit,post=postscript

the_end:

END