; Calcuate Energy Input to Thermal Plasma

; Initialize - 0 for windows, 1 for postscript (portrait), 2 for postscript (landscape)
; type ps=0, 1, or 2 on command line
checkvar, ps, 0
initialize, ps

!p.charsize=1.9
!p.multi = [0,1,5]

restgenx, file='ospex_results_oct28_1therm_pileup7.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 = 12. & e = n_elements(tarr) - 7.
stime = anytim(tarr(s),/vms)
etime = anytim(tarr(e),/vms)

utb = min(tarr)
; Plot emission measure
utplot,tarr(s:e)-utb,em(s:e),utb, $
   /ylog,ytitle='EM (10!u49!n cm!u-3!n)', $
   color=0,background=255, xticklen = 1, xtickname = ' ', xstyle = 1,$
   xtitle = ' ', title = 'Plots from RHESSI data  fit with one thermal component', ymargin = [0,2], yrange = [10.,100.]
;xyouts,0.3,0.88,'Start Time = ' + stime, color=0,/normal
;xyouts,0.3,0.85,'End Time = ' + etime, color=0,/normal

; Plot temperature in MK
utplot,tarr(s:e)-utb,t(s:e),utb, $
ytitle='Temperature (MK)', $
   color=0,background=255,xticklen = 1,$
   xtitle = ' ', ymargin = [0,0],yrange = [20., 40.], xstyle = 1

;-----------------------------------------------------------
v = find_v('c:\documents and settings\lhaga\my documents\2003 october 28', filenames = 'flux_output_oct28_12-25.txt',tarr )

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




;Plot thermal energy U = 3 n k T V erg
; 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)

utplot,tarr(s:e)-utb,U(s:e), utb, $
/ylog,ytitle='Thermal E (ergs)', $
   color=0,background=255,xticklen = 1,$
   xtitle = ' ', ymargin = [0,0], yrange = [5.e31,1.e32], xstyle = 1



;y = findgen(15600,2)
;A = [1.2e31,1.3e30]
;w = [5200.,20000.]

;y1 = A(0) * exp(-tarray/w(0))
;y2 = A(1) * exp(-tarray/w(1))

;outplot,tarray,y1,color=0
;outplot,tarray,y2,color=0

;xyouts,0.4,0.85,charsize=1.5, $
;   'Zero Intercepts = '+num2str(A(*)), /normal,color=0

;xyouts,0.4,0.75,charsize=1.5, $
;   'Time Constants (s) = '+num2str(t(*)), /normal,color=0

; Plot dU/dt - rate of change of thermal energy with time
; First, smooth energy array, U, with a boxcar of length, bcl = an odd number of time intervals
bcl = 5
Usmooth = smooth(U,bcl)
Udot = deriv(tarr,Usmooth)
Udot = smooth(Udot,bcl)


;utplot,tarr(s:e)-utb,Udot(s:e),utb, $
;   ytitle='Udot', $
;   color=0,background=255,xticklen = 1,xtitle = ' ', $
;   ymargin = [0,0],yrange=[-4.e28,3.e28], xstyle = 1

l=0
outplot, tarr(s:e)-utb, udot(s:e)*l, utb, color = 10

; Plot radiative loss rate
Lrad = loss_rad_chianti(em,t*1.e6)
utplot,tarr(s:e)-utb,Lrad(s:e),utb, color = 0, $
   ytitle='E loss rate (ergs/s)', $
   /ylog,yrange = [1.e25,1.e29], $
   background=255,xticklen = 1,$
   xtitle = ' ', ymargin = [0,0], xstyle = 1

;LradI = Ilrad(U, Lrad)
;LradIP = Ilrad_peak(lradi,u)
;LradIM = Ilrad_max(lradi)


; Plot conductive loss rate
; Loop cross-sectional area in cm^2
; Loop length in cm
area =  4.e18
length = 6.e9

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', $

;   ,background=255, ymargin = [4,0], yrange = [1.e22,1.e24], /ylog, $
;   xticklen = 1,xtitle = ' '

;xyouts,0.2,0.20, $
;'Area = ' + num2str(area) + ' cm!u2!n, Length = ' + num2str(length) + ' cm', $
;   color=0,/normal

xyouts,0.4,0.3, $
'black=lrad, red=lcond, green=lrad+lcond+udot', $
   color=0,/normal,charsize=1.5

; Plot rate of total energy input to thermal plasma
; Pin = Udot + Lrad + Lcond erg s^-1

    Pinn = Udot + Lrad + Lcond &$
    pin = pinn > 0

outplot,tarr(s:e)-utb,Pin(s:e),utb, color=7
;ytitle='Conductive energy loss rate in ergs s!u-1!n', $

; Integrate Pin with wrt time to get the total energy input to the thermal plasma

;IPin = integrate(Pin)

jmax = N_ELEMENTS(U)

Ipin = fltarr(jmax)
Ipin[0] = U[0]
FOR j = 1, jmax-2 DO BEGIN
Ipin[j] = Ipin[j-1] + (udot[j-1] + lcond[j-1] + lrad[j-1]) * 20.
ENDFOR

utplot, tarr(s:e)-utb,IPin(s:e),utb, ytitle='E Input (ergs)', $
   color=0,background=255, ymargin = [4,0], yrange = [3.e31,4.e31],$
   xticklen = 1, xstyle = 1;,xtitle = ' '

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

end