PRO dustem_fit_polarisation_example,model=model $

                                    ,sed_file=sed_file $
                                    ,Nitermax=Nitermax $
                                    ,postscript=postscript $
                                    ,fits_save_and_restore=fits_save_and_restore $
                                    ,help=help $
                                    ,wait=wait $
                                    ,verbose=verbose

  ;+

; NAME:
;    dustem_fit_polarisation_example  

;

; PURPOSE:

;    This is an example of how to fit observational SEDs (with
;    measurements in Stokes IQU) using 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_polarisation_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save_and_restore=][,/help,/wait,/verbose]
;

; INPUTS:
;    None

;

; OPTIONAL INPUT PARAMETERS:
;    None

;

; OUTPUTS:
;    None

;

; OPTIONAL OUTPUT PARAMETERS:

;    Plots, Results save structure in binary FITS table format
;

; ACCEPTED KEY-WORDS:
;    model = specifies the interstellar dust mixture used by DustEM

;           'G17_ModelA' model A from Guillet et al (2018). 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). Default.

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

;    postscript = if set, final plot is saved as postscript in the
;                current working directory
;    Nitermax = maximum number of fit iterations. Default is 5.
;    fits_save_and_restore = if set, saves results in a fits file, restore the file and plot restored results
;    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

; 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_polarisation_example
;    dustem_fit_polarisation_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits'
;    dustem_fit_polarisation_example,model='DBP90'

;

; MODIFICATION HISTORY:
;    Written by JPB Apr-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_polarisation_example'
  goto,the_end
END

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

  use_model='G17_MODELD'    ;Default is one of the Guillet et al (2017) models since they treat dust polarisation properties

ENDELSE

known_mdls=['MC10','DBP90','DL01','WD01_RV5P5B','DL07','J13','G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD'] 
pol_mdls=['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
test_model = where(pol_mdls eq use_model,ct)
if ct eq 0 then begin
   message,'The only models with polarisation are G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
   stop
end

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

;; ;=== AN EXAMPLE FOR G17_MODELA/B/C/D
;; ;=== Here we fit the dust abundances of the model, the dust
;; ;=== polarisation and a plug-in:
;; ;=== (i) synchrotron emission (also polarized)
;; ;=== The free parameters in the fit are all lower-bounded at zero. The
;; ;=== intensity of the dust-heating radiation field is fixed to 1.0*G0.
pd = [ $
     '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction

     '(*!dustem_params).grains(1).mdust_o_mh',$          ;amorphous carbon mass fraction   
     '(*!dustem_params).grains(2).mdust_o_mh', $         ;amorphous silicate mass fraction

     'dustem_plugin_modify_dust_pol_2',  $               ;This will set the polarization angle
     'dustem_plugin_synchrotron_2', $                    ;Synchrotron amplitude at 10 mm
     'dustem_plugin_synchrotron_3', $                    ;Synchrotron polarization fraction
     'dustem_plugin_synchrotron_4'  $                    ;Synchrotron polarization angle
     ]                   

;--- the SED in example_SED_3.xcat was created with the following
;    model parameters (pd_truth)
pd_truth=[7.8000E-04,7.8000E-04,7.8000E-04, 10. ,1.e-2, 0.3, 45.]
iv=pd_truth+[1.e-4,-1.e-4,1.e-4,-5.,0.,0.1,10] ; initialize DustEMWrap with values shifted from the true solution

Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
llims=replicate(0.,Npar)
fpd=[ '(*!dustem_params).G0']  & fiv=[1.0]                       

Nfix=n_elements(fpd)
if n_elements(fiv) ne Nfix then begin
   message,'Number of fixed parameters (fpd) does not equal number of initial values of fixed parameters (fiv)',/info
   stop
end

;== INITIALISE DUSTEM

;dustem_init,mode=use_model,grain_keywords=['logn-chrg-spin','?','plaw-pol'],polarization=1

dustem_init,model=use_model,polarization=1

!dustem_verbose=1
!dustem_show_plot=1
!dustem_nocatch=1

;=== READ EXAMPLE SED DATA
dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
file=dir+'example_SED_3.xcat'
IF keyword_set(sed_file) THEN file=sed_file
sed=read_xcat(file,/silent)

;== SET THE OBSERVATIONAL STRUCTURE
;== sed is passed twice -- the first occurrence is the SED that you
;== wish to fit, the second occurrence is the SED that you wish to visualise. 
dustem_set_data,sed,sed

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

;== INITIALIZE ANY PLUGINS

dustem_init_plugins,pd,fpd

;== INITIALIZE ANY FIXED PARAMETERS FOR PLUGINS
if Nfix gt 0 then dustem_init_fixed_params,fpd,fiv

;=== RUN fit
tol=1.e-16
xtol=1.e-16
Nitermax=30; just an initialization

xrange=[1.,2.e4]
yrange=[1e-4,1.e3]

loadct,13

t1=systime(0,/sec)

title='FIT POLARISATION EXAMPLE'    

res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []',title=title)
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
dustemwrap_plot,*(*!dustem_fit).CURRENT_PARAM_VALUES,dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (Final fit)'
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_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
;=== 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_and_restore,dustem_spectra_st=dustem_spectra_st
  ;==== plot result taken from the saved fits table
  res=*(*!dustem_fit).CURRENT_PARAM_VALUES
  dustemwrap_plot,res,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  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
  end
ENDIF

message,'Finished dustem_polarisation_example',/info

stop

the_end:

END