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

;COMMENT: for some reson polarization=1 doesn't work for ALL of the models except for the two last A and B 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_create_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_create_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_create_continuum_2']                ;Intensity of NIR continuum
              ]  
        
;         pd = [ $
;             '(*!dustem_params).gas.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_create_continuum_2', $
;             'dustem_create_stellar_population_O63',$ ; distance to the O3 stellar population
;             'dustem_create_stellar_population_O65',$ ; number of stars of the O3 stellar population
;             'dustem_create_stellar_population_B43',$ ; distance to the B4 stellar population
;             'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population
;             ]
        
        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_create_continuum_2', $      ;intensity of NIR continuum
;               
;               ;'!dustem_isrf_star_add[0].distance', $        ;distance of first stellar contribution to ISRF
;               '!dustem_isrf_star_add[0].amplitude', $        ;amplitude of first stellar contribution to 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' $
;               ]
;         ;initial parameter values for parameters to be fitted   
;         iv =   [ $  
;           1.0, $
;           0.002,   $ ;intensity of NIR continuum  
;           ;0.003,   $ ;distance of first stellar contribution to ISRF  
;           1.e-3,   $ ;amplitude of first stellar contribution to ISRF  
;           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  
;           ]
          
;NEW LINES

;========COMMENTS=======
;amplitudes haven't been used so they're by default set to 1.
;using only the below stellar populations produces amplitudes that exceed 1. I do not understand the physical meaning of this.
         
;         pd = [ $
;               '(*!dustem_params).gas.G0', $      ;G0
;               'dustem_create_continuum_2', $      ;intensity of NIR continuum
;               'dustem_create_stellar_population_O63',$ ; distance to the O3 stellar population
;               ;'dustem_create_stellar_population_O64',$ ; amplitude of the O3 stellar population
;               'dustem_create_stellar_population_O65',$ ; number of stars of the O3 stellar population
;               'dustem_create_stellar_population_B43',$ ; distance to the B4 stellar population
;               ;'dustem_create_stellar_population_B44',$ ; amplitude of the B4 stellar population
;               'dustem_create_stellar_population_B45',$ ; number of stars of the B4 stellar population
;               '(*!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' $
;               ]  
              
        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_params).grains(3).mdust_o_mh', $    ;amCBEx
              '(*!dustem_params).grains(4).mdust_o_mh', $
              ;'dustem_create_continuum_2', $      ;intensity of NIR continuum
              'dustem_create_stellar_population_O63',$ ; distance to the O6 stellar population
              
              ;'dustem_create_stellar_population_O64',$ ; amplitude of the O6 stellar population
              
              ;'dustem_create_stellar_population_O65',$ ; number of stars of the O6 stellar population
              'dustem_create_stellar_population_B43'$ ; distance to the B4 stellar population
              
              ;;'dustem_create_stellar_population_B44',$ ; amplitude of the B4 stellar population
              
              ;'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population
              ]        
 
          
        ;initial parameter values for parameters to be fitted   
;         iv =   [ $  
;           1.,     $ ;G0
;           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 big grains?
;           0.002,   $ ;intensity of NIR continuum    
;           10.,$
;           ;;1.,$
;           1.,$
;           10.,$
;           ;;1.,$
;           1.$  
;           ]          
;           
        
        iv =   [ $  
          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 big grains?
          ;0.002,   $ ;intensity of NIR continuum    
          10.,$
          ;;1.,$
          ;1.,$
          10.$
          ;;1.,$
          ;1.$  
          ]          
        
        
            
     
        Npar=n_elements(pd)
        
        ;ulimed=[0   ,0   ,0   ,0   ,0   ,0   , 0   ,0   ,0   ,0   ,0   ,0   ,0] 
        ;ulimed=[0   ,0   ,0   ,0   ,0   ,0   , 0   ,0   ,0   ,0] 

        
        ;llimed=[1   , 1   ,1   ,1   ,1   ,1   ,1   ,1   ,1   ,1   ,1   ,1   ,1]
        ;llimed=[1   , 1   ,1   ,1   ,1   ,1   ,1   ,1   ,1   ,1]
        
        ;llims =[0.  ,0.   ,0.  ,0.   , 0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.    ,0.   ,0.   ,0.]
        ;llims =[1.  ,0.   ,10.  ,0.   , 0.   ,0.   ,0.   ,0.   ,0.    ,0.   ,0.]
        
        ;lower limiting the stellar distances to 1pc avoid singularities
        ;[1e-9  ,
        ;llims =[1e-5  ,0.   ,0.    ,0.   ,0.   ,0.   ,0.   ,1.  ,0.     , 1.   ,0.   ]  
        
        ;llims =[0.   ,0.    ,0.   ,0.   ,0.   ,0.   ,1.  ,0.     , 1.   ,0.   ]
        llims =[0.    ,0.   ,0.   ,0.   ,0.   ,1.  , 1.   ]
        
        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_create_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
        ;dustem_init,model=use_model,/pol

         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_create_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_create_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

ENDCASE

;dustem_init,model=use_model,pol=polarization

!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'

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.instru EQ 'FIRAS',count)
IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec



;=== SET THE OBSERVATION STRUCTURE
st=dustem_set_data(sed=spec)

;== 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.
IF keyword_set(itermax) THEN use_Nitermax=itermax

t1=systime(0,/sec)
yr=[1e-4,100]
xr=[1,6e4]
;legend_xpos=0.6
;legend_ypos=0.8
tit='Spectral Energy Distribution (Running)'
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

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

the_end:

END