PRO dustem_variational_study_readme,model=model,help=help,png=png,postscript=postscript

;+
; NAME:
;    dustem_variational_study_readme
; PURPOSE:
;    This is an example of how to run a variational study on the dust model parameters with the dustem wrapper.
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_variational_study_readme[,model=][,/help][,/png][,/postscript]
; 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
;    postscript = 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_variational_study_readme,model='DBP90',/postscript
;    dustem_variational_study_readme,model='COMPIEGNE_ETAL2010',/postscript
; MODIFICATION HISTORY:
;    Written by J.P. Bernard July 15 2021
;    see evolution details on the dustem gitlab
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@irap.omp.eu) in case of problems.
;-

;This Readme shows how to run dustem from the wrapper

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_variational_study_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 wrapper
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(postscript) 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=!dustem_soft_dir
    ELSE: dir='les_DAT/'
  ENDCASE
ENDIF
;== The following structure contains the inputs to the model
st_model=dustem_read_all(dir_in)

;These are the values of G0 to be explored
g0s=[1.,2.,5.,10.,100.,1000.]
n_g0s=n_elements(g0s)

win=0L
window,win,title='DUSTEMWrapper Variational study'

CASE use_model OF
  'DBP90': pd = ['(*!dustem_params).G0'] ;This is the parameter description of G0 for models without gas parameters
  'COMPIEGNE_ETAL2010': pd = ['(*!dustem_params).gas.G0'] ;This is the parameter description of G0 for models with gas parameters
  'AJ13' : pd = ['(*!dustem_params).G0'] ;This is the parameter description of G0 for models without gas parameters
ENDCASE

FOR i=0L,n_g0s-1 DO BEGIN
  iv = [g0s[i]]
  dustem_init_parinfo,pd,iv
  params=[g0s[i]]
  st=dustem_run(params)
  ;do the plot
  cleanplot
  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+' G0='+strtrim(g0s[i],2)
  yr=[1e-11,1.e-4]
  xr=[1,5e3]
  ;stop
  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
    file_ps=dir_ps+use_model+'_variational_study_nuinuem_G0='+strtrim(g0s[i],2)+'.ps'
    spawn,'rm '+file_ps
    previous_device=!d.name
    set_plot,'PS'
    device,filename=file_ps
    !p.charsize=0.8
    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

ENDFOR

;stop

the_end:

END