PRO dustem_myisrf_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_myisrf_example  
;
; PURPOSE:
; This routine is an example of how to generate an alternative ISRF
; and then fit an observational SED (StokesI only) with DustEM and
; DustEMWrap using this non-default ISRF. The objective is to
; illustrate how to use dustem_plugin_modify_isrf 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 
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:
;    dustem_myisrf_example[,model=][,sed_file=][,isrf_file=][,postscript=][,Nitermax=][,fits_save=][,/help,/wait,/verbose]
;
; INPUTS:
;  
; 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.  
;    isrf_file = text file describing ISRF
;    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 = 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
;    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_myisrf_example,isrf_file='my_isrf.xcat',nitermax=10,fits_save='/tmp/mysavefile.fits'
;
; 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_myisrf_example'
  goto,the_end
END

IF keyword_set(model) THEN BEGIN
  use_model=strupcase(model)
ENDIF ELSE BEGIN
  use_model='DBP90'    ;default dust model for this example
ENDELSE

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

use_verbose=0
use_polarization=0   ; initialize Dustemwrap in no polarization mode 
use_window=2          ; default graphics window number to use for plotting the results
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

;=== AN EXAMPLE FOR DBP90
;=== Here we fit the PAH abundance of the model, and the amplitude
;=== of a user-supplied ISRF. A synchrotron component is included
;=== The other grain parameters are fixed to their default values.
;=== The free parameters are all lower-bounded at zero.
;=== Note that if '(*!dustem_params).G0' is not set explicitly as a fixed parameter, then
;=== the ISRF used by DustEM will include a standard Mathis field with
;=== amplitude, G0=1.
use_model='DBP90' ; you should specify this above, or in the command line
pd = [ $
     '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
     'dustem_plugin_modify_isrf_1', $           ; Amplitude of user-defined ISRF
     'dustem_plugin_synchrotron_2']             ;Synchrotron amplitude at 10 mm
iv =   [5.4e-4, 0.1, 0.01]
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
llims=replicate(1.e-15,Npar)
; To suppress the Mathis field altogether
;fpd=['(*!dustem_params).G0']
;fiv=[1.e-12]


;; ;=== AN EXAMPLE FOR DL07 
;; ;=== Here we fit the PAH abundances of the model, and the amplitude
;; ;=== of a user-supplied ISRF. A synchrotron component is included
;; ;=== The other grain parameters are fixed to their default values.
;; ;=== The free parameters are all lower-bounded at zero.
;; use_model='DL07' ; 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_plugin_modify_isrf_1', $                   ; Ampitude of user-defined ISRF
;;      'dustem_plugin_synchrotron_2']                     ;Synchrotron amplitude at 10 mm
;; iv =   [5.4e-4, 5.4e-4, 0.1, 0.01]
;; Npar=n_elements(pd)
;; ulimed=replicate(0,Npar)
;; llimed=replicate(1,Npar)
;; llims=replicate(1.e-15,Npar)
;; ; To suppress the Mathis field altogether
;; ;fpd=['(*!dustem_params).G0']
;; ;fiv=[1.e-12]


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,polarization=use_polarization
!dustem_nocatch=1
!dustem_verbose=use_verbose
IF keyword_set(noobj) THEN !dustem_noobj=1

;== Fill a structure with default inputs to the model. This
;== includes the default ISRF, which we want here just for plotting a
;comparison to the ISRF that we construct.
dir_in=!dustem_soft_dir  
st_model=dustem_read_all(dir_in)

;=== CREATE A DIFFERENT ISRF USING DUSTEM_CREATE_RFIELD AND OUTPUT TO FILE
if not keyword_set(isrf_file) then begin
   mywaves=[]
   my_isrf_file='./myisrf_3bb_8000K.dat'
   myisrf=dustem_create_rfield([8000,8000,8000],wdil=[1.d-11,1.d-11,1.d-11],isrf=0,fname=my_isrf_file,x=mywaves)
;   myisrf=dustem_create_rfield([8000],wdil=[1.d-12],isrf=0,fname=my_isrf_file,x=mywaves)
;myisrf=dustem_create_rfield([3500,4000,5500,6500,8000],wdil=[1.d-13,1.d-13,1.d-13,1.d-13,1.d-13],isrf=1,fname='./myisrf_mathis_5bb.dat',x=mywaves)
;myisrf=dustem_create_rfield([3500,4000,5500,6500,8000],wdil=[1.d-13,1.d-13,1.d-13,1.d-13,1.d-13],isrf=2,fname='./myisrf_habing_5bb.dat',x=mywaves)
;myisrf=dustem_create_rfield([3500,4000,5500,6500,8000],wdil=[1.d-13,1.d-13,1.d-13,1.d-13,1.d-13],isrf=0,fname='./myisrf_5bb.dat',x=mywaves)
   isrf_file=my_isrf_file
endif


;=== PLOT AS A SANITY CHECK
xtit=textoidl('\lambda (\mum)')
ytit=textoidl('log ISRF [4 \pi I_\nu (erg/cm^2/s/Hz)]')
tit='DUSTEM MYISRF EXAMPLE'
yr=[min(st_model.isrf.isrf) < min(myisrf),5.*(max(st_model.isrf.isrf)>max(myisrf))]
xr=[0.001,1e6]
window,use_window,xs=600,ys=400,tit=tit & use_window=use_window+2
cgplot,st_model.isrf.lambisrf,st_model.isrf.isrf $
       ,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog $
       ,title=tit,xtit=xtit,ytit=ytit,color=cgcolor('black'),/nodata
oplot,st_model.isrf.lambisrf,st_model.isrf.isrf,col=cgcolor('black'),thick=2
oplot,mywaves,myisrf,col=cgcolor('red'),thick=2
al_legend,/top,/right,clear=0,box=0 $
          ,['default ISRF','my ISRF'] $
          ,linsize=0.5,lines=0 $
          ,color=[cgcolor('black'),cgcolor('red')] $
          ,thick=2,charsize=1.3

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

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


;== 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
; ;== ADJUSTED DURING THE FIT
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

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

xr = [1.,5e5]
yr = [5e-8,1.00e6]
tit='MYISRF EXAMPLE'
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
xtit=textoidl('\lambda (\mum)')

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 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,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'
  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) 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 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,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,/info
     wait,wait
  ENDIF
ENDIF

message,'Finished dustem_myisrf_example',/info


the_end:

END