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

;This routine is an example of how to fit observational SEDs with
;DustEM and DustEMWrap. The objective is to illustrate how to use DustEMWrap
; (and not to do science -- the fit obtained by running this example is
; likely to be poor!)
;  
;For this example, the code uses the SED in the file example_SED_1.xcat,
;which is distributed in the Data/EXAMPLE_OBSDATA/ directory 
;
;The example SED has Stokes I photometric data points from SPITZER
;IRAC and MIPS and IRAS. Examples illustrating running DustEMWrap to
;fit data with spectral data, polarisation data and extinction data
;are provided in other _readme routines in the src/idl/ directory.

;;  ***COMMENT AH***
;; Remove or update the next few lines -->>
;; No spectrum data points.
;;SPECTRUM data points can be included and the corresponding filter
;;filed must read SPECTRUM. Note that its is note necessary
;;to use the .xcat file format, and data SED can be provided
;;manually, but the observation structure must have the structure shown below.
;;For the example to work, the DustEM and DustEMWrap packages must have
;;been configured and installed succesfully.
;;See dustem_cvs_readme.txt for install instructions).
;; *** END COMMENT AH***
  
;+
; NAME:
;    dustem_fit_polarisation_example  
; PURPOSE:
;    This is an example of how to fit SEDs with 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]
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    Plots, Results save structure
; ACCEPTED KEY-WORDS:
;    model = specifies the interstellar dust mixture used by DustEM
;           'MC10' model from Compiegne et al 2010 (default)
;           '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/SEDs/sample_SED.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
; 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
ENDELSE

dustem_define_la_common

;=== initialise dustem
dustem_init,mode=use_model,kwords=['logn-chrg-spin','plaw-pol','plaw-pol'],/pol

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

filename=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/example_sed3.xcat'
sed=read_xcat(filename)

dustem_set_data,sed

;=== Set which parameters you want to fit
;===== This example adds synchrotron with a polarization angle of 45 degree
pd = [ $
             '(*!dustem_params).G0', $                   ;G0
             '(*!dustem_params).grains(0).mdust_o_mh',$  ;PAH0 mass fraction
             '(*!dustem_params).grains(1).mdust_o_mh',$  ;PAH0 mass fraction   
             '(*!dustem_params).grains(2).mdust_o_mh', $    ;PAH1 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
             ]                   
p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04, 10.,1.e-2, 0.3, 45.]
iv=p_truth+[0.,1.e-4,1.e-4,1.e-4,-5.,0.,0.1,10]      ;initial values shifted from solution
;Final parameters :        1.0000000   0.00075715537   0.00075280676   0.00092480417     0.010006239      0.29981309       45.000163


Npar=n_elements(pd)  
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
llims=replicate(0,Npar)

dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims;,tied=tied
dustem_init_plugins,pd

;=== 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 EXAMPLE'    ;does not work if title is undefined ....
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)

the_end:

END