PRO dustem_fit_spectro_example,model=model $
                               ,sed_file=sed_file $
                               ,Nitermax=Nitermax $
                               ,fits_save_and_restore=fits_save_and_restore $
                               ,ucfactor=ucfactor $
                               ,wait=wait $
                               ,noobj=noobj $
                               ,help=help $
                               ,verbose=verbose
  
;+
; NAME:
;    dustem_fit_spectro_example  
;
; PURPOSE: This routine is an example of how to fit an observational
; spectrometer data in emission (StokesI) with DustEM and
; DustEMWrap.
;  
; See the DustEMWrap User Guide for more information.
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:
;    dustem_fit_spectro_example[,model=][sed_file=][ext_file=][,Nitermax=][,fits_save_and_restore=][,/help,/wait,/verbose,/noobj]
;
; INPUTS:
;    None
;
; OPTIONAL INPUT PARAMETERS:
;    None
;
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUT PARAMETERS:
;    Plots, results structure in binary FITS table format
;
; ACCEPTED KEY-WORDS:
;    model = specifies the interstellar dust mixture used by DustEM
;           'MC10' model from Compiegne et al 2010 
;           '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_file = string naming the path to text file in .xcat format that
;          describes the observational SED. If not set, the file
;          'Data/EXAMPLE_OBSDATA/example_SED_2.xcat' is used.  
;    Nitermax = maximum number of fit iterations. Default is 5.
;    ucfactor = Multiplicative factor for uncertainties on
;               spectrometer data. Default is 1.
;    fits_save_and_restore = if set, save the fit results in a binary
;               FITS file. The code then restore this file and plots
;               the results using the saved results information.
;    help      = if set, print this help
;    wait      = if set, wait this many seconds between each step of
;                the code (for illustration purposes)
;    verbose      = if set, subroutines will run in verbose mode
;    noobj     = if set, runs with no object graphics
;
; 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:
;
; EXAMPLES
;    dustem_fit_spectro_example
;    dustem_fit_spectro_example,Nitermax=10,fits_save_and_restore='/tmp/mysavefile.fits'
;    dustem_fit_spectro_example,model='DBP90'
;
; MODIFICATION HISTORY:
;    Written by AH Oct 2022
;    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_spectro_example'
  goto,the_end
END

IF keyword_set(model) THEN BEGIN
  use_model=strupcase(model)
ENDIF ELSE BEGIN
  use_model='WD01_RV5P5B'    
ENDELSE

known_mdls=['MC10','DBP90','DL01','WD01_RV5P5B','DL07','J13','G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD'] 
test_model = where(known_mdls eq use_model,ct)
if ct eq 0 then begin
   message,'ISM dust model '+use_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

use_polarization=0.   ;default is no polarization in models
use_window=2       ; default graphics window number to use for plotting the results
use_verbose=0
if keyword_set(verbose) then use_verbose=1
use_ucfactor=1.
if keyword_set(ucfactor) then use_ucfactor=ucfactor
use_Nitermax=5        ; maximum number of iterations for the fit
IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax

dustem_define_la_common

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

;=== Set the (model-dependent) parameters that you want to fit (pd),
;=== their initial values (iv)
;=== and whether they are bounded (ulimed,llimed,llims,ulims).
;=== Fixed parameters (fpd) and their values (fiv) are also set here.
;=== Refer to the DustEM and DustEMWrap User guides for an explanation
;=== of the physical meaning of dust model and plug-in parameters, and
;=== how to specify them.

;=== Example is provided for WD01_RV5P5B
;=== If you choose a different dust model, you will need to adjust the
;=== parameter information accordingly.

pd = [ '(*!dustem_params).G0', $ ; multiplicative factor of the radiation field
     '(*!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

iv = [2., 5.4e-4, 5.4e-4, 1.8e-4, 1.65e-3, 7.e-3] ; starting guesses

Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar) 
llims=replicate(1.e-15,Npar)

;=== Fixed parameters
fpd=[]
;=== Initial parameter values for fixed parameters
fiv=[]
        


if keyword_set(wait) then begin
   message,'Finished setting dust model and plug-in parameters: '+use_model,/info
   wait,wait
end

;== INITIALISE DUSTEM
dustem_init,model=use_model,polarization=use_polarization
!dustem_nocatch=1
!dustem_verbose=use_verbose
IF keyword_set(noobj) THEN !dustem_noobj=1
!EXCEPT=2 ; for debugging


;=== READ EXAMPLE EMISSION DATA TO FIT
dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
file=dir+'example_SED_2.xcat'
if keyword_set(sed_file) then file=sed_file 
sed=read_xcat(file,/silent)

;=== TWEAK THE UNCERTAINTIES ON THE SPECTROMETER DATA TO ADJUST THEIR
;=== RELATIVE WEIGHT IN THE FIT
fidx=where(sed.filter eq 'SPECTRUM',fct)
if fct gt 0 then sed[fidx].sigmaii = use_ucfactor*sed[fidx].sigmaii

;== SET THE OBSERVATIONAL STRUCTURE
dustem_set_data, sed,sed

;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT 
;== AND ACTIVATE ANY PLUGINS 
dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims,polarization=use_polarization

if keyword_set(wait) then begin
   message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info
   wait,wait
end


;== INFORMATION TO RUN THE FIT
tol=1.e-10  ;fit tolerance

;=== INFORMATION TO MAKE THE PLOT
;=== _x/_extinction means extinction plots, _m/_emissions means emission plots
xr_m = [1.,5e5]
yr_m = [5e-8,1.00e6]
tit_m='Spectral Energy Distribution' 

;===  RUN THE FIT
t1=systime(0,/sec)
res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
                      ,/xlog,/ylog,xr_m=xr_m,yr_m=yr_m,xr_x=xr_x,yr_x=yr_x,xtit=xtit,ytit=ytit,tit_m=tit_m $
                      ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
                      ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot)
t2=systime(0,/sec)

if keyword_set(wait) then begin
   message,'Finished running DustEMWrap, using Niters: '+strtrim(string(use_Nitermax),2),/info
   message,'Time taken [sec]: '+sigfig(t2-t1,2,/sci),/info
   wait,wait
end



IF keyword_set(fits_save_and_restore) THEN BEGIN
   message,'Writing out results structure: '+fits_save_and_restore,/info
   dustem_write_fits_table,filename=fits_save_and_restore,help=help
   dustem_read_fits_table,filename=fits_save_and_restore,dustem_st=dustem_spectra_st
                                ;==== plot result taken from the saved fits table
   res=*(*!dustem_fit).CURRENT_PARAM_VALUES
   IF !dustem_noobj THEN BEGIN
      dustemwrap_plot_noobj,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)'
   ENDIF ELSE BEGIN
      dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)'
   ENDELSE
   IF keyword_set(wait) THEN BEGIN
      message,'Saved the results as FITS in the file: '+fits_save_and_restore+', and made a plot using the data in this file',/info
      wait,wait
   ENDIF
ENDIF

message,'Finished dustem_fit_spectro_example',/info

the_end:

END