pro energy_input_102803_48

; Calcuate Energy Input to Thermal Plasma

; Input

	;Save file from GOES Workbench
   	path = 'C:\Documents and Settings\lhaga\My Documents\GOES Files for Brian\10_28\'
	file = 'idlsave_goes031028_48.dat'
	restore, 'idlsave_goes031028_48.dat'

		; Time range to plot
	plot_start_time = min(tarray)+utbase
	plot_start_time = anytim(plot_start_time,/atime)
	plot_end_time = max(tarray)+utbase
	plot_end_time = anytim(plot_end_time,/atime)

	;Define flare start and end times in yymmdd hh:mm:ss.sss
	flare_start_time = min(tarray)+utbase
	flare_start_time = anytim(flare_start_time,/atime)
	flare_end_time = max(tarray)+utbase
	flare_end_time = anytim(flare_end_time,/atime)


	;Define time for preflare background
	back_pre = tarray[10]

	A = 6.e18 ; Source area in cm^2
	V = A^(3./2.) ; Source volume in cm^3

	;Output file name written to the working directory
	out_file = 'GOES_28OCT2003_48.txt'

	;Set ps = 0 for windows output, 1 for postscript (portrait), 2 for postscript (landscape)
	ps=0

;End Input

;Save current plot parameters and revert to the default values
	old_device = !d.name
	store_plotvar
	cleanplot,/silent


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

;!p.multi = [remaining,columns,rows, , ]
!p.multi = [0,1,5]
;Make X-axis scale exactly match range of times
!x.style = 1
!p.charsize = 2.0
DEVICE, SET_FONT='Times', /TT_FONT

x_pixels = 600 & y_pixels = 700
if ps eq 0 then window,10, title = 'Thermal Parameters vs. Time', xsize = x_pixels, ysize = y_pixels

restore,path + file

; UTBASE is the time in seconds since 1-jan-1979
; tarray has the interval start times in seconds of day after UTBASE

delta_t = tarray[1] - tarray[0]
n_tarray = n_elements(tarray) - 1

;Calculate index of plot start and end times.
s = calc_t_index(plot_start_time, tarray, t_in_s, delta_t, utbase)
if s lt 0 or s gt n_tarray then s = 0
e = calc_t_index(plot_end_time, tarray, t_in_s, delta_t, utbase)
if e lt 0 or e gt n_tarray then e = n_tarray
trange = [UTBASE+tarray(s), UTBASE+tarray(e)]
!x.range = trange

stime = anytim(utbase+tarray(s),/atime)
etime = anytim(utbase+tarray(e),/atime)

;Flare start and end times in yymmdd hh:mm:ss.sss
; and in the index for tarray - fsti and feti.

fsti = calc_t_index(flare_start_time, tarray, t_in_s, delta_t, utbase)
if fsti lt 0 or fsti gt n_tarray then fsti = 0
feti = calc_t_index(flare_end_time, tarray, t_in_s, delta_t, utbase)
if feti lt 0 or feti gt n_tarray then feti = n_tarray

;Time for preflare background - back-pre
backi = calc_t_index(back_pre,tarray,t_in_s,delta_t,utbase)

;Define blank label for X-axis
xx=[' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ']

;Plot GOES ch. 1 & 2 fluxes in ergs s^-1 from watts m^-2
;Distance to Sun (1 AU) = 1.496 E+13 cm
au_in_cm = 1.496e+13
; 1 watt = 1 joule per second = 10^7 erg s^-1
; 1 watt m^-2 = (2.0 * !pi * au**2 * 1.e7) / 1.e4 erg s^-1
wpsm2eps = 2.0 * !pi * au_in_cm^2 * 1.e3

flux = yclean * wpsm2eps

ymax = max(yclean[s:e,*],min = ymin)
if ymin lt 1.e-8 then ymin = 1.e-8
if ymax gt 1.e-2 then ymax = 1.e-2
utplot, tarray[s:e],yclean[s:e,0], utbase, $
	/ylog, ytitle='Flux in watt m!u-2!n', $
	color=0,background=255, xticklen = 1, xtickname = xx,$
	xtitle = ' ', ymargin = [0,2], yrange = [ymin,ymax], $
	timerange=trange
outplot, tarray[s:e],yclean[s:e,1], $
	color = 2

; Plot emission measure

;Special correction for flare on 28-oct-2003
;Correct emission measures where channel 1 drops out.
; Make them equal to the previous value.
;emis(where(yclean lt 1.e-7)) = emis(where(yclean lt 1.e-7) - 1)
;emis(where(yclean lt 1.e-7)) = emis(where(yclean lt 1.e-7) - 1)

ymax = max(smooth(emis(s:e),10,/NAN),min = ymin)
utplot,tarray(s:e),emis(s:e),utbase, $
	/ylog,ytitle='Emission Measure in 10!u49!n cm!u-3!n', $
	color=0,background=255, xticklen = 1, xtickname = xx,$
	xtitle = ' ', ymargin = [0,0], yrange = [0.1,ymax], $
	timerange=trange

c_size = 1.
xyouts,0.4,0.95,'Start Time = ' + stime,color=0,/normal,charsize = c_size
xyouts,0.4,0.93,'End Time = ' + etime,color=0,/normal,charsize = c_size

; Plot temperature in MK

;Special correction for flare on 28-oct-2003
;Correct temperatures where channel 1 drops out
; Make them equal to the previous value.
;tempr(where(yclean lt 1.e-7)) = tempr(where(yclean lt 1.e-7) - 1)
;tempr(where(yclean lt 1.e-7)) = tempr(where(yclean lt 1.e-7) - 1)

ymax = max(smooth(tempr(s:e),10,/NAN),min = ymin)
if ymax gt 50. then ymax = 50.0
if ymin lt 0.0 then ymin = 0.0
utplot,tarray(s:e),tempr(s:e),ytitle='Temperature in MK', $
	color=0,background=255,xticklen = 1, xtickname = xx,$
	xtitle = ' ', ymargin = [0,0],yrange = [5.0,ymax+1.] , $
	timerange=trange

;Special correction for flare on 28-oct-2003
; Correct temperatures where channel 2 is saturated from interval 2515 to 2703
;slope = (26.8788 - 19.9531)/188.0
;for i = 0, 188 do tempr(i+2515) = 26.8788 - i * slope

outplot,tarray(s:e), tempr(s:e), color = 2

;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 * tempr * sqrt(emis * V)
	;Time of maximum thermal energy in plasma, Umax
	Umax = max(U[fsti:feti])
	;Array index of Umax
	iUmax = where(U eq Umax)
	iUmax = iUmax[0]

ymax = max(smooth(U(s:e),10,/NAN),min = ymin)
utplot,tarray(s:e),U(s:e),/ylog,ytitle='Thermal Energy U in ergs', $
	color=0,background=255,xticklen = 1, xtickname = xx,$
	xtitle = ' ', ymargin = [0,0], yrange = [ymin,ymax], $
	timerange=trange

xyouts,0.5,0.3,'Volume = ' + num2str(V,format='(e10.1)') + ' cm!u3!n', $
	color=0,/normal,charsize = c_size

; Skip exponential decay plots
GOTO, Label1
;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.2,0.55,charsize=1.0, $
	'Zero Intercepts = '+ num2str(A(*)), /normal,color=0, charsize = c_size

xyouts,0.2,0.52,charsize=1.0, $
	'Time Constants (s) = '+ num2str(t(*)), /normal,color=0, charsize = c_size

Label1:

; 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)

ymax = max(smooth(Udot(s:e),10,/NAN),min = ymin)
utplot,tarray(s:e),Udot(s:e),ytitle='dU/dt in ergs s!u-1!n', $
	color=0,background=255, $
	xticklen = 1, $ ;xtitle = ' ',  xtickname = xx,$
	ymargin = [4,0],yrange=[-5.e27,ymax], $
	timerange=trange

!p.multi = [0,1,3]
if ps eq 0 then window,11, title = 'Thermal Energies vs. Time', xsize = x_pixels, ysize = y_pixels

; Calculate radiative loss rates
	Lrad = loss_rad_chianti(emis,tempr*1.e6)

; Calcualate conductive loss rate

; Calculate loop cross-sectional area in cm^2 & loop length in cm
; Loop cross-sectional area in arcseconds x arcseconds
; = ribbon width x length
arcsec2cm = 7.253e7 ; cm per arcseconds
ribbon_width = 10 ;arcseconds
ribbon_length = 100 ;arcseconds
area_sqarcsec =  ribbon_width * ribbon_length ;area in square arcseconds
area = area_sqarcsec * (arcsec2cm)^2 ;area in square cm

loop_length = 50 ;arcseconds
length = loop_length * arcsec2cm

	Lcond = loss_cond(area,length,tempr*1.e6)

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

ymax = max([smooth(Pin(s:e),10,/NAN), smooth(Lcond(s:e), 10,/NAN)], min = ymin)
ymin = 1.e25
utplot,tarray(s:e),Pin(s:e),color = 0, $
	ytitle='Energy loss rate in ergs s!u-1!n', $
	/ylog,yrange = [ymin, ymax], $
	background=255,xticklen = 1,$
	xtitle = ' ',  xtickname = xx, ymargin = [0,2], $
	timerange=trange

outplot,tarray(s:e),lrad(s:e),color=4
outplot,tarray(s:e),lcond(s:e),color=2

outplot,tarray(s:e),Udot(s:e),color=10
outplot,tarray(s:e),-Udot(s:e),color=8

outplot, tarray[s:e], flux[s:e], color = 13

;Get SORCE data to plot
;	path = 'C:\My_Documents\My_HESSI\2003\oct-nov\10-28\SORCE\'
;t = rd_ascii(path + 'times.txt') ;t is in hh:mm:ss for 28 October 2003
;ts = anytim(t) ;Puts t in seconds of day
;dt = ts(1) - ts(0) ;delta-t in seconds
;n = rd_ascii(path + 't_numbers.txt') ;n is the flux in watts/m^2
;n = float(n)

;change from watts/m^2 to J/s
;D = 1.495e11 ; 1AU in cm
;n = n * 2. * !pi * D^2 ;n is now in Joules / s
;change from J/s to rgs /s
;n = n * 1.e7 ; n is now in erg/s
;n = n - 1.9070e33 ;n is now the increase in flux from the flare in erg/s

;Plot SORCE data
;outplot,ts,n,color = 9

;bcl = 11
;nsmooth = smooth(n,bcl)

;utbase = 7.8330240e8

xpos=0.4 & ypos = 0.96 & dypos = 0.02
xyouts, xpos,ypos, 'Total Eneregy Loss Rate - black', color=0, /normal, charsize = c_size
;xyouts, xpos, ypos, 'SORCE Increase in TSI - light blue', color = 9,  /normal, charsize = c_size
xyouts, xpos,ypos - dypos, 'Radiative Loss Rate - brown', color=4, /normal, charsize = c_size
xyouts, xpos,ypos - 2.*dypos, 'Conductive Loss Rate - red', color=2, /normal, charsize = c_size
xyouts, xpos+0.05,ypos - 3.*dypos, $
	'Area = ' + num2str(area,format='(e10.1)') + ' cm!u2!n  ' + $
	'Length = ' + num2str(length,format='(e10.1)') + ' cm', $
	color=2, /normal, charsize = c_size
;xyouts,  xpos+0.05,ypos - 4.*dypos, $
;	'Length = ' + num2str(length,format='(e10.1)') + ' cm', $
;	color=2, /normal, charsize = c_size
xyouts, xpos,ypos - 4.*dypos, 'Udot - blue', color=10, /normal, charsize = c_size
xyouts, xpos+0.13,ypos - 4.*dypos, '-Udot - green', color=8, /normal, charsize = c_size
xyouts, xpos,ypos - 5.*dypos, 'GOES X-ray flux', color=13, /normal, charsize = c_size

; Integrate lrad, lcond, and Pin wrt time to get the total energy input to the thermal plasma

; Zero parameters before and after flare
lrad[0:fsti] = 0.0 & lrad[feti:n_tarray] = 0.0
lcond[0:fsti] = 0.0 & lcond[feti:n_tarray] = 0.0
Pin[0:fsti] = 0.0 & Pin[feti:n_tarray] = 0.0
Udot[0:fsti] = 0.0 & Udot[feti:n_tarray] = 0.0
flux[0:fsti,*] = 0.0 & flux[feti:n_tarray,*] = 0.0

;Add fluxes from the two GOES channels together and subtract pre-flare background fluxes.

flux0= flux[*,0] - flux[backi,0]
flux1= flux[*,1] - flux[backi,1]
flux = flux[*,0] + flux[*,1] - flux[backi,0] - flux[backi,1]

;Determine length of each time interval
;n_t = n_elements(tarray)
delta_t = tarray[1:n_tarray] - tarray[0:n_tarray-1]
delta_t = [delta_t, delta_t(n_tarray-2)]

i_lrad = integrate(lrad) * delta_t
i_lcond = integrate(lcond) * delta_t
i_Pin = integrate(Pin) * delta_t
i_Udot = integrate(Udot) * delta_t + U(0)
i_flux = integrate(flux) * delta_t
i_flux0 = integrate(flux0) * delta_t
i_flux1 = integrate(flux1) * delta_t
;i_n = integrate(n(32:52)) * dt ;SORCE integrated radiated energy in ergs

ymax = max([i_lrad,i_lcond,i_Pin],min = ymin,/NAN)
ymin = 1.e28

utplot,tarray[s:e],i_Pin[s:e],ytitle='Energy Input to Thermal Plasma in ergs', $
	color=0,background=255, ymargin = [1,1],  yrange = [ymin,ymax], /ylog, $
	xtitle = ' ',  xtickname = xx,  $
	xticklen = 1, timerange=trange
outplot,tarray[s:e],i_lrad[s:e], color=4
outplot,tarray[s:e],i_lcond[s:e], color=2
;outplot,tarray[s:e],i_Udot[s:e], color = 12
outplot,tarray[s:e], U[s:e] - U[fsti], color = 10
outplot,tarray[fsti:feti], i_flux[fsti:feti], color = 13
;outplot,ts[32:52],i_n, color = 9

;Total Energy Input to Thermal Plasma
Etotal =  I_Pin(feti)
Rtotal = i_lrad(feti)
xpos = 0.4 & ypos = 0.62 & dypos = -0.02
xyouts,xpos,ypos, $
	'Peak energy of plasma = ' + num2str(Umax,format='(e10.1)') + ' ergs', $
	color = 0, /normal, charsize = c_size
xyouts,xpos,ypos + dypos, $
	'Total energy radiated from plasma = ' + num2str(Rtotal,format='(e10.1)') + ' ergs', $
	color = 0, /normal, charsize = c_size
xyouts,xpos, ypos + 2. * dypos, ' up to ' + flare_end_time, color = 0, /normal, charsize = c_size
;xyouts,xpos, Ypos + 3. * dypos, 'Total SORCE radiated energy = ' + num2str(i_n(20),format='(e10.1)' ) + ' ergs', $
;color = 0, /normal, charsize = c_size

;Plot L_total / L_x ratio
;if ps eq 0 then window,12, title = 'Ratios vs. Time', xsize = x_pixels, ysize = y_pixels
utplot,tarray[fsti+1:feti-1],lrad[fsti+1:feti-1]/flux[fsti+1:feti-1],ytitle='L!dtotal!n / L!dx!n ', $
	color=0,background=255, ymargin = [4,0], yrange = [0,200],  $ ;/ylog, $
	xticklen = 1, timerange=trange ;, xtitle = ' ',  xtickname = xx
outplot,tarray[fsti+1:feti-1],i_lrad[fsti+1:feti-1]/i_flux[fsti+1:feti-1], color=2
	;ytitle='Ratio of Time Integrated Fluxes - IL!dtotal!n / IL!dx!n ', $
	 ;ymargin = [4,0],  $ ;background=255, yrange = [ymin,ymax], /ylog, $
	;xticklen = 1, timerange=trange
xyouts, 0.25, 0.30, 'Radiated flux / GOES flux', color = 0, /normal, charsize = c_size
xyouts, 0.25, 0.28, 'Time integrated radiated flux / time integrated GOES flux', $
	color = 2, /normal, charsize = c_size
xyouts, 0.25, 0.26, 'L!dtotal!n / L!dx!n up to ' + flare_end_time + ' = ' + $
	num2str(i_lrad(feti-1)/i_flux(feti-1), format = '(f6.1)'), $
	color = 2, /normal, charsize = c_size

;if ps eq 0 then window,12, title = 'Thermal Energies vs. Time', xsize = x_pixels, ysize = y_pixels

;Print parameters at specific times.

close,/all
print, 'Open ', out_file
openw, 1, out_file

printf,1, flare_start_time, ' to ', flare_end_time
printf,1, V, format = '("Assumed source volume = ", E10.2, " cm^3")'
for i = iUmax, feti, feti - iUmax DO begin
	printf,1, '	  '
	if i eq iUmax then begin
		printf, 1, 'Energies up to time of maximum thermal energy (t_Umax) = ', anytim(utbase + tarray(iUmax),/vms), ' UT'
		printf, 1,  Umax, alog10(Umax), format = '("Maximum thermal energy (Umax) = ", E10.2, "  ergs, log10 = ", F5.2)'
	endif
	if i eq feti then begin
		printf, 1, 'Energies from start of flare at ', flare_start_time, ' UT to end of flare at ', flare_end_time, ' UT'
		printf, 1, Etotal, alog10(Etotal), format = '("Total Energy Input to Thermal Plasma = ", E10.2, " ergs, log10 = ",  F5.2)'
	endif
	;Print radiated and conducted energies
	printf, 1,  i_lrad(i), alog10( i_lrad(i)), format = '("Radiated Energy = ", E10.2, "  ergs, log10 = ", F5.2)'
	printf, 1,  i_lcond(i), alog10(i_lcond(i)), format = '("Conducted Energy = ", E10.2, "  ergs, log10 = ", F5.2)'
	printf, 1, i_flux(i), alog10(i_flux(i)), format = '("Lx = ", E10.2, " ergs, log10 = ", F5.2)'
	printf, 1, i_flux0(i), alog10(i_flux0(i)), format = '("Lx low = ", E10.2, " ergs, log10 = ", F5.2)'
	printf, 1, i_flux1(i), alog10(i_flux1(i)), format = '("Lx high = ", E10.2, " ergs, log10 = ", F5.2)'
	printf, 1, i_lrad(i)/i_flux(i), format = '("Ratio L_total/L_x = ", F5.2)'
endfor

	Total_energy =  Umax + i_lrad(iUmax) + i_lcond(iUmax)
	printf, 1,  Total_energy, alog10(Total_energy), format = '("Total Energy input to thermal plasma = ", E10.2, "  ergs, log10 = ", F5.2)'
	printf, 1,  'from ', flare_start_time, ' through ', flare_end_time

;Close all open files
print, 'Close ', out_file
close,/all

;Close device for postscript plots
if ps gt 0 then device,/close

;Restore to original state for future plots
	restore_plotvar
	set_plot,old_device

end