PRO dustem_myisrf_example,model=model $
                         ,sed_file=sed_file $
                         ,Nitermax=Nitermax $
                         ,postscript=postscript $
                         ,fits_save_and_restore=fits_save_and_restore $
                         ,wait=wait $
                         ,verbose=verbose $
                         ,help=help

;+
; NAME:
;    dustem_fit_myisrf_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 the modify_isrf plugin 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_2.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_myisrf_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 structure in binary FITS table format
;
; ACCEPTED KEY-WORDS:
;    model = specifies the interstellar dust mixture used by DustEM
;           'MC10' model from Compiegne et al 2010 
;           '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
;    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_2.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, save the fit results in a binary
;               FITS file. The code then restore this file and plots
;               the results using the saved results information.
;    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_stellarpopisrf_example
;    dustem_fit_stellarpopisrf_example,Nitermax=1,fits_save_and_restore='/tmp/mysavefile.fits'
;    dustem_fit_stellarpopisrf_example,model='DBP90'
;
; MODIFICATION HISTORY:
;    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_stellarpopisrf_example'
  goto,the_end
END

IF keyword_set(model) THEN BEGIN
  use_model=strupcase(model)
ENDIF ELSE BEGIN
  use_model='DL07'    ;Default is last dustem model
ENDELSE


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
;=== Refer to the DustEM and DustEMWrap userguides for an explanation
;    of the different grain types

CASE use_model OF
      'DBP90':BEGIN
         pd = [ $
             '(*!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', $    ;amCBEx
             'dustem_plugin_modify_isrf_1']  ;amplitude of user-defined ISRF (here G0 because mathis is chosen)
        rv =   [4.3e-4, 4.7e-4,6.4e-3,100.] ;rv because we're generating fake data and try to retrieve the original values
        
        iv = rv+[5.00E-4,5.00E-4,5.00E-4,-30] ;###setting the initial parameter vector
        
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)         
        
        
;         fpd=[ $
;           'dustem_plugin_stellar_population_O7V4',$ ;number of stars of stellar population' ,  $     ;multiplicative factor to total ISRF
;           'dustem_plugin_stellar_population_O5V3'$ ; stellar distance
;             ]
;         ;initial parameter values for fixed parameters
;         fiv=[ $
;           2. , $        ;number of stars of pop O7
;           45$      ;dustance of pop O5
;             ]

     END
      'DL01':BEGIN
         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
             'dustem_plugin_modify_isrf_1'] ;Amplitude of user-defined ISRF
             
        rv =   [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,50]
        iv =   rv+[5.e-4, 5.e-4,5.e-5,5.e-3,5.e-3,100]
        
;         fpd=[ $
;             'dustem_plugin_stellar_population_O5V3'$ ; stellar distance
;             ]
;         ;initial parameter values for fixed parameters
;         fiv=[ $
;           45$      ;fixed stellar distance
;             ]

        
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
      END
      'WD01_RV5p5B':BEGIN
;;  ***COMMENT AH***
;; we need to implement this, or remove         
         message, 'WD01 model not yet implemented in DustEMWrap',/info
;;  ***END COMMENT AH***
       END
      'DL07':BEGIN ;takes at least 16 iterations to find minima.
          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
              'dustem_plugin_modify_isrf_1',   $ ;amplitude of user-defined ISRF
              'dustem_plugin_stellar_population_O5V3'   $ ;distance to stellar pop O5
              ] 
                     
        rv =   [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,5.,1.];,10,1.,10.,1.]
        iv =   rv+ [5.e-4, 5.e-4,5.e-4,5.e-3,5.e-3,5.,2]
        
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        
        fpd=[ $
            'dustem_plugin_stellar_population_O5V4',$ ; stellar number of pop O5V
            '(*!dustem_params).G0'                  $
            ]
        ;initial parameter values for fixed parameters
        fiv=[ $
          30., $      ;fixed stellar population number
          10.  $
            ]
        
      END
      'MC10':BEGIN
        
        pd = [ $
              'dustem_plugin_modify_isrf_1', $  ;amplitude of user-defined ISRF
              '(*!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
              ]
        ;initial parameter values for parameters to be fitted   
        rv =   [ $  
          50, $    
          7.8e-4,  $ ;mass fraction of  PAH0  
          7.8e-4,  $ ;mass fraction of  PAH1  
          1.65e-4, $ ;mass fraction of  amCBEx
          1.45e-3, $ ;mass fraction of  amCBEx
          7.8e-3   $ ;mass fraction of aSilx
          ]
          
          iv = rv + [-30,5.E-4,5.E-4,5.E-4,5.E-3,5.E-3]
          
;           fpd = ['dustem_plugin_stellar_population_O5V3'] ;distance to stellar population O5V
;           fiv  = [25]
              
     
        Npar=n_elements(pd)
        
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        
      END
     'J13':BEGIN
        pd = [ $
             '(*!dustem_params).grains(1).mdust_o_mh',$     ;CM20 mass fraction
             '(*!dustem_params).grains(0).mdust_o_mh',$     ;CM20 mass fraction
             '(*!dustem_params).grains(2).mdust_o_mh', $    ;aPyM5
             '(*!dustem_params).grains(3).mdust_o_mh', $    ;aPyM5
             'dustem_plugin_modify_isrf_1']                 ;amplitude of user-defined ISRF
                           
        rv =   [7.8e-4, 7.8e-4,1.65e-4,1.45e-3,35]
        iv = rv+ [5.e-4,5.e-4,5.e-4,5.e-3,35]
        
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        
;         fpd = ['dustem_plugin_stellar_population_O5V3'] ;distance to stellar population O5V
;         fiv  = [5]
      
      END
      'ELSE':BEGIN
         message,'model '+model+' unknown',/continue
         message,'Known models are MC10,DBP90,DL01,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
         stop
      END
ENDCASE


if keyword_set(wait) then begin
   message,'Finished setting dust model and plug-in parameters: '+use_model,/info
   wait,wait
end


dustem_init,model=use_model

!dustem_nocatch=1
!dustem_verbose=1
;!dustem_show_plot=1
!EXCEPT=2 ;so I can locate the plotting error...
;!dustem_dim=1 ; this option is to dim the stellar population ISRF with the current Dustem extinction
;We're fitting total optical depths so the dimmed-ISRF scenario is a reasonable assumption.
;Because of this the final parameter values are slightly lower than the initial (real) values.
;But this is because the first run was not extinct with the current dustem extinction 
;since it comes prior to the dustem run. Meaning the first ISRF has to be extinct to see if the paramters are retrieved.

!dustem_isrf_file=ptr_new(!dustem_soft_dir+'data/ISRF_MATHIS.DAT')

;=== READ EXAMPLE DATA: EXTINCTION

;NB:HERE WE ARE NOT FITTING THE DATA. WE'RE JUST USING THE EXISTING FILTERS.

dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/'
file=dir+'example_SED_2.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


;=== GENERATE EXAMPLE DATA: EXTINCTION

;=== initializing IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro
;sed.StokesI=1.
;setting the errors (well equating all the tags to the sigmaii array but this is an initialization)
for i=4l,n_tags(sed)-1 do begin
    sed.(i) = sed.sigmaii
endfor

;###first call to dustem_set_data so that the !dustem_data tags that are used in the 'compute_' procedures are set

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
; 
; ;== INITIALIZE ANY PLUGINS
; dustem_init_plugins, pd,fpd=fpd 
; 
; dustem_init_fixed_params,fpd,fiv

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

;Generating the fake emission data
sed.StokesI = dustem_compute_sed(rv,st=st)

;== SET THE OBSERVATIONAL STRUCTURE
dustem_set_data, sed,sed


if keyword_set(wait) then begin
   message,'Finished initializing DustEMWrap, including plugins and fixed parameters',/info
   wait,wait
end

;== RUN THE FIT
tol=1.e-18
use_Nitermax=5  ;maximum number of iterations. 
IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax

xr_m = [1.,5e5]
yr_m = [5e-8,1.00e6]


tit='Spectral Energy Distribution'
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
xtit=textoidl('\lambda (\mum)')
;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 
t1=systime(0,/sec)
res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
                      ,/xlog,/ylog,xr_m=xr_m,yr_m=yr_m,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 keyword_set(postscript) THEN BEGIN
    dir_ps='./'
    set_plot,'PS'
    ps_file=dir_ps+postscript
    device,filename=ps_file,/color
ENDIF

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

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
  ;stop
  dustem_read_fits_table,filename=fits_save_and_restore,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,dustem_spectra_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)'
  ENDIF ELSE BEGIN
    dustemwrap_plot,res,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_and_restore+', and made a plot using the data in this file',/info
     wait,wait
  ENDIF
ENDIF

message,'Finished dustem_fit_intensity_example',/info


the_end:

END