PRO dustem_fit_sed_readme,postcript=postcript,model=model,help=help,png=png,itermax=itermax

;This Readme describes how to fit SEDs with dustem
;It runs on the SED stored into the file sample_SED.xcat
;in this directory. Remember that the goal here is not necessarily to
;obtain a good fit in the end, but to illustrate the method.
;The provided SED has only photometric data points from SPITZER
;IRAC and MIPS and IRAS. No spectrum data points.
;SPECTRUM data points can be included and the corresponding filter
;filed must read SPECTRUM. Note that its is note necessary
;to use the .xcat file format, and data SED can be provided
;manually, but the observation structure must have the structure shown below.
;Obviously, the dustem package must have been installed succesfully (see
;dustem_cvs_readme.txt for install instructions).

;+
; NAME:
;    dustem_fit_sed_readme
; PURPOSE:
;    This is an example of how to fit SEDs with the dustem wrapper.
;    It is meant to be an example to follow when writing your own
;    programs using the dustem IDL wrapper.
;    It is not meant to reproduce the result in Compiegne et al 2010
;    The SED used here is in sample_SED.xcat
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:

;    dustem_fit_sed_readme,postcript=postcript,model=model,help=help

; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:

;    model = Selects one of the dust mixture used by dustem

;           'COMPIEGNE_ETAL2010' from Compiegne et al 2010 (default)
;           'DBP90' from Desert et al 1990
;           'DL01' from Draine & Li 2001
;           'DL07' from Draine & Li 2007
;    postcript = if set plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps
;    help      = If set print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem fortran code must be installed
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES

;    dustem_fit_sed_readme,model='COMPIEGNE_ETAL2010'

; MODIFICATION HISTORY:
;    Written by J.P. Bernard April 1st 2011
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_fit_sed_readme'
  goto,the_end
ENDIF

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

ENDIF ELSE BEGIN

  use_model='THEMIS';G17_MODELA';'COMPIEGNE_ETAL2010'    ;Default is last dustem model

ENDELSE

IF keyword_set(png) THEN BEGIN

  dir_png=!dustem_dat+'/Figures/'

  force_mkdir,dir_png
ENDIF

polarization=0.   ;default is no polarization in models

;dustem_init,model=use_model,pol=polarization

;=== Set which parameters you want to fit

CASE use_model OF

      'DBP90':BEGIN
         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_params).grains(2).mdust_o_mh', $    ;amCBEx

             'dustem_plugin_continuum_2']                ;Intensity of NIR continuum

        iv =   [1.0, 4.3e-4, 4.7e-4,6.4e-3,0.001]

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

     END
      'DL01':BEGIN
         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_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_continuum_2']                ;Intensity of NIR continuum

        iv =   [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001]

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

      END
      'DL07':BEGIN

          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_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_continuum_2']                ;Intensity of NIR continuum

        iv =   [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001];,10,1.,10.,1.]

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

        ;llims =[1e-9  ,0.   ,0.   ,0.    ,0.   ,0.   ,0.   ,1.  ,0.   , 1.   ,0.   ]

      END
      'COMPIEGNE_ETAL2010':BEGIN

        
        ;ORIGINAL LINES
        

        ;parameter description of parameters to be fitted

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

              'dustem_plugin_continuum_2', $      ;intensity of NIR continuum
              'dustem_plugin_synchrotron_2',$

              '(*!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' $
              ]
        ;initial parameter values for parameters to be fitted   
        iv =   [ $  
          1.0, $

          0.002,   $ ;intensity of NIR continuum

          0.01,$    

          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  
          ]

        Npar=n_elements(pd)

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

        llims=replicate(0.,Npar)

        ;=== Fixed parameters

        ;parameter description of parameters to be set to a non-default value

        ;fpd=[ $

          ;'!dustem_isrf_star_add[0].amplitude', $        ;amplitiude of first stellar contribution to ISRF
          ;'!dustem_isrf_star_add[0].distance', $        ;amplitiude of first stellar contribution to ISRF

          ;'(*!dustem_params).gas.G0' $     ;multiplicative factor to total ISRF

          ;'dustem_create_continuum_2' $     ;intensity of NIR continuum

          ;  ]

        ;initial parameter values for fixed parameters

        ;fiv=[ $

          ;1.0, $      ;amplitiude of first stellar contribution to ISRF
          ;1.0, $      ;distance of first stellar contribution to ISRF

        ;  1. $        ;multiplicative factor to total ISRF

          ;1.e-6 $      ;intensity of NIR continuum

         ;   ]

      END

     'AJ13':BEGIN
        pd = [ $
             '(*!dustem_params).G0', $      ;G0
             '(*!dustem_params).grains(1).mdust_o_mh',$     ;PAH1 mass fraction
             '(*!dustem_params).grains(0).mdust_o_mh',$  ;PAH0 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', $      ;aSil

             'dustem_plugin_continuum_2']                ;Intensity of NIR continuum

        iv =   [1.0, 7.8e-4, 7.8e-4,1.65e-4,1.45e-3,7.8e-3,0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)

      END

      'G17_MODELA':BEGIN

         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_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_continuum_2']                ;Intensity of NIR continuum

        iv =   [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)

        polarization=1

      END
      'G17_MODELB':BEGIN
         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_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_continuum_2']                ;Intensity of NIR continuum

        iv =   [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001]
        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        polarization=1
      END
      'THEMIS':BEGIN
         pd = [ $

             '(*!dustem_params).G0', $   ;G0

             '(*!dustem_params).grains(0).mdust_o_mh',$  
             '(*!dustem_params).grains(1).mdust_o_mh',$     
             '(*!dustem_params).grains(2).mdust_o_mh', $  

             'dustem_plugin_continuum_2']                ;Intensity of NIR continuum

        iv =   [1.0,5.4e-4,2.33e-3,8.27e-3,0.001]

        Npar=n_elements(pd)
        ulimed=replicate(0,Npar)
        llimed=replicate(1,Npar)
        llims=replicate(0.,Npar)
        polarization=1
      END

ENDCASE

dustem_init,model=use_model,kwords=['plaw-ed','logn','logn'];,pol=polarization

!dustem_nocatch=1

!dustem_verbose=1
!dustem_show_plot=1

;=== Read sample SED
;=== Composite SED from Compiegne et al 2010, gathered by C. Bot
dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
file=dir+'Gal_composite_spectrum.xcat'

;file=dir+'Gal_composite_spectrum.xcat.old'

spec=read_xcat(file,/silent)

;ind=where(spec.error EQ 0.,count)
;IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec
ind=where(spec.sigmaII EQ 0.,count)
IF count NE 0 THEN spec(ind).sigmaII=(0.2*spec(ind).StokesI)^2

ind=where(spec.instru EQ 'FIRAS',count)

;IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec
IF count NE 0 THEN spec(ind).sigmaII=(0.2*spec(ind).StokesI)^2

;=== SET THE OBSERVATION STRUCTURE

st=dustem_set_data(st_sed_fit=spec);sed=spec)
;stop

;== SET THE FITTED PARAMETERS
dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims

;dustem_init_fixed_params,fpd,fiv

dustem_init_plugins, pd ;initializing the scopes of the plugins to be read later

;=== RUN fit

tol=1.e-16

;use_Nitermax=50 ;maximum number of iteration. This is the criterium which will stop the fit procedure , NOTA BENE: 3 iterations aren't enough to fit the entirety of the newly added stellar population parameters.
use_Nitermax=5 ;maximum number of iteration. This is the criterium which will stop the fit procedure , NOTA BENE: 3 iterations aren't enough to fit the entirety of the newly added stellar population parameters.

IF keyword_set(itermax) THEN use_Nitermax=itermax

t1=systime(0,/sec)

yr=[1.00e-4,1.00E2]
xr=[1.00E0,6.00e4]

;legend_xpos=0.6
;legend_ypos=0.8

tit='Spectral Energy Distribution'

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

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)

t2=systime(0,/sec)

;=== SAVE FIT RESULTS

file_out='/tmp/DUSTEM_fit_example.sav'

dustem_save_system_variables,file_out
message,'Saved '+file_out,/continue

;stop

;======================================
;====You can exit IDL here and re-enter
;======================================

file='/tmp/DUSTEM_fit_example.sav'
dustem_restore_system_variables,file

;=== Plot best fit

; tit='Spectral Energy Distribution (Saved)'
; ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
; xtit=textoidl('\lambda (\mum)')
; errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
; chi2=(*!dustem_fit).chi2
; rchi2=(*!dustem_fit).rchi2
; 
; window,2
; 
; ;=== RESTORE FIT RESULTS
; res=*(*!dustem_fit).current_param_values
; chi2=(*!dustem_fit).chi2
; rchi2=(*!dustem_fit).rchi2
; errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
; 
; ;=== Plot best fit
; loadct,13
; IF keyword_set(postcript) THEN BEGIN
;   set_plot,'PS'
; ;  ps_file=getenv('DUSTEM_SOFT_DIR')+'/Docs/Figures/'+'Last_dustem_fit.ps'
;   ps_file=!dustem_wrap_soft_dir+'Docs/Figures/'+'Last_dustem_fit.ps'
;   device,filename=ps_file,/color
; ENDIF
; 
; ;dustem_sed_plot,(*!dustem_fit_params),ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2
; dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog,legend_xpos=legend_xpos,legend_ypos=legend_ypos
; 
; 
; IF keyword_set(postcript) THEN BEGIN
;   device,/close
;   set_plot,'X'
;   message,'wrote '+ps_file,/info
; ENDIF
; IF keyword_set(png) THEN BEGIN
;   ;file_png=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_fit.png'
;   file_png=dir_png+'Last_dustem_fit_sed_readme.png'
;   write_png,file_png,tvrd(/true)
;   message,'Wrote '+file_png,/info
; ENDIF

tit='Spectral Energy Distribution (Final run)'
dustemwrap_plot,res,stp,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit

message,'dustem_mpfit_sed executed in '+strtrim(t2-t1,2)+' sec',/info

the_end:

END