; 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,6]

restore,'idlsave_goes020421_long_interval.dat'

; Calculate start time

; Time range to plot
s = 0. & e = n_elements(tarray) - 1
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='Emission Measure in 10!u49!n cm!u-3!n', title = 'plots from GOES data ',$
   color=0,background=255, xticklen = 1, xtickname = ' ',$
   xtitle = ' ', ymargin = [0,2], yrange = [0.001,100.]

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 in MK', $
   color=0,background=255,xticklen = 1,$
   xtitle = ' ', ymargin = [0,0],yrange = [5, 20]

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

c = rd_ascii('flux_output_april21_12-25_c50_changes_removed_noholes.txt')
c = c[1:*]
b = str2cols(c)
;b = strcompress(b,/remove_all)
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)) + 2.8e28

for i=0, n_elements(areas)-1 do begin
	i1 = value_locate(vec, stime_sec[i])
	i2 = value_locate(vec, etime_sec[i])
;	print, i1,i2, im_vol[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

   U = 1.309e15 * tempr * sqrt(emis * V)


utplot,tarray[s:e],U[s:e],utbase, $
   /ylog,ytitle='Thermal Energy U in ergs', $
   color=0,background=255,xticklen = 1,$
   xtitle = ' ', ymargin = [0,0], yrange=[1.e28,1.e32]

;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 = 51
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=[-4.e28,4.e28]

; Plot radiative loss rate

Lrad = loss_rad_chianti(emis,tempr*1.e6)
utplot,tarray(s:e),Lrad(s:e),color = 0, $
   ytitle='Energy loss rate in ergs s!u-1!n', $
   /ylog,yrange = [1.e21,1.e29], $
   background=255,xticklen = 1,$
   xtitle = ' ', ymargin = [0,0]
; Plot conductive loss rate
; Loop cross-sectional area in cm^2
; Loop length in cm
area = 12.e18
length = 3.6e10

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

; 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.
Ipin[j] = Ipin[j-1] + (lrad[j-1]) * 3.
ENDFOR


;utplot,tarray(s:e),IPin(s:e),ytitle='Energy Input to Thermal Plasma in ergs', $
utplot,tarray(s:e),IPin(s:e),ytitle='Integrated Radiated Energy in ergs', $
   color=0,background=255, ymargin = [4,0],  yrange = [1.e28,1.e32], $
   xticklen = 1;,xtitle = ' '

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

end