PRO dustem_fit_sed_orion,ps=ps,png=png

dustem_define_la_common

;=== initialise dustem

use_mode='G17_MODELC';'DBP90';'G17_MODELA';'THEMIS';;'USER_MODEL';'THEMIS';;'THEMIS';'G17_MODELA';'THEMIS';'NY_MODELA';'G17_MODELA';'DBP90';'USER_MODEL';'THEMIS';'NY_MODELA';'G17_MODELA';'DBP90';'COMPIEGNE_ETAL2010'
;,/pol
dustem_init,mode=use_mode,kwords=['logn-chrg-spin','plaw-pol','plaw-pol'],/pol;,kwords=['plaw-ed','logn','logn'];,/pol;kwords=['plaw-ed-chrg-spin'],/pol;,kwords=['plaw-ed'];kwords=['plaw-ed-chrg-spin'];,/pol;
;
;stop
!dustem_verbose=0
!dustem_show_plot=1
!dustem_nocatch=1

;specify level
lev=['300.000','400.000','500.000','600.000','800.000','1000.00','1200.00','1500.00','3000.00','5000.00','10000.0','15000.0']
indxlev=6 ; meaning level 400, 1 is for 600 etc...;THIS IS FOR THE PLOTTING OF ONE LEVEL ONLY
msk=['mask1','mask2','mask3'] ;this is used to be able to plot the new data that JP generated.
;1,2,3 correspond to the North-West, Center and South-East part of the OrionA integral shape filament
msk=msk(1) ;mask2 center of ORION_A 
plckprd=['PR3','PR4','SROLL2']
plckprd=plckprd(1) ;PR4
polangmatch=['ang5','ang10','ang15']
polangmatch=polangmatch(1) ;best match. Lowest difference.
co_corr = ['noCO0','noCO1','noCO2']
co_corr = co_corr(2)

;opening corresponding sed file
smooth_str='_smoothedLFI3'

; dir_seds='/Users/ilyeschoubani/Downloads/SC110012_V2p2_PRO_5000001211_2/Orion/seds/' ;new SEDS have a problem
; ;dir_seds='/Users/ilyeschoubani/Downloads/SC110012_V2p2_PRO_5000001211/Orion/seds/'
dir_seds = '/Users/ilyeschoubani/Downloads/SC110012_V2p2_PRO_5000001210/Orion/seds/'

;with the mask - new batch still has the error
;file_sed=dir_seds+'extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'_'+msk+'_'+co_corr+'.xcat'
;without the mask - new batch does not have the error
file_sed=dir_seds+'extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'_'+co_corr+'.xcat'
;;THIS WAS PRIOR;file_sed=dir_seds+'extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'.xcat'

data_struct=read_xcat(file_sed,/silent) ;This is the new general structure


;stop
;
;===reading and applying the HCD:
openr,unit,file_sed,/get_lun
str=''
WHILE not eof(unit) DO BEGIN
  readf,unit,str
  if strtrim(strmid(str,13,3),2) EQ '857' then *!dustem_HCD=float(strtrim(strmid(str,21),2))/5.25*1.00E+26  
ENDWHILE
close,unit
free_lun,unit
;===================

;stop
;generating fake errors for SED and Q/USED - here 10%
;stop

data_struct.sigmaII = (double((data_struct.StokesI*0.1)^2)) > 1d-20 ;this is false. Squaring it does not make sense. The square root is the one that should be taken 

;NOT THE REAL REASON NOW;because of a bigger error in the variance computation in the xcat file. CHIPASS has an error value for StokesQ and StokesU but not for their associated flux.
;CHIPASS filter isn't well defined
;I looked it up and found a bandpass of 64 MHz (in 1024 channels).
data_struct(-1).sigmauu = la_undef()
data_struct(-1).sigmaqq = la_undef()
data_struct(-1).sigmaii = la_undef()
data_struct(-1).stokesu = la_undef()
data_struct(-1).stokesq = la_undef()
data_struct(-1).stokesi = la_undef()

data_struct(-1).psi = la_undef()
data_struct(-1).sigma_psi = la_undef()

; data_struct(-1).stokesi = la_undef()
; tt=where(data_struct.sigmauu EQ -39.8637)
; ;stop
; 
; data_struct[tt].sigmauu = la_undef()
; data_struct[tt].sigmaqq = la_undef()

;data_struct.sigma_smallp = (double((data_struct.smallp*0.1)^2)) > 1d-20
idsigq=where(data_struct.sigmaQQ NE la_undef())
idsigu=where(data_struct.sigmaUU NE la_undef())
;stop
;......
;DON'T FORGET TO UNCOMMENT THIS.  

;data_struct[idsigq].sigmaQQ = (double((data_struct[idsigq].StokesQ*0.1)^2)) > 1d-20

;data_struct[idsigu].sigmaUU = (double((data_struct[idsigu].StokesU*0.1)^2)) > 1d-20

;#######FOR THE FITTING OF Q AND U
;locating the negative values in stokesQ data so that filters are omitted from the fit
indq=where(data_struct.stokesQ gt 0 and data_struct.stokesQ ne la_undef())
rmvq=(data_struct(indq).filter)

;locating the negative values in StokesU data so that filters are omitted form the fit
indu=where(data_struct.stokesU lt 0 and data_struct.stokesU ne la_undef())
rmvu=(data_struct(indu).filter)

rmvpsi=[rmvq,rmvu]
;removing LFI2 manually 
remove,1,rmvpsi

;stop
;BECAUSE SPASS band has no errors:
;NOW IT DOES
; (data_struct(-2).sigma_largeP)=(data_struct(-3).sigma_largeP)
; (data_struct(-2).sigma_smallp)=(data_struct(-3).sigma_smallp)
; (data_struct(-2).sigmaqq)=(data_struct(-3).sigmaqq)
; (data_struct(-2).sigmauu)=(data_struct(-3).sigmauu)
; (data_struct(-2).sigma_psi)=(data_struct(-3).sigma_psi)

;stop
fit_struct=data_struct
;stop
;uncomment this block afterwards - because you are now just fitting the sed.

; fit_struct=dustem_mask_data(fit_struct,omit=list(rmvq,rmvu,rmvpsi,rmvpsi,rmvpsi),data_set=['QSED','USED','PSI','POLSED','POLFRAC']) ;works
; rmvq=(data_struct(-2).filter)
; 
; 
; fit_struct=dustem_mask_data(fit_struct,omit=list(rmvq,rmvq,rmvq,rmvq,rmvq),data_set=['QSED','USED','PSI','POLSED','POLFRAC'])
; 
; 

; fit_struct=dustem_mask_data(fit_struct,omit=rmvq,data_set=['QSED'])
; fit_struct=dustem_mask_data(fit_struct,omit=rmvu,data_set=['USED'])

;stop
;since we're using q and u I will have to null polsed and polfrac. Otherwise set_data will return an error
;fit_struct.largeP=la_undef()
;fit_struct.smallp=la_undef()

;doing the same for the errors as dustem_check will notice some erros have undefined associated data points.
;fit_struct.sigma_largeP=la_undef()
;fit_struct.sigma_smallp=la_undef()

;because I am now fitting only the sed
; fit_struct.sigmaQQ = la_undef()
; fit_struct.sigmaUU = la_undef()
;stop
;show structure
show_st=data_struct
;stop
; show_st.largeP=la_undef()
; show_st.sigma_largeP=la_undef()
; show_st.stokesQ=la_undef()
; show_st.sigmaQQ=la_undef()
; show_st.stokesU=la_undef()
; show_st.sigmaUU=la_undef()

;#################################################

;params_TH=[5.44E-05,5.60E-04,1.22E-04,3.0,4.13E+06,2.25E+00];,0.30,53]

params_G17C=[8.85E-06,1.25E-05,5.30E-03,100.0,3.00E3,6.43,1.00E+07,2.73E+00,58]

;params_DBP=[1.06E-05,7.00E-04,2.32E-03,6.43,1.00E+07,2.57E+00]

;=== 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', $
             ;'(*!dustem_params).grains(3).mdust_o_mh', $
             ;'(*!dustem_params).grains(4).mdust_o_mh', $
             '(*!dustem_params).gas.Tgas',$
             '(*!dustem_params).gas.nh',$
             'dustem_plugin_stellar_population_O6V3',$      
             ;'dustem_plugin_mbbdy_isrf_1',$ 
             ;'dustem_plugin_mbbdy_isrf_2', $
             ;'dustem_plugin_mbbdy_1',$ 
             ;'dustem_plugin_mbbdy_2', $
             ;'dustem_plugin_mbbdy_3', $
             ;'dustem_plugin_mbbdy_isrf_3', $
             ;'dustem_plugin_mbbdy_isrf_4', $
             ;'dustem_plugin_synchrotron_1', $                 
             ;'dustem_plugin_synchrotron_2', $
             ;'dustem_plugin_synchrotron_3', $                 
             ;'dustem_plugin_synchrotron_4', $
             'dustem_plugin_freefree_1', $      ;4--->6           
             'dustem_plugin_freefree_2', $ ;5
             ;'dustem_plugin_freefree_3', $
             ;'dustem_plugin_modify_dust_pol_1',$
             'dustem_plugin_modify_dust_pol_2'$
             ;'dustem_plugin_mdfy_spinning_pol_1'$
             ;'dustem_plugin_continuum_2', $                 
             ;'dustem_plugin_continuum_1' $
            ;'dustem_plugin_mdfy_spinning_pol_1', $
            ;'dustem_plugin_mdfy_spinning_pol_2' $
              ]                   


iv=params_G17C

Npar=n_elements(pd)  

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

ulimed[6]=1

;uncomment
;ulimed[6]=1

;ulimed[9]=1
ulims[6]=1.00E+07
;ulims[8]=5.0E-03

;uncomment
;ulims[6]=9.00E-01;becomes 9

;stop

llimed=replicate(1,Npar);[1]
;llims=[1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20];,1.00e-20,1.00e-20]
llims=replicate(5.3e-6,Npar);[10]

;llims[0]=1.0e-9
;llims[3]=1.0
   
llims[6]=1.00E+04
llims[7]=1.00E-04

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

;   fpd = [ $
;                 '(*!dustem_params).grains(1).mdust_o_mh'$         
;                 ;'(*!dustem_params).grains(1).mdust_o_mh'$ 
;                 ;'(*!dustem_params).grains(4).mdust_o_mh'$   
;                 ]                   
;   fiv = [1.00E-20] 
; 

; dustem_init_fixed_params, fpd, fiv

;=== SET THE OBSERVATION STRUCTURE

ind=where(data_struct.wave GE 0.1)

st_ie=dustem_set_data(st_sed_fit=fit_struct[ind],st_sed_show=show_st[ind],rchi2_weight=rchi2_weight,f_HI=f_HI) ; comment if second line is to be used
;stop

;=== RUN fit
tol=1.e-10
xtol=1.e-10
Nitermax=3;0; just an initialization

xr=[1.,5e5]
yr=[5e-8,1.00e6]

loadct,13
;!y.range=[1e-8,100] ;This is to adjust plot range from outside the routine
yyy="restrng=textoidl('[MJy/sr]_{"+strmid(smooth_str,9)+'|'+plckprd+"}')"
yyy=execute(yyy)
title=textoidl('OrionA I_{240\mum} > ')+strtrim(lev[indxlev],2)+restrng

sttr='          I_{\nu} (MJy/sr) for N_{H}='+strmid(string(*!dustem_HCD,format='(1E10.2)'),2)+'H/cm^2'
ytit=textoidl(sttr)

;xtit=''
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,_extra=_extra);
t2=systime(0,/sec)

;help,(*!dustem_fit),/STRUCTURE


;=== SAVE FIT RESULTS

;NB: THIS LINE IS WRONG BECAUSE I AM NOT USING THE ANGLE CORRELATION CRITERIA
sedid='extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'_masked_'+polangmatch+'_'+msk
file='/tmp/DUSTEM_lastfit_orionA_'+sedid+'.sav'
dustem_save_system_variables,file
message,'Saved '+file,/continue

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

;=== Plot best fit

; win=1 ; who wrote this? 
; 
; window,win & win=win+1


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

dustemwrap_plot,res,stp,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=title
;print, 'DONE!!! IT WORKED'


; 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
; 
; ps_file=!dustem_wrap_soft_dir+'Docs/'+'Last_dustem_fit.ps'
; 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_file,png=png,fpol=fpol


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

the_end:

END