PRO dustem_fit_sed_polsed_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_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_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_polsed_readme'
  goto,the_end
ENDIF

;=== initialise dustem
;dustem_init,mode='COMPIEGNE_ETAL2010',/pol
use_mode='G17_MODELA'

dustem_init,mode=use_mode,/pol

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

;stop

;=== 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']
filters=['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)
polsed=sed ; comment if polext is to be used
;polext=sed ;comment if polsed is to be used

;UNCOMMENT THIS BLOCK AFTER THE TESTING

; ;stokes fake parameters - comment if stokes extinction parameters are to be used
; qsed=sed
; used=sed

; ;stokes extinction fake parameters - comment if stokes parameters are to be used
; qext=sed
; uext=sed

st=dustem_set_data(ext=ext,polsed=polsed,polext=polext,sed=sed,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI,qsed=qsed,used=used,qext=qext,uext=uext)


;=== 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_create_synchrotron_2',$ 
;              'dustem_create_stokes_1',$ ;Polarization angle psi_ref_dust
;              'dustem_create_stokes_2',$;Polarization angle psi_ref_synch
;              'dustem_create_stokes_3'$;Polarization fraction of synchrotron spectrum
;              ;'dustem_create_stext_1'$;Polarization angle in extinction PSI_ext
;               ]                   

;TESTING THE ISEF PLUGIN
; 
; pd = [ $i
;              ;'(*!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 
;              'dustem_create_stellar_population_O63',$ ; distance to the O6 stellar population
;              ;'dustem_create_stellar_population_O65',$ ; number of stars of the O6 stellar population
;              'dustem_create_stellar_population_B43'$ ; distance to the B4 stellar population
;              ;'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population 
;               ]                   

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
      'dustem_create_stellar_population_O63', $ ; distance to the O6 stellar population
      'dustem_create_stellar_population_O65', $ ; number of stars of the O6 stellar population
      'dustem_create_stellar_population_B43', $ ; distance to the B4 stellar population
      'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population
        ] 
 ;p_truth =   [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001,13.,3.,10.,1.]
 ;p_truth =   [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001,13.,3.,10.,1.]   
  p_truth =   [1.8e-4,2.33e-3,8.27e-3,0.001,13.,3.,10.,1.]


;UNCOMMENT THIS WHOLE BLOCK AFTER THE TESTING

; p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,70,20,0.2]    
; ;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,25]    


;UNCOMMENT THIS BLOCK AFTER THE TESTING
;p_truth=[1.,7.8000E-04,7.8000E-04,2E-03,5e-3,70,20,0.2]

;TESTING THE ISRF PLUGIN

;p_truth=[1.,7.8000E-04,7.8000E-04,2E-03,2e-3,13.,3.,9.,2.]
;p_truth=[7.8000E-04,7.8000E-04,2E-03,2e-3,13.,3.];,9.,2.]
;p_truth=[7.8000E-04,7.8000E-04,2E-03,10.,8.]
;p_truth=[10.,8.]


;=== set initial value of parameters you want to fit
iv=p_truth   ;exact solution



;UNCOMMENT THIS BLOCK AFTER THE TESTING

; iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5,5e-2,5,5,0.2];shifted from solution
; ;iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5,5];shifted from solution

;TESTING THE ISRF PLUGIN

;iv=p_truth+[0,3.e-5,1.e-5,-2.e-5,9e-4,-3.,-2,2.,1.]
;iv=p_truth+[3.e-5,1.e-5,-2.e-5,9e-4,-3.,-2];,2.,1.]
;iv=p_truth+[3.e-5,1.e-5,-2.e-5,-2,-2]
;iv=p_truth+[-2,-2]

;iv=p_truth+[0.,5.4e-5, 5.4e-5,1.8e-5,2.33e-4,8.27e-4,1.0e-4,-2.,-2,3.,3.]
;iv=p_truth+[5.4e-5, 5.4e-5,1.8e-5,2.33e-4,8.27e-4,1.0e-4,-2.,-2,3.,3.]

iv=p_truth+[1.8e-5,2.33e-4,8.27e-4,1.0e-4,-2.,-2,3.,3.]


;UNCOMMENT THIS BLOCK AFTER THE TESTING

; 
; ulimed=[0   ,0   ,0   ,0   ,0   ,0   , 0   ,0] ;Polarization angles exist in the range [0-180] degress
; 
; ;ulimed=[0   ,0   ,0   ,0   ,0] ;Polarization angles exist in the range [0-180] degress
; 
; 
; llimed=[1   , 1   ,1   ,1   ,1   ,1   ,1   ,1]
; ;llimed=[1   , 1   ,1   ,1   ,1]
; 
; llims =[0.  ,0.   ,0.  ,0.   , 0.   ,0.   ,0.   ,0.]
; ;llims =[0.  ,0.   ,0.  ,0.   ,0.]


;TESTING THE ISRF PLUGIN

Npar=n_elements(pd)  
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
;llims =[1e-5  ,0.   ,0.    ,0.   ,0.   ,1.  ,0.   ,1.   ,0.   ]

;llims =[0.   ,0.    ,0.   ,0.   ,1.   ,0.   ];,1.   ,0.   ]
;llims =[0.   ,0.    ,0.   ,1.   ,1.   ]
;llims =[1.   ,1.   ]

;llims = [1.0e-5   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,1.   ,1.   ,1.   ,1.   ]
;llims = [0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,1.   ,1.   ,1.   ,1.   ]

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

;== 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_init_fixed_params,fpd,fiv
    
;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
;(*!dustem_params).grains[4].TYPE_KEYWORDS='pol-plaw-ed'


dustem_init_plugins, pd

help,*!dustem_fit

help,!run_pol

sed.spec=dustem_compute_sed(p_truth,sst,_extra=extra)   ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
sed.error=(sed.spec*0.1);>0.1

;UNCOMMENT THIS BLOCK AFTER THE TESTING
;polsed.spec=dustem_compute_polsed(p_truth,sstp,_extra=extra,dustem_qsed,dustem_used)   ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code

polsed.spec=dustem_compute_polsed(p_truth,sstp,_extra=extra)

polsed.error=(polsed.spec*0.1)


; ;creating fake polext
;   
; polext.spec=dustem_compute_polext(p_truth,ssttp,dustem_qext,dustem_uext)
; polext.error=(polext.spec*0.1)
  


;UNCOMMENT THIS BLOCK AFTER THE TESTING

; ;SET qsed and used data - comment if extinction polarization data is to be used
; qsed.spec=dustem_qsed
; qsed.error=qsed.spec*0.1
; 
; 
; used.spec=dustem_used
; used.error=used.spec*0.1


; ;SET qext and uext data - comment if polarizarion data is to be used
; qext.spec=dustem_qext
; qext.error=qext.spec*0.1
; 
; 
; uext.spec=dustem_uext
; uext.error=uext.spec*0.1
; ; 

;=== SET THE OBSERVATION STRUCTURE

ind=where(sed.wave GE 0.1)

;UNCOMMENT THIS BLOCK AFTER THE TESTING
; st=dustem_set_data(polsed=polsed[ind],sed=sed[ind],qsed=qsed[ind],used=used[ind]) ; comment if second line is to be used
; ;st=dustem_set_data(sed=sed[ind],polext=polext[ind],qext=qext[ind],uext=uext[ind])


st=dustem_set_data(polsed=polsed[ind],sed=sed[ind])






;=== RUN fit
tol=1.e-16
xtol=1.e-16
Nitermax=30; 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-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.

;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