; flare_energy_x13_x27.pro
; 
; Used to make Figure 1 of the Xray-IE paper
;
; header = 'C:\Users\dsiskind\Documents\science\LWS_flares\'
file1='20031028_halfA'
file2='20170910_7nm'
suffix = '.sav'

;
  flarefile1= file1 + suffix
  flarefile2 = file2 + suffix
  print,flarefile1,'  ', flarefile2
restore,flarefile1
wave1=wavelength(1:139)     ;through out the 0.25 value
flux1=flare_spectrum
imax1=669
peakflux1=flux1(imax1,1:139)
restore,flarefile2
wave2=wavelength
flux2=flare_spectrum
imax2=966
peakflux2=flux2(imax2,0:138)
;
energy_ev=12397/wavelength    ;photon energy in ev
energy_erg = 1.6e-12*energy_ev   ;photon energy in erg
energy_joule = 1e-7*energy_erg    ;photon energy in Joule
;
maxflux1 = 1e4*peakflux1              ;convert flux to photons/s/m2
maxflux2 = 1e4*peakflux2
;   
  flare_joule1 = maxflux1*energy_joule      ;yields (J/s/m2 or W/m2)
  flare_joule2 = maxflux2*energy_joule
powerbin1=fltarr(7)                   ;units of power: W/m2
powerbin2=fltarr(7)

powerbin1(0)=total(flare_joule1(0:18))
powerbin2(0)=total(flare_joule2(0:18))

xray1 = total(flare_joule1(0:22))
xray2 = total(flare_joule2(0:22))

print,'integrated energies > 1 keV after 1 year (in kJ)'
print,3.1e7*xray1/1000,3.1e7*xray2/1000,3.1e7*xray3/1000
for ibin=0,5 do begin
  print,ibin,wave1(19+ibin*20),wave1(38+ibin*20)
  powerbin1(ibin+1) = total(flare_joule1(19 + ibin*20:38+ibin*20))
  powerbin2(ibin+1) = total(flare_joule2(19+ibin*20:38+ibin*20))
 
endfor
wavebin=findgen(7)+0.5
!x.title='Wavelength (nm)'
!y.title='W/m!u2!n'
plot_io,wavebin,powerbin1,psym=10,yr=[1e-6,1e-2],xr=[0,7],xs=1,thick=2
oplot,wavebin,powerbin2,linestyle=2,psym=10,thick=2
oplot,wave1/10.,flare_joule1
oplot,wave2/10.,flare_joule2,linestyle=2
;
end