; 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.multi = [0,1,5]

!p.charsize = 1.9

restore,'idlsave_goes031023.dat'

; Calculate start time

; Time range to plot
s = 300. & e = n_elements(tarray) - 200
stime = anytim(utbase+tarray(s),/atime)
etime = anytim(utbase+tarray(e),/atime)

; Plot emission measure
utplot,tarray(s:e),emis(s:e),utbase, $
   /ylog, ytitle='EMeasure (10!u49!n cm!u-3!n)', title = ' Plot from GOES data ', $
   color=0,background=255, xticklen = 1, xtickname = ' ',$
   xtitle = ' ', ymargin = [0,2.], yrange = [0.01,15.], xstyle = 1

;xyouts,0.3,0.95,'Start Time = ' + stime,color=0,/normal
;xyouts,0.3,0.85,'End Time = ' + etime,color=0,/normal

; Plot temperature in MK
utplot,tarray(s:e),tempr(s:e),utbase, $
   ytitle='Temperature (MK)', $
   color=0,background=255,xticklen = 1,$
   xtitle = ' ', ymargin = [0,0], yrange = [5, 25], xstyle = 1

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

c = rd_ascii('flux_output_oct23_12-25_d45678.txt')
c = c[1:*]
b = str2cols(c)
box = fix(reform(b[0,*]))
areas = float(reform(b[5,*]))
q = where(box eq 0)
b0 = b[*,q]
box = (box[q])
areas = (areas[q])
box = fix(reform(b0[0,*]))
areas = float(reform(b0[5,*]))
statimes = reform(b0[1,*] + ' ' + b0[2,*])
endtimes=reform(b0[1,*] + ' ' + b0[3,*])
stime_sec = anytim(reform(statimes))
etime_sec = anytim(reform(endtimes))

goestimes = tarray+utbase
vec = goestimes
vol = (4./3.) * (1. / sqrt (!pi))
im_vol = vol * areas^(1.5)*72498958.3^3

v = fltarr(n_elements(tarray)) + 8.8e27

for i=0, n_elements(areas)-1 do begin
	i1 = value_locate(vec, stime_sec[i])
	i2 = value_locate(vec, etime_sec[i])

	 v[i1:i2] = im_vol[i]

endfor


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


;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
;reads ascii file and returns it as a text array

bcl=51
vsmooth = smooth(v,bcl)
   U = 1.309e15 * tempr * sqrt(emis * Vsmooth)


utplot,tarray[s:e],U[s:e],utbase, $
   /ylog,ytitle='Thermal E (ergs)', $
   color=0,background=255,xticklen = 1,$
   xtitle = ' ', ymargin = [0,0], yrange=[1e29,1e31], xstyle = 1

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

y1 = A(0) * exp(-tarray/t(0))
y2 = A(1) * exp(-tarray/t(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 = 100
Usmooth = smooth(u,bcl)
Udot = deriv(tarray,Usmooth)
Udot = smooth(Udot,bcl)
;utplot,tarray(s:e),Udot(s:e),ytitle='Udot', $
;   color=0,background=255,xticklen = 1,xtitle = ' ', $
;   ymargin = [0,0],yrange=[-1.e28,2.e28], xstyle = 1

; Plot radiative loss rate

Lrad = loss_rad_chianti(emis,tempr*1.e6)
utplot,tarray(s:e),Lrad(s:e),color = 0, $
   ytitle='E loss rate (ergs/s)', $
   /ylog,yrange = [1.e22,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 = 6.3e18
length = 2.1e9

Lcond = loss_cond(area,length,tempr*1.e6)
outplot,tarray(s:e),Lcond(s:e),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.23, $
'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

   Pin = Udot + Lrad + Lcond
outplot,tarray(s:e),Pin(s:e),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]) * 3.
ENDFOR


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

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

end