;Program to use results obtained with SPEX. (Sui et al, Spectroscopy 2nd steps)

;Assumed following command has been issued in a SPEX session:
; idl, spex_save_data , /all, 'spex_2mar14_0137-0152_6-250_keV.sav'

restore,'spex_2mar14_0137-0152_6-250_keV.sav' 
!p.multi=[0,6,5]

loadct, 2  & red=92 & green=30 & blue=220 & white=255 ;4
if (!d.name NE 'PS') then window, 10, xsize = 1200, ysize = 1000
for interval=0,29 do begin
 
;Compute the photon flux (=(obsi-backi)/convi)
 phflux = (obsi-backi)/convi
 
;Compute the error of the photon flux (=sqrt(eobsi^2 + ebacki^2) /convi):
errflux = sqrt(eobsi^2 + ebacki^2) /convi
 
;Compute the background flux (=backi/convi)
backflux = backi/convi
 
;Compute the mean value of the energy edges for each bin
; print, edges[*, 0:2]
edge_products, edges, mean = eph
print, eph[0:2]
 
;Get the fitting parameters for the first spectral fit:
fitpara = apar_arr[*, interval]
print, fitpara
 
;Check the UT time of that interval and compare to the time history window
uttime = anytim(xselect[*, interval] + utbase, /ecs)
print, uttime

;Compute the thermal photon flux:
thermflux = f_vth(edges, fitpara[0:1])
help, thermflux
 
;Set the cutoff energy and low-energy slope that were used for the power-law fit but not saved
a_cutoff = [1, 1.5]
 
;Add these paramters to the array of the other fit parameters that were saved.
fitpara = [fitpara, a_cutoff]
print, fitpara
 
;Calculate the nonthermal flux
nonthermflux = f_bpow(eph, fitpara[2:7])
help, nonthermflux
 
;Plot the spectrum
 ti=uttime[0]
;Plot the thermal spectrum
plot, eph, thermflux, /xlog, /ylog, xtitle = 'Photon Energy (keV)', $
ytitle = 'Flux (photons s!u-1!n cm!u-2!n keV!u-1!n)', xrange = [10.0, 100.], $
yrange=[1.0e-2, 1.0e6], xstyle = 1, ystyle = 1, /nodata,TIT=TI
oplot, eph, thermflux, lineStyle = 1, color = blue, thick = 4
 
;Overplot the nonthermal spectrum
oplot, eph, nonthermFlux, lineStyle = 0, thick = 2, color = green
;Overplot the total model spectrum
oplot, eph, thermflux + nonthermflux, linestyle = 0, color =white , thick = 1; 2;4
;Overplot the measured photon fluxes for the same interval
oplot, eph, phflux[*, interval], psym = 1, color = blue, symsize = 0.2, thick = 1;2
;Overplot the error bars on the measured flux
oplot_err, eph, phflux[*, interval], yerr = errflux[*, interval], psym = 1, symsize = 1, color = white ; 144
;Overplot the background spectrum
oplot, eph, backflux[*, 0], psym = 5, symsize = 0.5, color = red ;208
 
endfor ; interval

end