PRO dustem_fit_modified_bb_readme,postcript=postcript,mode=mode,help=help

;This Readme describes how to fit SEDs with a modified Black-Body
;It runs on the SED stored into the file sample_SED.xcat
;Remember that the goal here is just to illustrate the method.
;Note that its is not 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 DustemWrap documentation for install instructions.

;+
; NAME:
;    dustem_fit_modified_bb_readme
; PURPOSE:
;    This is an example of how to fit SEDs with a modified Black-Body
;    using the dustem wrapper.
;    It is meant to be an example to follow when writing your own
;    programs using the dustem IDL wrapper.
;    The SED used here is in sample_SED.xcat
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_fit_modified_bb_readme,postcript=postcript,help=help
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    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 idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    dustem_fit_modified_bb_readme
; MODIFICATION HISTORY:
;    Written by J.P. Bernard June 9th 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_modified_bb_readme'
  goto,the_end
ENDIF

;=== initialise dustem
dustem_init
!dustem_verbose=1
!dustem_show_plot=1

;=== Read sample SED
;=== 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'
;sed=read_xcat(file,/silent)
;=== Or make a fake one out for given filters
;filters=['IRAS4','HFI1','HFI2','HFI3']
filters=['IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3']
Nfilt=n_elements(filters)
sed=dustem_initialize_sed(Nfilt)
sed.filter=filters
sed.wave=dustem_filter2wav(filters)
sed.instru=dustem_filter2instru(filters)

;VG : dustem_set_data turned into a function to be used indifferently for sed, ext, polext
;dustem_set_data,sed  ;This is to set filters in !dustem_data.sed. Needed by dustem_compute_gb_sed
;!dustem_data.sed=dustem_set_data(sed) ;This is to set filters in !dustem_data.sed. Needed by dustem_compute_gb_sed
st=dustem_set_data(sed=sed)

;p_truth=[0.5D0,22.D0,1.8D0]
p_truth=[0.5,19.,1.8]
sed.spec=dustem_compute_gb_sed(p_truth,_extra=extra)
sed.error=sed.spec*0.05


;=== SET THE OBSERVATION STRUCTURE

ind=where(sed.wave GE 100.)
;VG : dustem_set_data turned into a function to be used indifferently for sed, ext, polext
;dustem_set_data,sed(ind)
st=dustem_set_data(sed=sed(ind))
;!dustem_data.sed=dustem_set_data(sed(ind))

;=== Set which parameters you want to fit
pd = ['Amplitude','T','beta']
;=== set initial value of parameters you want to fit
;iv =   [1.,20.,2.]
;iv =   [0.5,22.D0,1.5D0]
iv=[0.5,22.,1.8]
ulimed=[0   , 0   ,0]
llimed=[1   , 1   ,1]
llims =[0.  ,0.   ,0.]

;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims


;=== RUN fit
;tol=1.e-20
tol=1.e-20
xtol=1.e-10
Nitermax=100               ;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,xtol=xtol,Nitermax=Nitermax,function_name='dustem_greybody_mpfit')
t2=systime(0,/sec)

;stop

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

;=== Plot best fit
yr=[1e-3,100]
xr=[50,4000]
tit='DUSTEM Modified BB 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)
errors=(*(*!dustem_fit).current_param_errors)
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,function_name='dustem_greybody_mpfit'

;stop

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

;=== RESTORE FIT RESULTS
;=== initialise dustem
dustem_init
file=!dustem_res+'DUSTEM_bb_fit_example.sav'
dustem_restore_sed_fit,file

;=== recover best fit values
res=*(*!dustem_fit).current_param_values
chi2=(*!dustem_fit).chi2
rchi2=(*!dustem_fit).rchi2
errors=*(*!dustem_fit).current_param_errors

;=== Plot best fit
yr=[1e-3,100]
xr=[50,4000]
tit='DUSTEM Modified BB 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=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_bb_fit.ps'
  device,filename=ps_file,/color
ENDIF
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,function_name='dustem_greybody_mpfit'
IF keyword_set(postcript) THEN BEGIN
  device,/close
  set_plot,'X'
  message,'wrote '+ps_file,/info
ENDIF

print,'dustem_mpfit_sed executed in ',t2-t1,' sec'

the_end:

END