PRO dustem_fit_sed_ext_stokesi_example,model=model $
                                       ,sed_file=sed_file $
                                       ,ext_file=ext_file $
                                       ,Nitermax=Nitermax $
                                       ,fits_save=fits_save $
                                       ,wait=wait $
                                       ,noobj=noobj $
                                       ,help=help $
                                       ,verbose=verbose

;+
; NAME:
;    dustem_fit_sed_ext_stokesi_example  
;
; PURPOSE: This routine is an example of how to fit an observational
; data in emission + extinction SED (StokesI only) with DustEM and
; DustEMWrap.
;  
; See the DustEMWrap User Guide for more information.
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:
;    dustem_fit_sed_ext_stokesi_example[,model=][sed_file=][ext_file=][,Nitermax=][,fits_save=][,/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. See userguide or dustem_test_model_exists.pro
;            for more details about available models in current
;            release. Default for this example is J13.
;    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_3.xcat' is used to
;          identify the filters, and the model is used to generate a
;          predicted SED.  
;    ext_file = string naming the path to text file in .xcat format that
;          describes the observational extinction SED. If not set, the file
;          'Data/EXAMPLE_OBSDATA/Mathis90Fitz99_DISM_NH20.xcat' is
;          used to identify the wavelength range, and the dust model
;          is used to generate the corresponding extinction. 
;    Nitermax = maximum number of fit iterations. Default is 5.
;    fits_save = 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_sed_ext_stokesi_example
;    dustem_fit_sed_ext_stokesi_example,Nitermax=1,fits_save='/tmp/mysavefile.fits'
;    dustem_fit_sed_ext_stokesi_example,model='MC10'
;
; 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_ext_stokesi_example'
  goto,the_end
END

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

exists=dustem_test_model_exists(use_model)
if exists ne 1 then $
   message,'Unknown dust model: '+use_model

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_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.

;; ;=== AN EXAMPLE FOR J13
;; ;=== Here we fit the dust abundances of the J13 model, the
;; ;=== intensity of the dust-heating radiation field as well as two plug-ins:
;; ;=== (i) free-free emission
;; ;=== (ii)continuum due to a blackbody 
;; ;=== The temperature of the blackbody is fixed to 1000K
;; ;=== The free parameters in the fit are all lower-bounded at zero.

;=== If you choose a different dust model, you will need to adjust the
;=== parameter information accordingly.

pd = [ $
     '(*!dustem_params).G0' $                                   ;G0
     ,'(*!dustem_params).grains(0).mdust_o_mh'$                  ;CM20 -- power law size distribution
     ,'(*!dustem_params).grains(1).mdust_o_mh'$                 ;CM20 -- logN size distribution
     ,'(*!dustem_params).grains(2).mdust_o_mh' $                ;aPyM5
     ,'(*!dustem_params).grains(3).mdust_o_mh' $                ;aOlM5
     ,'dustem_plugin_continuum_2']                              ;Intensity at peak of the BB continuum
true_vals =   [1.5,1.8e-3, 6.e-4, 3.5e-3, 3.55e-3, 0.003]         
iv =   [1.2,1.7e-3, 6.3e-4, 2.55e-3, 2.55e-3, 0.001]         
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
llims=replicate(1.e-15,Npar)
fpd=[ 'dustem_plugin_continuum_1'] ; Temperature of the BB
fiv=[1000.]                       
    


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 TEMPLATE EXTINCTION DATA 
;=== IF THE USER DOES NOT PROVIDE AN EXTINCTION CURVE FILE , WE READ MATHIS TO INITIALISE THE OBSERVATIONAL STRUCTURE.
dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
use_ext_file=dir+'Mathis90Fitz99_DISM_NH20.xcat' 
IF keyword_set(ext_file) THEN use_ext_file=ext_file
ext=read_xcat(use_ext_file,/silent)
IF not keyword_set(ext_file) THEN BEGIN
   message, 'You have not provided an extinction curve. I will generate one from the '+use_model+' dust model',/info
   ext.EXT_I=1.
ENDIF   



;=== READ TEMPLATE EMISSION DATA 
;=== IF THE USER DOES NOT PROVIDE A SED FILE , WE READ AN EXAMPLE_SED FILE TO INITIALISE THE OBSERVATIONAL STRUCTURE.
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)
IF not keyword_set(sed_file) THEN BEGIN
   message, 'You have not provided an SED. I will generate one from the '+use_model+' dust model',/info
   sed.StokesI=1.
ENDIF   


;=== HERE WE SET THE DATA SO THAT THE COMPUTE_ FUNCTIONS BELOW HAVE
;=== THE NECESSARY FIELDS AND INFORMATION
dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext 

;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT 
;== AND ACTIVATE ANY PLUGINS (HENCE WE NEED TO DO THIS BEFORE THE COMPUTE_ CALLS)
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

;== GENERATE THE EXTINCTION DATA USING THE MODEL AND TRUE VALUES
if not keyword_set(ext_file) then $
   ext.EXT_I = dustem_compute_ext(true_vals,st=st)

;== GENERATE THE EMISSION DATA USING THE MODEL AND TRUE VALUES
IF not keyword_set(sed_file) THEN $
   sed.StokesI = dustem_compute_sed(true_vals,st=st)

;=== RESET THE OBSERVATIONAL STRUCTURE IF WE HAVE GENERATED DATA USING
;=== A MODEL
if not keyword_set(ext_file) or  not keyword_set(sed_file) THEN BEGIN
   dustem_set_data, m_fit=sed,m_show=sed,x_fit=ext,x_show=ext
   if keyword_set(wait) then begin
      message,'Finished generating data from the model '+use_model,/info
      wait,wait
   end
end


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

;=== INFORMATION TO MAKE THE PLOT
;=== _x/_extinction means extinction plots, _m/_emissions means emission plots
xr_x = [0.01,30]
yr_x = [5.00E-8,10]
xr_m = [1.,5e5]
yr_m = [5e-8,1.00e6]
tit_m='Spectral Energy Distribution' 
tit_x='Dust Optical Depth' 

;===  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,tit_x=tit_x $
                      ,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


;=== MAKE THE FINAL PLOT
IF keyword_set(postscript) THEN BEGIN
    dir_ps='./'
    set_plot,'PS'
    ps_file=dir_ps+postscript
    device,filename=ps_file,/color
ENDIF

IF !dustem_noobj THEN BEGIN 
  dustemwrap_plot_noobj,(*(*!dustem_fit).CURRENT_PARAM_VALUES),st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (Final fit)',tit_x=tit_x+' (Final fit)' 
ENDIF ELSE BEGIN
  dustemwrap_plot,(*(*!dustem_fit).CURRENT_PARAM_VALUES),st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (Final fit)',tit_x=tit_x+' (Final fit)' 
ENDELSE

IF keyword_set(postscript) THEN BEGIN
  set_plot,'X'
  device,/close
  message,'Wrote '+ps_file,/info
ENDIF

if keyword_set(wait) then begin
   message,'Made the plot of the final results',/info
   wait,wait
end


IF keyword_set(fits_save) THEN BEGIN
   message,'Writing out results structure: '+fits_save,/info
   dustem_write_fits_table,filename=fits_save,help=help
;=== At this point, you could erase all dustem system variables, or exit idl... all the
;=== information needed to recover the results and remake the plots has been saved in the FITS table

;;    dustem_read_fits_table,filename=fits_save,dustem_st=dustem_spectra_st
;; ;==== plot result taken from the saved fits table
;;    res=*(*!dustem_fit).CURRENT_PARAM_VALUES
;;   ;If the user prefers non-OOP plotting then reset !dustem_noobj=1 here since dustem_init is called in reads_fits_table.
;;    IF !dustem_noobj THEN BEGIN
;;       dustemwrap_plot_noobj,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)',tit_x=tit_x+' (From Saved FITS file)'
;;    ENDIF ELSE BEGIN
;;       dustemwrap_plot,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,tit_m=tit_m+' (From Saved FITS file)',tit_x=tit_x+' (From Saved FITS file)'
;;    ENDELSE

   IF keyword_set(wait) THEN BEGIN
      message,'Saved the results as FITS in the file: '+fits_save,/info
      wait,wait
   ENDIF
ENDIF

the_end:
message,'Finished dustem_fit_sed_ext_stokesi_example',/info

END