;plots thermal energy vs time

restgenx, file='ospex_results_nov3_multitherm.geny', results

a0 = results.spex_summ_params[15,*]
a0 = reform(a0)
aone = results.spex_summ_params[16,*]
aone = reform(aone)
a2 = results.spex_summ_params[17,*]
a2 = reform(a2)
a3 = results.spex_summ_params[18,*]
a3 = reform(a3)

;start and end time in sec
start_time = results.spex_summ_time_interval[0,*]
start_time = reform(start_time)
end_time = results.spex_summ_time_interval[1,*]
end_time = reform(end_time)
tarr = (start_time + end_time)/2.

; Time range to plot
s = 3 & e = n_elements(tarr) - 1
stime = anytim(tarr(s),/vms)
etime = anytim(tarr(e),/vms)

utb = min(tarr)

;runs program find_v to find the volume
v = find_v('c:\documents and settings\lhaga\my documents\2003 november 3 hour9', filenames = 'flux_output_nov3_12-25_d345678_c50_changed.txt',tarr )

; plot thermal energy
; U = 3 n k T V erg
; nV = SQRT(EM * V)
; n = 10^11cm^-3   n10 = n/10^10
; DEM 10^49 cm^-3 keV^-1
; erg/keV = 1.6022*10^-9 erg/keV
; k = 1.38e-16

U = 1.92 * 10.^31 * (a0 / (a3-2)) * ((aone/2.)^(2-a3) - (a2/2.)^(2-a3)) / 10.

utplot,tarr(s:e)-utb,U(s:e), utb, $
/ylog,ytitle='Thermal Energy U in ergs', yrange = [10.e28, 10e32]

end