checkvar, ps, 0
initialize, ps

restgenx, file='ospex_results_april21_leah.geny', results

;emission measure in cm^-3
EMa = results.spex_summ_params[0,*]
em = (EMa)
em = reform(em)

;temp in K
Ta = results.spex_summ_params[1,*]
;from keV to MK
t = (Ta) * 11.6
t = reform(t)

;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 = 0 & e = n_elements(tarr) - 1
stime = anytim(tarr(s),/vms)
etime = anytim(tarr(e),/vms)

utb = min(tarr)

v = find_v('c:\documents and settings\lhaga\my documents\2002 april 21', filenames = 'flux_output_april21_6-12_c50.txt',tarr )

; Boltzmann constant - k = 1.38e-16 erg K^-1
   k = 1.38e-16

; i.e., U = 4.14e-16 T(K) sqrt(EM(cm^-3) * V(cm^-3)) ergs
; U = 4.14e-16 tempr * 1.e6 (K) sqrt(emis*1.e49(cm^-3) * V(cm^3)) ergs
; U = 4.14 e-10 * sqrt(10) * 1.e24 *tempr * sqrt(emis * V) ergs
; U = 1.309e15 * tempr * sqrt(emis * V) erg
; V = source volume in cm^3

   U = 1.309e15 * t * sqrt(em * V)

bcl = 11
Usmooth = smooth(U,bcl)
Udot = deriv(tarr,Usmooth)
Udot = smooth(Udot,bcl)

; Plot radiative loss rate
Lrad = loss_rad_chianti(em,t*1.e6)
utplot,tarr(s:e)-utb,Lrad(s:e),utb, $
   ytitle='Energy loss rate in ergs s!u-1!n', title = 'April 21, 2002    Red = conductive loss,  Black = radiative loss,  Pink = enthlapy',$
   /ylog,position=[0.15,0.10,0.95,0.90], yrange=[1e23,1e32]

; Plot conductive loss rate
; Loop cross-sectional area in cm^2
; Loop length in cm
area = 1.2e19
length = 0.5e9

Lcond = loss_cond(area,length,t*1.e6)
outplot,tarr(s:e)-utb,Lcond(s:e), utb, color=2
;ytitle='Conductive energy loss rate in ergs s!u-1!n', $

a = find_a('c:\documents and settings\lhaga\my documents\2002 april 21', filenames = 'flux_output_april21_6-12_c50.txt',tarr )

terma = sqrt(em*v)
termb = sqrt(em*v)
last_term = (terma-termb)/20

;equation from SMM in ergs/sec
Lenthalpy = 2.5 * a * Ks * T * 1.e6 * last_term

;equation from antonucci et al. paper
;Lenthalpy = 5 * Kb * sqrt((EM*1.e49) / V) * T * 1.e6 * 2.e7 * 2 * a

outplot,tarr(s:e)-utb,Lenthalpy(s:e), utb, color=12

if ps gt 0 then device,/close
set_plot,'win'

end