PRO xspec_write_pha, event_es, arf_file, pha_example_file, output_pha_file
;
; Fill an XSPEC pha file with an event histogram
; Hardly optimized but it works!
;
; INPUTS:
;    event_es            array of event energies in keV
;    arf_file            arf file to read in to get energy bins
;                           assuming the arf file is square
;    pha_example_file    a template pha file to provide header
;                           information - e.g., created by XSPEC fakeit
;    output_pha_file     file name for the output file
;
; 6/4/97 dd

; Get the energy bin limits from arf file
arf = mrdfits(arf_file,1)
e_lo = arf.energ_lo
e_hi = arf.energ_hi

; Create the histogram array to fill
nhist = n_elements(e_lo)
hist = FLTARR(nhist)

; Loop through the events and fill the histogram
; The e_lo and e_hi can have gaps between them so I
; just use the e_hi and count the event in the bin
; just below the first one above its energy.
;
print, ' xspec_pha: binning '+STRING(n_elements(event_es))+' events ...'
ie = LONG(0)
while ie LT n_elements(event_es) do begin
  ; which bin?
  hi_bins = where(e_hi GT event_es(ie), nfound)
  if nfound GE 1 then begin
    bin = hi_bins(0)-1 > 0
    hist(bin) = hist(bin) + 1.0
  end else begin
    print, '* No bin selected for event index '+STRING(ie,FORMAT='(I10)')+'
    print, '   E_event = '+ STRING(event_es(ie))
  end
  ie = ie + 1
end

; Read in the (faked) pha file header and extension header
; and write it out
pri_header = headfits(pha_example_file,EXTEN=0)
ext_header = headfits(pha_example_file,EXTEN=1)

; Header values could be changed here if desired

; Create the output FITS file
fxwrite, output_pha_file, pri_header

; Create the binary table extension
; for the histogram data...
fxbcreate, unit, output_pha_file, ext_header, errmsg=errmsg

; And write the data...
; 
; write the channels
col = 1  + 0*indgen(nhist)
row = 1 + indgen(nhist)
channel = row
for id=0,nhist-1 do begin
  fxbwrite, unit, channel(id), col(id), row(id)
end
; write the counts
; The pha channels are in dispersion distance and so
; are flipped in order, so flip them when outputing hist():
col = 2  + 0*indgen(nhist)
row = 1 + indgen(nhist)
for id=0,nhist-1 do begin
  fxbwrite, unit, hist((nhist-1)-id), col(id), row(id)
end

; close it up
fxbfinish, unit

RETURN
END