PRO dustem_fit_ext_pol_example,model=model $

                               ,ext_file=ext_file $
                               ,Nitermax=Nitermax $

                                 ,postscript=postscript $
                              ,fits_save=fits_save $

                               ,wait=wait $
                               ,noobj=noobj $
                               ,help=help $
                               ,verbose=verbose

;+
; NAME:
;    dustem_fit_ext_pol_example  
;
; PURPOSE:
; This routine is an example of how to fit an observational extinction SED

; (StokesI,Q,U) with DustEM and DustEMWrap.

;  
; See the DustEMWrap User Guide for more information.
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:

;    dustem_fit_ext_pol_example[,model=][ext_file=][,postscript=][,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_init.pro for more details about
;            available models.

;    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 define the observational structure, and a synthetic
;          observation is generated using the dust model.

;    Nitermax = maximum number of fit iterations. Default is 5.

;    postscript = if set, final plot is saved as postscript file

;    fits_save = if set, save the fit results in a binary

;               FITS file. 

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

;    dustem_fit_ext_pol_example,Nitermax=1,fits_save='/tmp/mysavefile.fits'

;    dustem_fit_ext_pol_example,model='G17_MODELD'

;
; MODIFICATION HISTORY:

;    Written by JPB Apr-2021

;    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_ext_pol_example'

  goto,the_end
END

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

  use_model='G17_MODELD'    

ENDELSE

exists=dustem_test_model_exists(use_model,/pol)
if exists ne 1 then $

   message,'Unknown/unpolarised 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=1        ; 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 G17_MODELD
;=== If you choose a different dust model, you will need to adjust the
;=== parameter information accordingly.

pd = [ $
     '(*!dustem_params).grains(0).mdust_o_mh',$         ;PAH0 mass fraction
     '(*!dustem_params).grains(1).mdust_o_mh',$         ;amCBE mass fraction
     '(*!dustem_params).grains(2).mdust_o_mh',$          ;aSil
     'dustem_plugin_modify_dust_polx_2' ]

true_vals = [7.1e-4, 1.3e-3,5.6e-3,35.] ; these are the true values that will be used to generate the data if no ext_file is specified

iv = true_vals+[5.00E-4,-7.00E-4,8.00E-4,5] ; this is the vector of starting guesses for the fit

Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar) & llimed[3]=0

llims=replicate(1.e-15,Npar)

use_polarization=1

;=== 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 OBSERVATIONAL 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+' model',/info
  ;;=== CREATE ENTRIES FOR ALL STOKES (OTHERWISE NEXT CALL TO DUSTEM_SET_DATA WILL COMPLAIN)
  ;;=== WE ASSUME THAT A USER-PROVIDED EXT CURVE WILL SATISFY ALL REQUIREMENTS
   ext.EXT_I=1.
   for i=4l,n_tags(ext)-1 do begin
      ext.(i) = ext.EXT_I/1.E10
   endfor
ENDIF   
 
;=== HERE WE SET THE DATA SO THAT THE COMPUTE_ FUNCTIONS BELOW HAVE
;=== THE NECESSARY FIELDS AND INFORMATION
;== note that m_fit and m_show do not exist -- this is deliberate
;== since we are not treating emission data

dustem_set_data, m_fit=m_fit,m_show=m_show,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 with parameters and plugins',/info
     wait,wait
  ENDIF

  ;stop

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

   toto = dustem_compute_stokext(true_vals,st=st)

   ext.EXT_Q = toto[0]
   ext.EXT_U = toto[1]

   ;stop

   
;== RESET THE OBSERVATIONAL STRUCTURE
;== note that m_fit and m_show do not exist -- this is deliberate
;== since we are not treating emission data

   dustem_set_data, m_fit=m_fit,m_show=m_show,x_fit=ext,x_show=ext

   IF keyword_set(wait) THEN BEGIN

      message,'Finished generating a synthetic EXT from the model '+use_model,/info

     wait,wait
  ENDIF

end

;stop

;== INFORMATION TO RUN THE FIT

tol=1.e-10

;=== INFORMATION TO MAKE THE PLOTS
;=== _x/_extinction means extinction plots, _m/_emissions means emission plots

xr_x=[0.01,30]

yr_fpolext=[1.00E-10,50]

;yr_fpolext=[1e-4,30]

tit='Dust Optical Depth'

;===  RUN THE FIT

;stop

t1=systime(0,/sec)

res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $

                      ,/xlog,/ylog,xr_x=xr_x,yr_fpolext=yr_fpolext,xtit=xtit,ytit=ytit,title=tit $

                      ,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='./'
    mydevice=!d.name
    set_plot,'PS'
;    ps_file=dir_ps+postscript
    ps_file=postscript
    device,filename=ps_file,/color
ENDIF

;stop

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

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

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

;; Moved the following to dustem_fitsio_example.pro
;; dustem_read_fits_table,filename=fits_save,dustem_st=dustem_st
  ;; ;==== plot result taken from the saved fits table
  ;; res=*(*!dustem_fit).CURRENT_PARAM_VALUES
  ;; !dustem_end=-1; subsequent calls to the fit. 0 is intialization an 1 is for the last iteration. It's initialized to 0 after the creation of the fits file.
  ;; IF !dustem_noobj THEN BEGIN
  ;;    dustemwrap_plot_noobj,res,st=dustem_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  ;; ENDIF ELSE BEGIN
  ;;   dustemwrap_plot,res,st=dustem_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  ;; ENDELSE

  IF keyword_set(wait) THEN BEGIN

     message,'Saved the results as FITS in the file: '+fits_save+', and made a plot using the data in this file',/info

     wait,wait
  ENDIF

ENDIF

the_end:

message,'Finished dustem_fit_ext_pol_example',/info

END