PRO dustem_fit_intensity_mbb_example,model=model $

                                    ,sed_file=sed_file $
                                    ,Nitermax=Nitermax $
                                    ,postscript=postscript $

                                    ,fits_save=fits_save $

                                    ,noobj=noobj $

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

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

; (StokesI only) with DustEM and DustEMWrap.

;  
; For this example, the code uses the SED in the file example_SED_1.xcat,
; which is distributed in the Data/EXAMPLE_OBSDATA/ directory 
;

; In addition, this example uses the plugin dustem_plugin_mbbdy.pro to
; fit the FIR part of the spectrum using a modified black-body rather
; than a physical ISM dust model for the big grain population.
;  

; 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_mbb_example[,model=][sed_file=][,postscript=][,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.  

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

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

;    fits_save = if set, save the fit 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_mbb_example

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

;    dustem_fit_intensity_mbb_example,model='DBP90'
;
; MODIFICATION HISTORY:

;    Written by AH 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_intensity_mbb_example'
  goto,the_end
END

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

  use_model='MC10'    ;Default is the MC10 model

ENDELSE

exists=dustem_test_model_exists(use_model)
if exists ne 1 then $
   message,'Unknown dust model'

use_polarization=0   ; initialize Dustemwrap in no polarization mode 
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

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

;=== Here we fit the

;=== intensity of the dust-heating radiation field as well as several plug-ins:

;=== (i) a modified black body  (initial values correspond to a
;          blackbody peaking in the far-infrared)

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

;=== We invoke the MC10 dust model, and fit the abundances of the PAH components. The BG 
;=== components are fixed to ~zero since we want to fit this part of the
;=== observed SED by a MBB.
;use_model='MC10' ; you should specify this when you run the routine!

pd = [ $
     '(*!dustem_params).G0', $                           ;G0

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

     'dustem_plugin_mbbdy_1', $                          ;Amplitude of mBB
     'dustem_plugin_mbbdy_2', $                          ;T of mBB

     'dustem_plugin_mbbdy_3']                            ;Emissivity index of mBB

iv =   [1.8, 7.8e-4, 7.8e-4, 4.e-5, 22., 1.8]

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

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

fpd = [ $

     '(*!dustem_params).grains(2).mdust_o_mh',$          ;amCBEx
     '(*!dustem_params).grains(3).mdust_o_mh',$          ;amCBEx
     '(*!dustem_params).grains(4).mdust_o_mh']           ;aSilx

fiv =   [1.e-12,1.e-12,1.e-12]

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,/show_plots

!dustem_nocatch=1
!dustem_verbose=use_verbose

IF keyword_set(noobj) THEN !dustem_noobj=1

;=== READ EXAMPLE SED DATA
dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'

file=dir+'example_SED_4.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 (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 -- 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

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, 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='MBB FIT 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

ENDIF

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

  set_plot,'X'

  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

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 following to the dustem_fitsio_example
   ;; 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 !dustem_noobj THEN BEGIN 
  ;;   dustemwrap_plot_noobj,res,st=dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  ;; ENDIF ELSE BEGIN
  ;;   dustemwrap_plot,res,st=dustem_spectra_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,/info

     wait,wait

  ENDIF

ENDIF

the_end:

message,'Finished dustem_fit_intensity_mbb_example',/info

END