PRO dustem_fit_polarization_example,model=model $
                                    ,sed_file=sed_file $
                                    ,Nitermax=Nitermax $
                                    ,postscript=postscript $
                                    ,fits_save=fits_save $
                                    ,help=help $
                                    ,noobj=noobj $
                                    ,wait=wait $
                                    ,verbose=verbose
  ;+
; NAME:
;    dustem_fit_polarization_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_polarization_example[,model=][sed_file=][,postscript=][,Nitermax=][,fits_save=][,/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. See userguide or dustem_test_model_exists.pro
;            for more details about available models in current release.
;    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 file
;    Nitermax = maximum number of fit iterations. Default is 5.
;    fits_save = if set, saves results in a FITS file.
;    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:
;
; EXAMPLES
;    dustem_fit_polarization_example
;    dustem_fit_polarization_example,Nitermax=1,fits_save='/tmp/mysavefile.fits'
;    dustem_fit_polarization_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_polarization_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 polarization properties
ENDELSE

exists=dustem_test_model_exists(use_model,/pol)
if exists ne 1 then $
   message,'Unknown/non-polarized dust model: '+use_model

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
;; ;=== polarization and a plug-in:
;; ;=== (i) synchrotron emission (also polarized)
;; ;=== The free parameters in the fit are lower-bounded at
;; ;=== ~zero. The synchrotron parameters are constrained betweem
;; ;=== upper and lower bounds
;; ;=== The intensity of the dust-heating radiation field and the
;; ;=== spectral index of the CRE energy spectrum are fixed.

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',  $               ;Dust 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
     ]                   
iv=[7.2e-4,7.2e-4,8.5e-4,-10.,0.015,0.25,35.]

Npar=n_elements(pd)
ulimed=[0,0,0,0,1,1,1]
llimed=[1,1,1,0,1,1,1]
ulims=[0,0,0,0,0.02,0.5,60.]
llims=[1.e-15,1.e-15,1.e-15,0,0.005,0.05,30.]

fpd=[ '(*!dustem_params).G0' , $ ; intensity of the dust-heating radiation field
      'dustem_plugin_synchrotron_1']  ;Synchrotron CR spectral index
fiv=[1.0,3.0]                       
use_polarization=1

;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,pol=use_polarization
!dustem_verbose=1
;!dustem_show_plot=1
!dustem_nocatch=1
IF keyword_set(noobj) THEN !dustem_noobj=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,m_fit=sed,m_show=sed



;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE
;== ADJUSTED DURING THE FIT
dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims,pol=use_polarization

;=== RUN fit
;number of iterations has already been specified in the beginning of the code.
tol=1.e-16
xtol=1.e-16
use_Nitermax=5; an initialization
IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax

xrange=[1.,2.e4]
yrange=[1e-4,1.e3]
;Set show_plot to 0 to hide plot
;Commented or set to 1 is the same since !dustem_show_plot (existing sysvar) is initialized to 1 in dustem_init
;show_plot = 0 


loadct,13
t1=systime(0,/sec)
title='FIT POLARIZATION EXAMPLE'
res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=use_Nitermax,xrange=xrange,/xstyle,yrange=yrange $
                      ,/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []' $
                      ,title=title,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

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

IF keyword_set(postscript) THEN BEGIN
  device,/close
  set_plot,mydevice
  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,title=title+' (From saved FITS file)'
  ;; ENDIF ELSE BEGIN
  ;;   dustemwrap_plot,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=title+' (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
  end
ENDIF


the_end:
message,'Finished dustem_polarization_example',/info

END