PRO dustem_fit_sed_readme,postcript=postcript,mode=mode,help=help

;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,mode=mode,help=help
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    mode = 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,mode='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(mode) THEN BEGIN
  use_mode=strupcase(mode)
ENDIF ELSE BEGIN
  use_mode='COMPIEGNE_ETAL2010'    ;Default is last dustem model
ENDELSE

;dbp=1
;=== initialise dustem
;dustem_init,/wrap
;dustem_init
dustem_init,mode=use_mode
!dustem_verbose=1
!dustem_show_plot=1

;=== Read sample SED
;;dir=getenv('DUSTEM_SOFT_DIR')+'/src/'
;dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
;file=dir+'sample_SED.xcat'
;spec=read_xcat(file,/silent)
;=== Composite SED from Compiegne et al 2010, gathered by C. Bot
;dir=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/SEDs/'
dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
file=dir+'Gal_composite_spectrum.xcat'
;file=dir+'Gal_composite_spectrum_wAKARI.xcat'
;stop
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
;=== Draine Model SED for U=1. by C. Bot
;;dir=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/SEDs/'
;dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
;file=dir+'SED_DraineModel_U1.0.xcat'
;spec=read_xcat(file,/silent)
;=== Fake SED made from Dustem outputs
;spec=dustem_compute_sed(p_dim,st=st,cont=cont,_extra=extra)

;=== Set which parameters you want to fit
;CASE getenv('DUSTEM_WHICH') OF
CASE !dustem_which OF
  'WEB3p8':BEGIN
;=== This is to use Desert model
;=== must have done dustem_init,/DBP90
    CASE use_mode 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]
;        iv =   [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)
        ;=== Fixed parameters
;        fpd=['(*!dustem_params).G0'] ;Fixed parameter description
;        fiv=[1.5]
;        dustem_init_fixed_params,fpd,fiv
     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]
;        iv =   [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)
        ;=== Fixed parameters
;        fpd=['(*!dustem_params).G0'] ;Fixed parameter description
;        fiv=[1.0]
;        dustem_init_fixed_params,fpd,fiv
      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
        iv =   [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001]
;        iv =   [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)
        ;=== Fixed parameters
;        fpd=['(*!dustem_params).G0'] ;Fixed parameter description
;        fiv=[1.0]
;        dustem_init_fixed_params,fpd,fiv
      END
      'COMPIEGNE_ETAL2010':BEGIN
;     PUTTING G0 last DOES NOT WORK BECAUSE OF ORDER ISSUE
        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]
;        iv =   [7.8e-4, 7.7e-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)
        ;=== Fixed parameters
;        fpd=['(*!dustem_params).G0'] ;Fixed parameter description
;        fiv=[1.0]
;        dustem_init_fixed_params,fpd,fiv
      END
    ENDCASE
  END
  ELSE: BEGIN
     pd = ['(*!dustem_params).grains(0).propmass',$  ;PAH mass fraction
        '(*!dustem_params).grains(1).propmass',$  ;VSG mass fraction
        '(*!dustem_params).grains(2).propmass',$  ;BG mass fraction
        'dustem_create_continuum_2', $            ;Intensity of NIR continuum
        'dustem_create_isrf_1']                   ;ISRF scaling factor
    ;=== set initial value of parameters you want to fit
    iv =   [7e-5, 8e-4,6.6e-3,0.001,4.]
    ;in this example, we enforce parameters to be >0 and Xisrf >3
    ;Caution, initial values must be within limits ...
    Npar=n_elements(pd)
    ulimed=replicate(0,Npar)
    llimed=replicate(1,Npar)
    llims=replicate(0.,Npar)
  END
ENDCASE


;=== SET THE OBSERVATION STRUCTURE

st=dustem_set_data(sed=spec)

;VG : dustem_set_data turned into a function to be used indifferently for sed, ext, polext
;dustem_set_data,spec
;defsysv, '!dustem_data', {sed:ptr_new()}    ;Data to fit

;== SET THE FITTED PARAMETERS
dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
;stop
;dustem_sort_params
;hprint,*(*!dustem_fit).param_descs
;hprint,*(*!dustem_fit).param_init_values
;hprint,*(*!dustem_fit).param_func
;=== RUN fit
tol=1.e-14
;tol=1.e-2
;Nitermax=20               ;maximum number of iteration. This is the criterium which will stop the fit procedure
Nitermax=2               ;maximum number of iteration. This is the criterium which will stop the fit procedure

loadct,13
!y.range=[1e-4,10]              ;This is to ajust plot range from outside the routine
t1=systime(0,/sec)
res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax,gtol=gtol)
;res=dustem_mpfit_sed(tol=tol,Nitermax=Nitermax) ;,func_name=func_name,cf_min=cf_min)
t2=systime(0,/sec)

;=== With FXD Model
;Final parameters :    0.00012425071   0.00022732277    0.0021897348   0.00047272697       4.4916225
;Iter      6   CHI-SQUARE =       76.345245          DOF = 6
;=== With MC Model
;Final parameters :    0.00012292305   0.00020672300    0.0020245400   0.00076911024       4.9751991
;Iter      6   CHI-SQUARE =       76.047516          DOF = 6
;=== With LV Model
;Final parameters :    0.00013302325    0.0010706767    0.0039392709   0.00074149811       4.7549752
;Iter     10   CHI-SQUARE =       90.478020          DOF = 6
;=== With LV Model using GRAIN_DBP.DAT
;Iter      7   CHI-SQUARE =       117.27229          DOF = 6
;Final parameters :    0.00011163016   0.00015163894    0.0017526653   0.00096639410       5.8391213

;stop

;=== SAVE FIT RESULTS
;file_out=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
file_out=!dustem_res+'DUSTEM_fit_example.sav'
dustem_save_sed_fit,file_out

;=== Plot best fit
yr=[1e-4,10]
xr=[1,2000]
tit='DUSTEM Example'
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

loadct,13
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

;stop

;======================================
;====You can exit IDL here and re-enter
;======================================
window,1

;=== RESTORE FIT RESULTS
;=== initialise dustem
dustem_init,mode=mode
;file=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
file=!dustem_res+'DUSTEM_fit_example.sav'
dustem_restore_sed_fit,file

;=== recover best fit values
;ind=where(*!dustem_fit_chi2 EQ min(*!dustem_fit_chi2))
;res=(*!dustem_fit_params);
res=*(*!dustem_fit).current_param_values
chi2=(*!dustem_fit).chi2
rchi2=(*!dustem_fit).rchi2
;model_sed=*!dustem_fit_sed
;errors=*!dustem_fit_params_error;
errors=*(*!dustem_fit).current_param_errors

;=== Plot best fit
yr=[1e-4,10]
xr=[1,2000]
tit='DUSTEM Example (Saved)'
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
xtit=textoidl('\lambda (\mum)')

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
IF keyword_set(postcript) THEN BEGIN
  device,/close
  set_plot,'X'
  message,'wrote '+ps_file,/info
ENDIF

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

the_end:

END