PRO dustem_fit_intensity_example,model=model $
                                ,sed_file=sed_file $
                                ,Nitermax=Nitermax $

                                ,fits_save=fits_save $

                                ,help=help $
                                ,wait=wait $

                                ,noobj=noobj $

                                ,verbose=verbose

;+
; NAME:
;    dustem_fit_intensity_example  

;

; PURPOSE:
; This routine is an example of how to fit an observational SED
; (StokesI only) 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
; IRAC, MIPS and IRAS. Examples illustrating running DustEMWrap to
; fit spectral data, polarisation data and extinction data
; are provided in other _example routines in the src/idl/
; directory. See the DustEMWrap User Guide for more information.

;

; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example

;

; CALLING SEQUENCE:

;    dustem_fit_intensity_example[,model=][sed_file=][,Nitermax=][,fits_save=][,/help,/wait,/verbose]

;

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

;    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_1.xcat' is used.  

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

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

;               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

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

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

;    dustem_fit_intensity_example,model='DBP90'

;

; 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_intensity_example'
  goto,the_end
END

;; ;===================================
;; ;=== SET THE DUST MODEL HERE
;; ;===================================

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

  use_model='DBP90'    ;Example with default keywords uses the DBP90 model

ENDELSE

use_polarization=0   ; here we initialize DustEMWrap in no polarization mode since we are only fitting Stokes I
exists=dustem_test_model_exists(use_model,/silent)

if exists ne 1 then $
   message,'Unknown dust model'

;; ;===================================
;; ;=== PARSE OTHER INPUTS AND SET SOME BASIC RUN PARAMETERS
;; ;===================================

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

use_window=0          ; default graphics window number to use for plotting the results

dustem_define_la_common

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

;=== Examples are provided for some of the dust models. 
;=== To try them, uncomment the model that you want to try and re-run

;; ;===================================
;; ;=== AN EXAMPLE FOR DBP90
;; ;===================================

;=== Here we fit the dust abundances of the DBP90 model, the

;=== intensity of the dust-heating radiation field. We also invoke a plug-in
;=== for that fits parameters of synchrotron emission

;=== The free parameters are all lower-bounded at zero.

;=== use_model='DBP90' ; you should specify this above (line 93)

;=== pd is the structure of free parameters

pd = [ $
     '(*!dustem_params).G0', $                           ;G0
     '(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
     '(*!dustem_params).grains(1).mdust_o_mh',$          ;VSG mass fraction

     '(*!dustem_params).grains(2).mdust_o_mh'];,$        ;BG mass fraction
;     'dustem_plugin_synchrotron_1', $                   ;Spectral index of CREs
;     'dustem_plugin_synchrotron_2']                     ;Synchrotron amplitude at 10 mm

;=== iv is the vector of initial values for the free parameters
;=== following lines are for a run without fitting the synchrotron
;=== i.e. [ initial values for G0,PAH0,VSG,BG]

iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3] 

; initial values vector for run including the synchrotron plugin
; [G0,PAH0,VSG,BG,alpha_CR,Amp_syn]
;iv =   [1.6, 2.2e-4, 5.7e-4, 3.4e-3, 2.7,0.01]

; initial values vector for run including the synchrotron plugin and

; fixing the G0 and BB continuum parameters (see below)
; [PAH0,VSG,BG,alpha_CR,Amp_syn]

;iv =   [2.2e-4, 5.7e-4, 3.4e-3, 2.7,0.01]

Npar=n_elements(pd)

ulimed=replicate(0,Npar) ; flag ON=1, OFF=0
llimed=replicate(1,Npar) ; flag ON=1, OFF=0
llims=replicate(1.e-15,Npar) ; lower limit value for each free parameter
; ulims=replicate(1.e15,Npar) ; upper limit value for each free parameter

fpd=[] & fiv=[]

; example using fixed parameters (ISRF and a NIR continuum)

; Uncomment the following lines to include these fixed parameters
;fpd = ['(*!dustem_params).G0', $    ;G0
;     'dustem_plugin_continuum_1', $                    ;Temperature of a BB
;     'dustem_plugin_continuum_2']                      ;Peak amplitude of a BB

; fixed values vector 
; [G0,T_BB,Amp_BB]
;fiv=[0.5,750,1.e-2]

;; ;==================================

;; ;=== EXAMPLES FOR OTHER PHYSICAL DUST MODELS

;; ;=== START BELOW HERE
;; ;==================================

;; ;===================================

;; ;=== AN EXAMPLE FOR DL07

;; ;===================================

;; ;=== Here we fit the dust abundances of the model and the
;; ;=== intensity of the ISRF (via gas.G0 since the model
;; ;=== includes spinning dust)

;; ;=== The free parameters are all lower-bounded at zero.

;; ;=== use_model='DL07' ; you should specify this above (line 93), or in the command line when you run the example
;; pd = [  '(*!dustem_params).gas.G0'], $                 ;G0

;;         '(*!dustem_params).grains(0).mdust_o_mh',$     ;PAH0 mass fraction
;;         '(*!dustem_params).grains(1).mdust_o_mh',$     ;PAH1 mass fraction
;;         '(*!dustem_params).grains(2).mdust_o_mh', $    ;Gra
;;         '(*!dustem_params).grains(3).mdust_o_mh', $    ;Gra

;;         '(*!dustem_params).grains(4).mdust_o_mh']      ;aSil

;; iv =   [1.5 ,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3]
;; 
;; ;=== Uncomment the following lines if you want to fix the ISRF instead of leaving it as a free parameter
;; ;pd = ['(*!dustem_params).grains(0).mdust_o_mh',$     ;PAH0 mass fraction            
;; ;        '(*!dustem_params).grains(1).mdust_o_mh',$     ;PAH1 mass fraction
;; ;        '(*!dustem_params).grains(2).mdust_o_mh', $    ;Gra
;; ;        '(*!dustem_params).grains(3).mdust_o_mh', $    ;Gra
;; ;        '(*!dustem_params).grains(4).mdust_o_mh']      ;aSil
;; ;iv =   [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3]
;; ;fpd=['(*!dustem_params).gas.G0']
;; ;fiv=[2.]
;; 
;; Npar=n_elements(pd)
;; ulimed=replicate(0,Npar)
;; llimed=replicate(1,Npar)
;; llims=replicate(1.e-15,Npar)

;; ;===================================
;; ;=== AN EXAMPLE FOR MC10
;; ;===================================

;; ;; ;=== Here we fit the dust abundances of the MC10 model, the
;; ;; ;=== intensity of the dust-heating radiation field as well as a plug-in:
;; ;; ;===  (i) continuum due to a blackbody 
;; ;; ;=== The intensity of the dust-heating radiation field is fixed to
;; ;; ;=== 1.5*G0 and the tmperature of the blackbody is fixed to 1200K
;; ;; ;=== The free parameters in the fit are lower-bounded at zero.

;; ;; ;=== use_model='MC10' ; you should specify this above, or in the command line

;; pd = [ $
;;      '(*!dustem_params).grains(0).mdust_o_mh'$      ;PAH0 mass fraction
;;      ,'(*!dustem_params).grains(1).mdust_o_mh' $    ;PAH1 mass fraction
;;      ,'(*!dustem_params).grains(2).mdust_o_mh' $    ;amCBEx
;;      ,'(*!dustem_params).grains(3).mdust_o_mh' $    ;amCBEx
;;      ,'(*!dustem_params).grains(4).mdust_o_mh' $    ;aSilx
;;      ,'dustem_plugin_continuum_2' $                 ;Intensity at peak of the continuum
;;      ]
;; iv =   [ 7.8e-4, 7.8e-4, 1.65e-4, 1.45e-3, 6.7e-3, 0.003]
;; Npar=n_elements(pd)
;; ulimed=replicate(0,Npar)
;; llimed=replicate(1,Npar)
;; llims=replicate(1.e-15,Npar)
;; fpd=[ $
;;     '(*!dustem_params).G0'   $    ; ISRF intensity 
;;     ,'dustem_plugin_continuum_1' $    ;temperature of blackbody the produces the continuum 
;;     ]
;; fiv=[1.5, 1200.]

;; ;===================================

;; ;=== 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 the free 
;; ;=== parameters of 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.

;; ;=== use_model='J13' ; you should specify this above, or in the command line

;; 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_freefree_1' $                              ;ionized gas temperature
;;      ,'dustem_plugin_freefree_2' $                              ;free-free amplitude

;;      ,'dustem_plugin_continuum_2']                              ;intensity at peak of the BB continuum

;; iv =   [1.2,1.7e-3, 6.3e-4, 2.55e-3, 2.55e-3, 7500., 0.4, 0.001]         

;; fpd=[ 'dustem_plugin_continuum_1'] ; Temperature of the BB
;; fiv=[1000.]                       

;; Npar=n_elements(pd)

;; ulimed=replicate(0,Npar)
;; llimed=replicate(1,Npar)

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

;;===================================
;;=== 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 EXAMPLE SED DATA
;;===================================

dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'

file=dir+'example_SED_1.xcat'

IF keyword_set(sed_file) THEN file=sed_file
sed=read_xcat(file,/silent)

if keyword_set(wait) then begin
   message,'Finished reading SED data: '+file,/info
   wait,wait
end

;;===================================

;;=== ADJUST THE UNCERTAINTIES FROM WITHIN THE CODE (FOR ILLUSTRATION)

;;===================================

ind=where(sed.sigmaII LT (0.2*sed.StokesI)^2,count)

IF count NE 0 THEN sed[ind].sigmaII=(0.2*sed[ind].StokesI)^2

;;===================================
;=== SET THE OBSERVATIONAL STRUCTURE
;;===================================
;== sed is passed twice in the call -- the first occurrence (m_sed) is the SED that you
;== wish to fit, the second occurrence (m_show) is the SED that you wish to visualise. 

dustem_set_data,m_fit=sed,m_show=sed

;;===================================
;=== INITIALISE DUSTEM WITH INITIAL VALUES AND LIMITS OF ALL PARAMETERS 
;;===================================

dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims

if keyword_set(wait) then begin

   message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info

   wait,wait
end

;;===================================
;;=== INFORMATION TO RUN THE FIT
;;===================================

tol=1.e-16     ;fit tolerence

;;===================================
;;=== INFORMATION TO MAKE THE PLOT
;;===================================

yr=[1.00e-4,1.00E2] ; y-axis limits
xr=[1.00E0,6.00e4] ; x-axis limits

tit='FIT INTENSITY EXAMPLE' ; plot title

ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title
xtit=textoidl('\lambda (\mum)') ; x-axis title

;;===================================
;;===  RUN THE FIT
;;===================================

t1=systime(0,/sec)
res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
                      ,/xlog,/ylog,xr=xr,yr=yr,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 !dustem_noobj THEN BEGIN

  dustem_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(wait) then begin
   message,'Made the plot of the final results',/info
   wait,wait
end

;;===================================
;;=== WRITE OUT THE DUSTEMWRAP FITTING RESULTS IN A BINARY FITS FILE
;;===================================

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 fitting results and remake any plots has been saved in the FITS table
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_intensity_example',/info

END