PRO dustem_fit_sed_qsed_used_readme,postcript=postcript,mode=mode,help=help

;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_qsed_used_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_polsed_readme
; 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_qsed_used_readme'
  goto,the_end
ENDIF


dustem_define_la_common

;=== initialise dustem
;dustem_init,mode='COMPIEGNE_ETAL2010',/pol
use_mode='THEMIS';'G17_MODELA' ; try with model D 

dustem_init,mode=use_mode,/pol

!dustem_verbose=1
!dustem_show_plot=1
help,!run_pol

;stop

!dustem_verbose=1
!dustem_show_plot=1


;filters=['DIRBE1','DIRBE2','DIRBE3','DIRBE4','IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','PILOT1','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)



;INITIALIZING IQU SEDS and associated errors
sed.StokesI=sed.wave*0
sed.StokesQ=sed.wave*0
sed.StokesU=sed.wave*0
sed.sigmaII=sed.wave*0
sed.sigmaQQ=sed.wave*0
sed.sigmaUU=sed.wave*0



;INITIALIZING THE !dustem_data tags that will be internally used 
st=dustem_set_data(sed=sed,rchi2_weight=rchi2_weight,f_HI=f_HI)



;=== Set which parameters you want to fit
pd = [ $
             '(*!dustem_params).G0', $      ;G0
             '(*!dustem_params).grains(0).mdust_o_mh',$  
             '(*!dustem_params).grains(1).mdust_o_mh',$     
             ;'(*!dustem_params).grains(1).mdust_o_mh',$    
             '(*!dustem_params).grains(2).mdust_o_mh',$    
             ;'dustem_plugin_synchrotron_2', $               ;Synchrotron amplitude
             ;'dustem_plugin_synchrotron_4', $               ;Synchrotron polarization angle
             'dustem_plugin_modify_dust_polarization_2',$ ;Polarization angle applied to the Fortran dust model via a core plugin
             'dustem_plugin_modify_spinning_polarization_1', $
             'dustem_plugin_modify_spinning_polarization_2' $
              ]                   

;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,68.,0.03,30.]    
;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,68.,0.03,30.]
;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04];,0.01,30.]



p_truth=[1.,0.17E-03,0.63E-03,0.51E-03,68.,0.03,30.]



;p_truth=[68.]


; fpd = [ $
;              '(*!dustem_params).G0', $      ;G0
;              '(*!dustem_params).grains(0).mdust_o_mh',$  
;              '(*!dustem_params).grains(1).mdust_o_mh',$     
;              '(*!dustem_params).grains(2).mdust_o_mh'$    
; ;              '(*!dustem_params).grains(3).mdust_o_mh',$    
; ;              'dustem_plugin_modify_dust_polarization_2', $ ;Polarization angle applied to the Fortran dust model via a core plugin
; ;              'dustem_plugin_modify_spinning_polarization_1' $
;                ]                   
; ; 
; fiv=[1.,7.8000E-04,7.8000E-04,1.6500E-04];,68.,0.01]


fixed=[1,1,1,1,0];0,1,1]
 
;dustem_init_parinfo,pd,p_truth;,fixed=fixed

;dustem_init_fixed_params,fpd,fiv

iv=p_truth+[0.,5.000E-04,5.000E-04,5.000E-03,-10.,1.00e-2,10.]      ;shifted from solution ,10.


;iv=p_truth+[0.,5.000E-04,5.000E-04,5.000E-04];,5.e-4,1.]      ;shifted from solution ,10.


Npar=n_elements(pd)  

ulimed=replicate(0,Npar);[0]

;ulimed=[0,0,0,0,0,1,1,1]
;ulims=[0.,0.,0.,0.,0.,100.,1.,100.]

llimed=replicate(1,Npar);[1]
llims=replicate(0,Npar);[10]

dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims;,fixed=fixed

;dustem_init_fixed_params,fpd,fiv

dustem_init_plugins,pd;, fpd=fpd



sed.StokesI=dustem_compute_sed(p_truth,ste)   ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
sed.sigmaII=double((sed.StokesI*0.3)^2)

toto=dustem_compute_stokes(p_truth,stee,dustem_qsed,dustem_used) ;this procedure also allows for the extraction of the spectra

sed.StokesQ=dustem_qsed
sed.sigmaQQ=(double((sed.StokesQ*0.03)^2)) > 1d-40
;tes=where(finite(sed.sigmaQQ) eq 0,coun)
;if coun ne 0 then sed.sigmaQQ(tes)=0.

sed.StokesU=dustem_used
sed.sigmaUU=(double((sed.StokesU*0.03)^2)) > 1d-40
;tes=where(finite(sed.sigmaUU) eq 0,coun)
;if coun ne 0 then sed.sigmaUU(tes)=0.


;======================================================================================================================



; llimed=[1];replicate(1,Npar)
; llims=[1e-12];replicate(0,Npar)

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


;dustem_init_plugins,pd,fpd=fpd



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


; 
; sed.StokesI=dustem_compute_sed(p_truth,ssti)   ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
; sed.sigmaII=(sed.StokesI*0.3)^2
; 
; 
; toto=dustem_compute_stokes(p_truth,sstqu,dustem_qsed,dustem_used) ;this procedure also allows for the extraction of the spectra
; 
; sed.StokesQ=dustem_qsed
; sed.sigmaQQ=(double(sed.StokesQ)*0.4^2) > 1d-50
; 
; 
; sed.StokesU=dustem_used
; sed.sigmaUU=(double(sed.StokesU)*0.4^2) > 1d-50


;=== SET THE OBSERVATION STRUCTURE

ind=where(sed.wave GE 0.1)

st_bidon=dustem_set_data(sed=sed[ind],rchi2_weight=rchi2_weight,f_HI=f_HI) ; comment if second line is to be used

;=== RUN fit
tol=1.e-40
xtol=1.e-40
Nitermax=23; just an initialization


;for i=0L,n_tags(!dustem_data)-1 do begin
;if ptr_valid(!dustem_data.(i)) then Nitermax+=1
;endfor
;Nitermax=4           ;maximum number of iteration. This is the criterion which will stop the fit procedure

xrange=[1.,2.e4]
yrange=[1e-7,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.

;help,(*!dustem_fit),/STRUCTURE

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

;stop

;======================================
;====You can exit IDL here and re-enter   -   THIS PART APPEARS TO BE FAULTY. NEEDS INVESTIGATION 
;======================================
file='/tmp/DUSTEM_polsed_fit_example.sav'
dustem_restore_system_variables,file

;=== Plot best fit
win=1
window,win & win=win+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 polarization 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/'+'Last_dustem_spinning_fit.ps';'/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,/pol

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