PRO dustem_fit_spin_test_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_spin_test_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_spin_test_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_spin_test_readme'
  goto,the_end
ENDIF

;=== initialise dustem
dustem_init,mode='COMPIEGNE_ETAL2010'
!dustem_verbose=1
!dustem_show_plot=1

;=== make a fake one out for given filters

filters=['DIRBE1','DIRBE2','DIRBE3','DIRBE4','IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']

Nfilt=n_elements(filters)
sed=dustem_initialize_sed(Nfilt)
sed.filter=filters
sed.wave=dustem_filter2wav(filters)
sed.instru=dustem_filter2instru(filters)
st=dustem_set_data(sed=sed)

;1.000000
;PAH0_MC10         25 logn-chrg-spin-zm           7.8000E-04  2.2400E+00  3.5000E-08  1.2000E-07  6.4000E-08  1.0000E-01
;PAH1_MC10         25 logn-chrg-spin-zm           7.8000E-04  2.2400E+00  3.5000E-08  1.2000E-07  6.4000E-08  1.0000E-01
;amCBEx            15 logn                        1.6500E-04  1.8100E+00  6.0000E-08  2.0000E-06  2.0000E-07  3.5000E-01
;amCBEx            25 plaw-ed                     1.4500E-03  1.8100E+00  4.0000E-07  2.0000E-04 -2.8000E+00  1.5000E-05  1.5000E-05  2.0000E+00
;aSilx             25 plaw-ed                     6.7000E-03  3.0000E+00  4.0000E-07  2.0000E-04 -3.4000E+00  2.0000E-05  2.0000E-05  2.0000E+00

goto,freefree_tag

;=== Set which parameters you want to fit
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', $    ;aSilx
             '(*!dustem_params).gas.G0', $                  ;Gas G0
             '(*!dustem_params).gas.Tgas', $                ;Gas temperature
             '(*!dustem_params).gas.nH', $                  ;Gas temperature

;             'dustem_create_freefree_1', $                 ;Gas temperature for free-free emission
;             'dustem_create_freefree_2', $                 ;H-alpha for free-free emission (R)

             'dustem_create_continuum_2']                   ;Intensity of NIR continuum

;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]    
;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,0.01,0.1]    
p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1]    

;=== set initial value of parameters you want to fit

;iv=p_truth   ;exact solution
;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,0.01,0.001]    ;shifted
iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001]    ;shifted

ulimed=[0   , 0   ,0,0,0,0,0,0,0]
llimed=[1   , 1   ,1,1,1,1,1,1,1]
llims =[0.  ,0.   ,0.,0.,0.,0.,0.,0.,0.]

goto,out_of_there

freefree_tag:
;=== Set which parameters you want to fit
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', $    ;aSilx
             '(*!dustem_params).gas.G0', $                  ;Gas G0
             '(*!dustem_params).gas.Tgas', $                ;Gas temperature
             '(*!dustem_params).gas.nH', $                  ;Gas temperature
             'dustem_create_freefree_1', $                 ;Gas temperature for free-free emission
             'dustem_create_freefree_2', $                 ;H-alpha for free-free emission (R)
             'dustem_create_continuum_2']                   ;Intensity of NIR continuum

;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]    
p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,1.e-11,0.1]    
;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1]    

;=== set initial value of parameters you want to fit
;iv=p_truth   ;exact solution
iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,1.e-11,0.001]    ;shifted
;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001]    ;shifted

ulimed=[0   , 0   ,0,0,0,0,0,0,0,0,0]
llimed=[1   , 1   ,1,1,1,1,1,1,1,1,1]
llims =[0.  ,0.   ,0.,0.,0.,0.,0.,0.,0.,0,0]

synchrotron_tag:
;=== Set which parameters you want to fit
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', $    ;aSilx
             '(*!dustem_params).gas.G0', $                  ;Gas G0
             '(*!dustem_params).gas.Tgas', $                ;Gas temperature
             '(*!dustem_params).gas.nH', $                  ;Gas temperature
             'dustem_create_synchrotron_1', $                 ;CR spectral index for synchrotron emission
             'dustem_create_synchrotron_2', $                 ;Amplitude for synchrotron emission (R)
             'dustem_create_continuum_2']                   ;Intensity of NIR continuum

;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]    
p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,3.,0.01,0.1]    
;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1]    

;=== set initial value of parameters you want to fit
;iv=p_truth   ;exact solution
iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,3.,0.01,0.1]    ;shifted
;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001]    ;shifted

ulimed=[0   , 0   ,0,0,0,0,0,0,0,0,0]
llimed=[1   , 1   ,1,1,1,1,1,1,1,1,1]
llims =[0.  ,0.   ,0.,0.,0.,0.,0.,0.,0.,0,0]

out_of_there:

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

(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'

help,*!dustem_fit

;stop

sed.spec=dustem_compute_sed(p_truth,sst,_extra=extra)
sed.error=(sed.spec*0.1);>0.1
;sed.error[*]=100.

;=== SET THE OBSERVATION STRUCTURE

;ind=where(sed.wave GE 5.)
ind=where(sed.wave GE 0.1)

st=dustem_set_data(sed=sed[ind])


;=== RUN fit
;tol=1.e-20
tol=1.e-10
xtol=1.e-10
Nitermax=5              ;maximum number of iteration. This is the criterium which will stop the fit procedure

xrange=[1.,2.e4]
yrange=[1e-4,1.e3]

loadct,13
;!y.range=[1e-8,100]              ;This is to ajust plot range from outside the routine
t1=systime(0,/sec)

res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []')

t2=systime(0,/sec)

print,(res-p_truth)/p_truth*100.

;stop

;=== SAVE FIT RESULTS

file_out='/tmp/DUSTEM_spinning_fit_example.sav'
dustem_save_system_variables,file_out
message,'Saved '+file_out,/continue

stop

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=yrange,xr=xrange,/ysty,/xsty, $
                res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog

;stop

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

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

window,1

xrange=[1.,2.e4]
yrange=[1e-4,1.e3]

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

;errors=*(*!dustem_fit).current_param_errors
errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)

;=== Plot best fit

tit='DUSTEM spinning dust 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_spinning_fit.ps'

  device,filename=ps_file,/color
ENDIF

dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
;dustem_sed_plot,*(*!dustem_fit).current_param_values, $
;                ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog

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