PRO dustem_stellarpopisrf_example,model=model $
                                 ,sed_file=sed_file $
                                 ,isrf_file=isrf_file $
                                 ,Nitermax=Nitermax $
                                ,postscript=postscript $
                                 ,fits_save=fits_save $
                                 ,wait=wait $
                                 ,noobj=noobj $
                                 ,verbose=verbose $
                                 ,help=help

;+
; NAME:
;    dustem_stellarpopisrf_example  
;
; PURPOSE:
; This routine is an example of how to fit an observational SED
; (StokesI only) with DustEM and DustEMWrap, and an ISRF that is due
; to a user-defined population of nearby main sequence stars.
;
; DustEMWrap reads information about the stellar spectral types (effective
; temperature, radius) from the EEM_dwarf_UBVIJHK_colors_Teff.txt file
; that is located in the Data/STELLARPOPS/ directory. This data file was
; authored by Prof. Erik Mamajek (see the file for more details).
;  
; For this example, we generate an SED using an input model and then
; launch the fit with a starting guess that has been shifted away from
; the true parameter values.
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:
;    dustem_fit_stellarpopisrf_example[,model=][sed_file=][isrf_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_2.xcat' is used.  
;    isrf_file = text file describing ISRF
;    postscript = if set, final plot is saved as postscript file
;    Nitermax = maximum number of fit iterations. Default is 5.
;    fits_save = if set, save the fit results in a binary
;               FITS file. 
;    postscript = if set, final plot is saved as postscript 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_stellarpopisrf_example
;    dustem_fit_stellarpopisrf_example,Nitermax=10,fits_save='/tmp/mysavefile.fits'
;    dustem_fit_stellarpopisrf_example,model='DBP90',isrf_file='./myisrf_habing.dat',sed_file='./mysed.xcat'
;
; 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='DBP90'    ;Default is last dustem model
ENDELSE

exists=dustem_test_model_exists(use_model)
if exists eq 0 then $
   message,'Unknown dust model'

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

;;===============================
;; example parameter initialisation 1
;; Here we attempt to fit the distance to a single O7V star
;; Other ISRF contributions are fixed to ~zero

pd = ['dustem_plugin_stellar_population_O7V3'];          ;distance to O7V star
true_vals =   [20]                ; true distance to the star
iv = true_vals+[4.]               ; distance we use as an initial guess

fpd=['dustem_plugin_stellar_population_O7V4' , $ ;number of O7V star (FIXED)
     '(*!dustem_params).G0'] ; Mathis field (FIXED TO ~ZERO)
fiv=[1.,1.e-12]

Npar=n_elements(pd)
ulimed=replicate(0,Npar) & ulimed[0]=1 
ulims=replicate(0,Npar) & ulims[0]=25.
llimed=replicate(1,Npar)
llims=replicate(1.e-15,Npar) & llims[0]=15.        

;;=== end example part 1
;;===============================

;; ;;===============================
;; ;;=== example parameter initialisation 2
;; ;;=== Here we aim to fit the distance to a single O7V star
;; ;;=== We also include a diffuse Habing ISRF that we generate below
;; ;;=== dustem_create_rfield

;; pd = [ 'dustem_plugin_stellar_population_O7V3' , $          ;distance to O7V star
;;      'dustem_plugin_modify_isrf_1']                         ; amplitude of user-ISRF
;; true_vals =   [20,1.]                ;  distance to the star and true Habing ISRF amplitude
;; iv = true_vals+[4.,-0.2]               ; distance and ISRF amplitude we use as initial guesses

;; fpd=['dustem_plugin_stellar_population_O7V4' , $ ;number of O7V star (FIXED)
;;      '(*!dustem_params).G0'] ; Mathis field (FIXED TO ~ZERO)
;; fiv=[1.,1.e-12]

;; Npar=n_elements(pd)
;; ulimed=replicate(0,Npar) & ulimed[0]=1 
;; ulims=replicate(0,Npar) & ulims[0]=25.
;; llimed=replicate(1,Npar)
;; llims=replicate(1.e-15,Npar) & llims[0]=15.        

;; mywaves=[]
;; my_isrf_file='./myisrf_habing.dat'
;; myisrf=dustem_create_rfield([0],isrf=2,fname=my_isrf_file,x=mywaves)
;; isrf_file=my_isrf_file

;; 
;;=== end example part 2
;; ;;===============================

;=== INITIALIZE DUSTEMWRAP

dustem_init,model=use_model,polarization=use_polarization

!dustem_nocatch=1
!dustem_verbose=1
IF keyword_set(noobj) THEN !dustem_noobj=1

if keyword_set(isrf_file) then begin
   message,'Setting ISRF component from file: '+isrf_file,/info
   !dustem_isrf_file=ptr_new(isrf_file)
endif

;=== GENERATE DATA:
;NB: HERE WE ARE READING AN SED FILE JUST TO SET-UP THE SED STRUCTURE
;AND FILTERS. WE REPLACE THE STOKES I and STOKES I UNCERTAINTY VALUES
;LATER USING THE MODEL ITSELF
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

;=== initializing IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro
for i=4l,n_tags(sed)-1 do begin
    sed.(i) = sed.sigmaii
endfor
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

;=== Generate the emission data using the dust model
sed.StokesI = dustem_compute_sed(true_vals,st=st)
sed.SigmaII = sed.StokesI*0.01

;== SET THE OBSERVATIONAL STRUCTURE
dustem_set_data, m_fit=sed, m_show=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-16

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


tit='FIT WITH STELLAR POP ISRF EXAMPLE'
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='./'
    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
  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

;; moved the following to dustem_fitsio_example.pro
   ;; 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,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 made a plot using the data in this file',/info
     wait,wait
  ENDIF
ENDIF

message,'Finished dustem_stellarpopisrf_example',/info


the_end:

END