PRO dustem_fit_sed_extsed_readme,postscript=postscript,model=model,help=help,itermax=itermax

;This Readme describes how to fit total intensity and polarization SEDs jointly
;Obviously, the dustem package must have been installed succesfully
;see DustemWrap documentation for install instructions.

;+
; NAME:
;    dustem_fit_sed_polsed_readme
; PURPOSE:
;    This is an example of how to fit total intensity and polarization SEDs jointly
;    It is meant to be an example to follow when writing your own
;    programs using the dustem IDL wrapper.
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_fit_sed_polsed_readme,postcript=postcript,mode=mode,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_sed_extsed_readme,/postscript

; MODIFICATION HISTORY:
;    Written by J.P. Bernard June 29th 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_extsed_readme'
  goto,the_end
ENDIF

;=== initialise dustem
;dustem_init,mode='COMPIEGNE_ETAL2010',/pol
use_model='COMPIEGNE_ETAL2010'
IF keyword_set(model) THEN use_model=model

dustem_init,mode=use_model

!dustem_verbose=1
!dustem_show_plot=1

IF keyword_set(postscript) THEN BEGIN
  ps_dir=!dustem_dat+'/Docs/Figures/'
  force_mkdir,ps_dir
ENDIF

;stop

;==== read emission SED

dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
file=dir+'Gal_composite_spectrum.xcat'
;file=dir+'Mathis1990_DISM.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

;==== read extinction SED

dir=!dustem_wrap_soft_dir+'/Data/EXTs/'
file=dir+'Mathis1990_DISM.xcat'
extsed=read_xcat(file,/silent)
ind=where(extsed.error EQ 0.,count)
IF count NE 0 THEN extsed[ind].error=0.2*spec(ind).spec


st=dustem_set_data(ext=extsed,sed=spec,rchi2_weight=rchi2_weight,f_HI=f_HI)

;=== Set which parameters you want to fit
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' $    ;amCBEx
              ]                   

p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04]    

;=== set initial value of parameters you want to fit
;iv=p_truth   ;exact solution
iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5]    ;shifted from solution

ulimed=[0   , 0   ,0   ,0 ]
llimed=[1   , 1   ,1   ,1 ]
llims =[0.  ,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-10
xtol=1.e-10

Nitermax=2              ;maximum number of iteration. This is the criterium which will stop the fit procedure

IF keyword_set(itermax) THEN Nitermax=itermax

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)

yr=[1e-4,10]
xr=[1,1.e5]
legend_xpos=0.6
legend_ypos=0.8

tit='DUSTEMWrapper Intensity SED (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,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit=xtit,ytit=ytit,legend_xpos=legend_xpos,legend_ypos=legend_ypos,title=tit)

t2=systime(0,/sec)

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

;=== SAVE FIT RESULTS
file_out='/tmp/DUSTEM_sedextsed_fit_example.sav'
dustem_save_system_variables,file_out
message,'Saved '+file_out,/continue

;======================================
;====You can exit IDL here and re-enter
;======================================
file='/tmp/DUSTEM_sedextsed_fit_example.sav'
dustem_restore_system_variables,file

;=== Plot best fit intensity SED
win=0L

window,win & win=win+1

cleanplot

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

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)

tit='DUSTEMWrapper polarization dust Example (Saved)'

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

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,legend_xpos=legend_xpos,legend_ypos=legend_ypos

loadct,13
IF keyword_set(postscript) THEN BEGIN
  previous_device=!d.name
  set_plot,'PS'
  !p.charsize=0.8
  ps_file=ps_dir+'Last_dustem_sed_fit.ps'
  device,filename=ps_file,/color

  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,legend_xpos=legend_xpos,legend_ypos=legend_ypos

  device,/close
  set_plot,previous_device
  message,'wrote '+ps_file,/info
ENDIF

;stop
;=== Plot best fit Extinction SED

loadct,13

window,win & win=win+1
cleanplot
xr=[0.3,40]
yr=[1.e-26,1.e-21]
dustem_extsed_plot,*(*!dustem_fit).current_param_values,pol=pol,xr=xr,yr=yr,/xsty,/ysty

IF keyword_set(postscript) THEN BEGIN
  previous_device=!d.name

  set_plot,'PS'

  !p.charsize=0.8
  ps_file=ps_dir+'Last_dustem_ext_ir_fit.ps'

  device,filename=ps_file,/color

  dustem_extsed_plot,*(*!dustem_fit).current_param_values,pol=pol,xr=xr,yr=yr,/xsty,/ysty

  device,/close
  set_plot,previous_device
  message,'wrote '+ps_file,/info

ENDIF

;stop

loadct,13
window,win & win=win+1
cleanplot
xr=[0.2,10]
yr=[0,3e-21]
dustem_extsed_plot,*(*!dustem_fit).current_param_values,pol=pol,xr=xr,yr=yr,/xsty,/ysty,/UV

IF keyword_set(postscript) THEN BEGIN
  previous_device=!d.name
  set_plot,'PS'
  !p.charsize=0.8
  ps_file=ps_dir+'Last_dustem_ext_uv_fit.ps'
  device,filename=ps_file,/color
  
  dustem_extsed_plot,*(*!dustem_fit).current_param_values,pol=pol,xr=xr,yr=yr,/xsty,/ysty,/UV
  

  device,/close

  set_plot,previous_device

  message,'wrote '+ps_file,/info
ENDIF

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

stop

the_end:

END