PRO dustem_fit_sed_readme,model=model $
                          ,sed=sed $
                          ,itermax=itermax $
                          ,postscript=postscript $
                          ,png=png $
                          ,save=save $
                          ,restore=restore $
                          ,help=help

;This routine is an example of how to fit observational SEDs 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 data with spectral data, polarisation data and extinction data
;are provided in other _readme routines in the src/idl/ directory.

;;  ***COMMENT AH***
;; Remove or update the next few lines -->>
;; No spectrum data points.
;;SPECTRUM data points can be included and the corresponding filter
;;filed must read SPECTRUM. Note that its is note necessary
;;to use the .xcat file format, and data SED can be provided
;;manually, but the observation structure must have the structure shown below.
;;For the example to work, the DustEM and DustEMWrap packages must have
;;been configured and installed succesfully.
;;See dustem_cvs_readme.txt for install instructions).
;; *** END COMMENT AH***
  
;+
; NAME:
;    dustem_fit_sed_readme  
; PURPOSE:
;    This is an example of how to fit SEDs with DustEMWrap.
;    It is intended as an example to follow when writing your own
;    programs to analyse data with DustEMWrap.
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
; CALLING SEQUENCE:
;    dustem_fit_sed_readme,model=model,sed=sed,save=save,postscript=postscript,png=png,itermax=itermax,help=help
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    Plots, Results save structure
; 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 = 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
;    png = if set, final plot is saved as png  in the
;          current working directory
;    save = if set, SED fits results saved as IDL .sav file in the
;           current working directory
;    itermax = maximum number of fit iterations. Default is 5.
;    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
; PROCEDURES AND SUBROUTINES USED:
;    *** COMMENT AH --> is this really NONE? ****
; EXAMPLES
;    dustem_fit_sed_readme,model='DBP90',sed='../obs/input_sed.xcat',png='example.png',save='SED_fitresults.sav'
; 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_sed_readme'
  goto,the_end
END

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

use_polarization=0.   ;default is no polarization in models
use_window=2       ; default graphics window number to use for plotting the results

;=== Set the (model-dependent) parameters that you want to fit
;=== 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',$     ;PAH1 mass fraction
             '(*!dustem_params).grains(2).mdust_o_mh', $    ;amCBEx
             '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)
     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)
      END
      'WD01_RV5p5B':BEGIN
;;  ***COMMENT AH***
;; we need to implement this, or remove         
         message, 'WD01 model not yet implemented in DustEMWrap',/info
;;  ***END COMMENT AH***
       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,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001];,10,1.,10.,1.]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
      END
      'MC10':BEGIN
        
        pd = [ $
              ;'(*!dustem_params).gas.G0', $      ;G0
             ; 'dustem_plugin_continuum_2', $      ;intensity of NIR continuum
             ; 'dustem_plugin_synchrotron_2',$
              '(*!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 parameter values for parameters to be fitted   
        iv =   [ $  
   ;       1.0, $
  ;       0.002,   $ ;intensity of NIR continuum
  ;        0.01,$    
          7.8e-4,  $ ;mass fraction of  PAH0  
          7.8e-4,  $ ;mass fraction of  PAH1  
          1.65e-4, $ ;mass fraction of  amCBEx
          1.45e-3, $ ;mass fraction of  amCBEx
          7.8e-3   $   ;mass fraction of aSilx
          ]
          
     
        Npar=n_elements(pd)
        
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        
        
        ;=== Fixed parameters
        ;parameter description of parameters to be set to a non-default value
        fpd=[ $
          '(*!dustem_params).gas.G0' ,  $     ;multiplicative factor to total ISRF
          'dustem_plugin_continuum_2' $     ;intensity of NIR continuum
            ]
        ;initial parameter values for fixed parameters
        fiv=[ $
          1. , $        ;multiplicative factor to total ISRF
          3.e-3$      ;intensity of NIR continuum
            ]
;;  ***END COMMENT AH***
      END
     'J13':BEGIN
        pd = [ $
             '(*!dustem_params).G0', $      ;G0
             '(*!dustem_params).grains(1).mdust_o_mh',$     ;PAH1 mass fraction
             '(*!dustem_params).grains(0).mdust_o_mh',$  ;PAH0 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', $      ;aSil
             'dustem_plugin_continuum_2']                ;Intensity of NIR continuum
        iv =   [1.0, 7.8e-4, 7.8e-4,1.65e-4,1.45e-3,7.8e-3,0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
      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',$     ;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)
        use_polarization=0
      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',$     ;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)
        use_polarization=1
     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',$     ;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)
        use_polarization=1
      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',$     ;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)
        use_polarization=1
     END
      'ELSE':BEGIN
         message,'model '+model+' unknown',/continue
         message,'Known models are MC10,DBP90,DL01,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
         stop
      END
ENDCASE



;== INITIALISE DUSTEM
;;  ***COMMENT AH***
;; do we need the use_polarization key word activated here? -->
;; len(kwords) must match Npar_dust (i.e. without plugins), otherwise
;; fails in dustem_read_grain -- should set this within the above case statement?
;; IC: yes good idea. Maybe by default set 'kwords' to the original dust keywords and let the user modify kwords(i) where i is the index of the dust species.
dustem_init,model=use_model;,kwords=['plaw-ed','logn','plaw-ed','plaw-ed','plaw-ed'];,polarization=use_polarization
;;  ***END COMMENT AH***
!dustem_nocatch=1
!dustem_verbose=1
!dustem_show_plot=1


;=== READ EXAMPLE DATA
dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
file=dir+'example_SED_1.xcat'
if keyword_set(sed) then file=sed
spec=read_xcat(file,/silent)

;;  ***COMMENT AH***
;; explain what we are doing here -->
ind=where(spec.sigmaII EQ 0.,count)
IF count NE 0 THEN spec[ind].sigmaII=(0.2*spec(ind).StokesI)^2
ind=where(spec.instru EQ 'FIRAS',count)
IF count NE 0 THEN spec[ind].sigmaII=(0.2*spec(ind).StokesI)^2
;;  ***END COMMENT AH***

;== SET THE OBSERVATIONAL STRUCTURE
st=dustem_set_data(m_fit=spec,m_show=spec);sed=spec)

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

stop
dustem_init_fixed_params,fpd,fiv


stop
;== RUN THE FIT
tol=1.e-16
use_Nitermax=5 ;maximum number of iterations. 
IF keyword_set(itermax) THEN use_Nitermax=itermax

yr=[1.00e-4,1.00E2]
xr=[1.00E0,6.00e4]
tit='Spectral Energy Distribution'
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
xtit=textoidl('\lambda (\mum)')

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)

;=== SAVE FIT RESULTS
;file_out='/tmp/DUSTEM_fit_example.sav'
if keyword_set(save) then begin
   dir_sav='./'
   dustem_save_system_variables,dir_sav+save
   message,'Saved fit results in '+dir_sav+save,/continue
   if not keyword_set(restore) then $
      restore=dir_sav+save
end
message,'The fit executed in '+strtrim(t2-t1,2)+' sec',/info

;======================================
;====You could exit IDL here. The remaining lines of code (essentially
;====plotting the results) would work by returning to this point
;====restoring the IDL output file that was created above.
;======================================

;file='/tmp/DUSTEM_fit_example.sav'
if keyword_set(restore) then begin
   file=restore
   dustem_restore_system_variables,file
end
   
;== PLOT THE FIT RESULTS RESTORED FROM .SAV FILE
; tit='Spectral Energy Distribution (Saved)'
; ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
; xtit=textoidl('\lambda (\mum)')
; errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
; chi2=(*!dustem_fit).chi2
; rchi2=(*!dustem_fit).rchi2
; 
; res=*(*!dustem_fit).current_param_values
; chi2=(*!dustem_fit).chi2
; rchi2=(*!dustem_fit).rchi2
; errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)

;=== PLOT FIT RESULTS AND SAVE TO GRAPHICS FILE IF REQUESTED
; window,use_window
; loadct,13
IF keyword_set(postscript) THEN BEGIN
;  dir_ps=!dustem_dat+'/Figures/'
  dir_ps='./'
;  force_mkdir,dir_ps
  set_plot,'PS'
  ps_file=dir_ps+postscript
  device,filename=ps_file,/color
ENDIF

dustemwrap_plot,res,stp,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit

; dustem_sed_plot,*(*!dustem_fit).current_param_values $
;                 ,ytit=ytit,xtit=xtit,title=tit $
;                 ,yr=yr,xr=xr,/ysty,/xsty $
;                 ,res=res,errors=errors,chi2=chi2,rchi2=rchi2 $
;                 ,/xlog,/ylog,legend_xpos=legend_xpos,legend_ypos=legend_ypos


IF keyword_set(postscript) THEN BEGIN
  device,/close
  set_plot,'X'
  message,'Wrote '+ps_file,/info
  stop
ENDIF
IF keyword_set(png) THEN BEGIN
;  dir_png=!dustem_dat+'/Figures/'
  dir_png='./'
;  force_mkdir,dir_png
  file_png=dir_png+png
  write_png,file_png,tvrd(/true)
  message,'Wrote '+file_png,/info
ENDIF

message,'Finished dustem_fit_sed_readme',/info

the_end:

END