pro hessi_figure4_inputs, e_si, si_counts, therm_fluence, edg2b, ph, cfront, crear, ctotal


;We just want the thermal component!
;read in the fluxes used in HESSI_SPEC.PRO, only to extract the thermal component
;CREATE THE INCIDENT SOLAR FLARE PHOTON FLUXES
;-------------------------------------------
;Read in Fluence file
;** Structure HESSI_FLUX, 7 tags, length=67840:
;   EM              FLOAT     Array(3335)
;   THERM1          FLOAT     Array(3335)
;   BPOW1           FLOAT     Array(3335)
;   NUCL1           FLOAT     Array(3335)
;   PI_NOT          FLOAT     Array(140, 2)
;   COMP            FLOAT     Array(3335)
;   README          STRING    'comp is fluence divided by 100 sec'
; THERM1 IS fluence of 1e50 cm-3 at 20 MegaKelvin over 100 seconds
read_hessi_fluxes, hflux
etherm = [findgen(100)*.1+1.,11+findgen(50)] 
flux_therm = interpol(hflux.therm1,hflux.em,etherm) * 4;fluence in 400 sec (photons/keV)
;add the power-law to create the si_flux
;then later rebin the thermal to the bins of the photons bins for the Ge's
si_fluence = flux_therm+ 4.*5685.*(30./etherm)^2.5 ;photons/keV/cm2
e_si = etherm

;Response was created in two separate Monte Carlo runs
;by Said Slassi, one at high energy starting at 300 keV
;at a fluence 4x the 27 April 1981 gamma-ray event
;Our model for comparison will be this fluence over
;400 seconds to obtain the fluxes
;
;The continuum fluences are:
;10keV<E<30keV, 5685(E/30keV)**-2.5 ph/keV.cm2
;30keV<E<76.17keV, 890(E/50keV)**-3.63 ph/keV.cm2
;E>76.17 keV, 0.4(E/1000keV)**-2.4 ph/keV.cm2
; These are over 100 seconds.  All of this is represented in
; his data files "*.spc", 5 of which we will use in this figure:
;        counts27Aprx4.spc
;        newfrontsonly.spc
;        newphotoninput.spc
;        newrearsonly.spc
;        photons27Aprx4.spc
; The file names are self-explanatory in light of the
; previous discussion.  Also, for this model the statistics
; are calculated for the source centered over the collimator
; axis resulting in an average transmission of 36% for perfect
; grids instead of the expected average of 25%.  This changes the
; important statistics at high energy only marginally.
; merge photon inputs files at 4x 27 April level
; all of these are fluences integrated over time and total detector area!
hessi_read_spc,'photons27Aprx4.spc',ema,edg2a,ph4
hessi_read_spc,'newphotoninput.spc',emb,edg2b,ph
;energies ema and emb are in MeV
ph=ph*4
nstart = n_elements(emb)-n_elements(ema)
ph(0,nstart)=ph4 
ph_therm = interpol(flux_therm, etherm/1000., emb)* 12*!pi*3.55^2 ;photons/keV/det_area
ph = ph + ph_therm

;-------------------------------------------
;calculate the counts in the GeDs from this thermal flux
;Grids
trans = 0.6^2 

;GE WINDOW ATTENUATION
read_seqfile,buff,loc_file( path=hessi_data_paths(),'hessi_windows.dat')
ge_atten=fltarr(100,2)
iline=(where(strpos(buff,'100 attenuation') ne -1))(0)+1
reads, buff(iline:*), ge_atten
ge_atten=transpose(ge_atten)
ge_fil = interpol( (ge_atten(1,*))(*), (ge_atten(0,*))(*)/1000., emb)  $
    < max(ge_atten(1,*))
;-------------------------------------------
trans = ge_fil * trans 



;merge response vectors too!
hessi_read_spc,'counts27Aprx4.spc',ema,edg2a,ctotal4
hessi_read_spc,'newrearsonly.spc',emb,edg2b,crear 
hessi_read_spc,'newfrontsonly.spc',emb,edg2b,cfront

crear = crear *4
cfront = cfront *4 + trans*ph_therm
ctotal = crear + cfront
ctotal(0,nstart) = ctotal4
crear = ctotal - cfront

;SI DETECTORS
;-------------------------------------------
si_scale = 0.6^2*12.
si_eff = 1-exp(-.2*xsec((e_si<999.)>1.01,14,'pe'))

;Extract the attenuation coefficients for the Si detectors
iline=(where(strpos(buff,'Energies') ne -1))(0)+1
si_energ=fltarr(140)  
reads, buff(iline:*), si_energ
iline=(where(strpos(buff,'Effective area cm2') ne -1))(0)+2
si_eff_fil  =fltarr(140)  
reads, buff(iline:*), si_eff_fil
si_fil_energ =si_eff_fil/( 1-exp(-.2*xsec(si_energ,14,'pe')))
;take sharp edges out of filters here,  filter should really be applied inside
;comp_spec_calc.pro, but for the purposes of the proposal figure we can do it here.
em_smooth = findgen(200)*.25 +1.
si_fil = smooth( interpol( si_fil_energ, si_energ, em_smooth),8)
si_fil = interpol(si_fil, em_smooth, e_si)

w50 = where(e_si ge 50.,n50)
if n50 ge 1 then si_eff(w50)=0.0
wlld_si = where( e_si lt 2.,nlld)
if nlld ge 1 then si_eff(wlld_si) = 0.0

si_eff_area = si_eff*si_scale*si_fil
;si_fluence has units of photons/cm2/keV
si_counts = si_fluence * si_eff_area
e_si = e_si / 1000. ;convert output energies to MeV, fluences are /keV

therm_fluence = si_fluence * 12*!pi*3.55^2 ;get the low energy flux in Ge units
end