checkvar,ps,0
initialize, ps

restgenx, file='ospex_results_oct23.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.

;-----------------------------------------------------------
v = find_v_oct23('C:\Documents and Settings\lhaga\My Documents\2003 October 23', filenames = 'flux_output_oct23_12-25_d45678_changed.txt',tarr )

;----------------------------------------------------------

; nV = SQRT(EM * V)
; 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)


jmax = N_ELEMENTS(U)

;********** UVCS & RHESSI DATA POINTS **********

;number of UVCS observations
iuvcs = 4

;UVCS start times in seconds from utbase
uvcs_t = strarr(iuvcs)
uvcs_ts = fltarr(iuvcs)
uvcs_t[0] = '23-Oct-2003 19:52:16.000'
uvcs_t[1] = '23-Oct-2003 19:54:27.000'
uvcs_t[2] = '23-Oct-2003 19:56:37.000'
uvcs_t[3] = '23-Oct-2003 19:58:46.000'
uvcs_ts = anytim(uvcs_t, /sec)

;computed transition region energy fluxes (erg/sec) from UVCS
uvcs_f = fltarr(iuvcs)
uvcs_f[0] = 0.404d+26
uvcs_f[1] = 2.328d+26
uvcs_f[2] = 2.956d+26
uvcs_f[3] = 2.676d+26

;find indices for RHESSI times corresponding to UVCS time intervals
rhessi_ind_l = fltarr(iuvcs)
rhessi_ind_h = fltarr(iuvcs)
rhessi_ind_l = value_locate(tarr, uvcs_ts)
rhessi_ind_h = value_locate(tarr, uvcs_ts + 120.)

;create array containing UVCS fluxes at GOES times
uvcs_flux = fltarr(jmax)
FOR i = 0, iuvcs-1 DO BEGIN
; if (i LE 2) then begin
;   indarr2 = where( (tarr GE uvcs_ts[i]+120) and (tarr LT uvcs_ts[i+1]) )
;    uvcs_flux[indarr2] = uvcs_f[i]
;endif
 indarr = WHERE( (tarr GE uvcs_ts[i]) AND (tarr LT uvcs_ts[i]+120.) )
 uvcs_flux[indarr] = uvcs_f[i]
ENDFOR

;radiative loss in ergs/s
Lrad = loss_rad_chianti(em,t*1.e6)

;compute the energy flux at each time tgoes + lrad  2 point derivative
udot = (shift(u,-1)-u)/(shift(tarr,-1)-tarr)

;therm_ein = integrate(udot)
;compute the time evolution of the energy input to the thermal plasma
therm_ein = fltarr(jmax)
therm_ein[0] = U[0]
FOR j = 1, jmax-2 DO BEGIN
therm_ein[j] = therm_ein[j-1] + (udot[j-1] + uvcs_flux[j-1]) * (tarr[j] - tarr[j-1])
ENDFOR

;********** Plot Results **********

;setprintversion
;pson, file = 'goes_th_energy_0723.ps', paper = 'LETTER'

;time limits for utplot
tmin = '23-Jul-2002 00:20:30.000'
tmin = anytim(tmin, /sec)
tmax = '23-Jul-2002 00:36:00.000'
tmax = anytim(tmax, /sec)

utb = min(tarr)
utplot, tarr - utb, U, utb, /ylog, xstyle = 1, ytitle = 'Energy (ergs)', $
    xtitle = 'Time (UT)', title = 'July 23, 2002   RHESSI and UVCS data',  labelpar = [10, 17], $
    nticks = 3, timerange = [tmin,tmax] ,$
    yrange = [1.0e26, 3.0e31], position = [0.09, 0.1, 0.95, 0.9],  /nodata, ystyle = 1, charsize = 1.2

outplot, tarr-utb, u, linestyle = 1, thick = 2

;Plots black line with Eloss added to U
;overplot the total energy injected into the thermal plasma
outplot, tarr-utb,therm_ein, linestyle = 0, thick = 2, color = 10

FOR ip = 0, iuvcs-1 DO BEGIN

;bottom lines
;overplot the UVCS fluxes
PlotS, [uvcs_ts[ip],uvcs_ts[ip]+120.]-utb, [uvcs_f[ip],uvcs_f[ip]], thick = 3

ENDFOR

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

;psoff
;clearprintversion

END