PRO dustem_fit_sed&polsed&polext_NY


;Thanks for this! You're the first user that can give me actual feedback.
;I realized that there is plenty of work that still needs to be done. 
;The wrapper is a war zone at the moment but we're getting there.
 
 
dustem_define_la_common ;Enables the wrapper to use a value called la_undef() 
;which is how the wrapper recognizes undefined values.
;This value equals -32768 and should be present in your xcat files for all the values that you wish to omit from the fit. 
;Setting unused/unwanted values/errors to zero is highly unrecommended (unless your case makes it reasonable to assume a null flux value. ie StokesQ = 0. for some band)
;I'm still working w/ JP on this so you shouldn't really worry about this. 
 
 
;☆☆☆DUSTEM INITIALIZATION☆☆☆
model='THEMIS' 
dustem_init,mode=model,/pol

!dustem_verbose=1
!dustem_show_plot=1

;☆☆☆READING THE DATA☆☆☆
dir_seds = '/Users/ilyeschoubani/Downloads/test_wrap_NY/'

;###README### 
;If the format of the SED is the new one (the current one), one just needs to read the xcat file. 
;Otherwise, a procedure called dustem_old2new_sed_format.pro does the trick. That is what I used for your SED xcat file.
;This new format has been adjusted so that the data pertaining to the polarized emission is included in the SED (primary/main) structure. 
;In other words, the polarization SED should also be there, which isn't our case since the polarization SED is in a seperate xcat file.
;Because of that, we will be reading the SED (which has the new format now) and fill the polarization-related tags (ie: P,p,their associated variances etc...) by reading their xcat files first (in our case we just have polsed). 

;☆☆☆Reading the SED xcat that has the new format☆☆☆
str_sed = dir_seds+'SED_HD21.xcat' ;It has the same name because the old one has 'old' attached to its string now (after running the aforementioned procedure).
sed = read_xcat(str_sed,/silent)

;☆☆☆Reading the POLSED xcat and including its data in the (primary/main) sed structure☆☆☆  

;☆☆☆reading polsed xcat file☆☆☆
str_polsed = dir_seds+'POLSED_HD21.xcat'
polsed = read_xcat(str_polsed,/silent)

;☆☆☆now adding its content to the SED main structure☆☆☆
match, sed.wave, polsed.wave, subsed, subpolsed
sed(subsed).LARGEP = polsed(subpolsed).spec
sed(subsed).sigma_LARGEP = polsed(subpolsed).error^2 ;The squared error value here is because the wrapper uses variances now. When preparing to fit,the wrapper computes the square root of this value so don't be alarmed as the error remains the same.

;###README: 
;The wrapper follows the same formalism when it comes to the handling of the extinction data.
;Meaning there is also a (primary/main) structure for extinction that also contains data pertaining to polarization (cross-section in the case of polext).
;Since we do not have extinction data (only Serkowski), we should read the Serkowski xcat file, and use its content to create a fake extinction (primary/main) procedure.
;The creation of this fake extinction procedure is necessary for dustem to be aware of the polext data.
;However I never implemented this before. By 'this' I mean the way the wrapper fits polext without the presence of ext.
;This means that your work has enabled us to think about a very important issue. 


;☆☆☆The first thing we do is read the Serkowski xcat file that has the new format☆☆☆
;###README:
;I used a procedure called dustem_old2new_sed_format.pro that does the job.
;This means that I have created the (primary/main) extinction structure with the Serkowski xcat file.
;The lines below are witness to the necessary changes that need to be made so that said structure includes the polext data.  

;☆☆☆Reading the POLEXT xcat file☆☆☆
str_polext = dir_seds+'Serkowski_HD21.xcat'
ext = read_xcat(str_polext,/silent)


;☆☆☆Modifying the extinction structure and placing the polext data inside it☆☆☆
ext.instru = 'EXTINCTION'
ext.EXT_P = ext.EXT_I
ext.SIGEXTP = ext.SIGEXTII
ext.EXT_I = la_undef(4)
ext.SIGEXTII = la_undef(4)

;=== Set which parameters you want to fit
pd = [ $
             '(*!dustem_params).gas.G0', $      ;G0
             '(*!dustem_params).grains(0).mdust_o_mh', $
             '(*!dustem_params).grains(1).mdust_o_mh', $         
             '(*!dustem_params).grains(2).mdust_o_mh' $
              ]                   


iv = [0.17E-02,0.63E-03,0.51E-02]

Npar=n_elements(pd)  

ulimed=replicate(0,Npar)

llimed=replicate(1,Npar)
llims=replicate(0,Npar)

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

dustem_init_plugins,pd

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


;=== RUN fit
tol=1.e-6
xtol=1.e-6
Nitermax=30; just an initialization

xr=[1.,5e5]
yr=[1e-2,1.00e8]

loadct,13

title=textoidl('TEST_SED')
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
t1=systime(0,/sec)
res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,xtit=xtit,ytit=ytit,title=title)
t2=systime(0,/sec)


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

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

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

loadct,13

dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=title,yr=yr,xr=xr,/ysty,/xsty,res=res,chi2=chi2,rchi2=rchi2,/xlog,/ylog,/pol,ps=ps,png=png

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

the_end:

END