pro fits_spectra, filename, ut, rate, erate=erate, livetime=livetime, $
    start_time=start_time, end_time=end_time, edges=edges, width=width, $
        det_id=det_id, dsel=dsel, addpar=addpar, write=write, printx=printx, $
        channels=channels, headers=headers, data_unit = data_unit, error=error

;on_error, 2                    ;return to caller
                  
;Distinguish between OGIP and BFITS keywords and standards

phead = headfits(filename)
instrument = fxpar(phead, 'INSTRUME')
mnemonic   = fxpar(phead,'MNEMONIC')

;We are either reading a BFITS or OGIP standard binary table
if (where_are(instrument,'BATSE'))(0) eq 0 and (where_are(mnemonic,'BFITS'))(0) eq 0 then $
    bfits = 1 else bfits = 0

;EXTENSION NAMES for energy calibration and time series data
ebounds = 'EBOUNDS' ;name for pha energy calibration
spectrum= 'SPECTRUM'    ;count (rate) data
edge_labels   = ['E_MIN','E_MAX']
if bfits then begin
   ebounds = 'BATSE_E_CALIB'    
   spectrum= 'BATSE BURST SPECTRA'
   edge_labels  = 'E_EDGES'
endif

error = 1
if keyword_set(printx) then printproc = 'printx' else printproc = 'print'

file = findfile(fcheck(filename,''),count = nfile)
if nfile ne 1 then begin
    call_procedure,printproc, "'"+filename+"' does not exist or is not unique, Error."
        return
endif

headers = phead

for extens = 1,2 do begin
  fxbopen, unit, filename, extens, header
  headers = temporary([headers,header])



  fxbfind, header, 'ttype',type_columns,types,n_found
  fxbfind, header, 'tunit',unit_columns,units,n_found
  types = strtrim(types,2)
  units = strtrim(units,2)


  case strtrim(fxpar( header, 'extname' ),2)  of 
    ebounds: begin  ;PHA energy edges
       if bfits then begin
          fxbread, unit, edges, edge_labels
      if (size(edges))(0) eq 2 then begin
        naxis2= (size(edges))(2)
        edges_2 = fltarr(2, (size(edges))(1) -1,  naxis2)
            for i=0,naxis2-1 do begin
            edge_products, edges(*,i), edges_2 = edg
            edges_2(0,0,i) = edg
        endfor
        edges=edges_2
      endif else edge_products, edges, edges_2=edges
      
       endif else begin
          fxbread, unit, channels, 'channel'
          fxbread, unit, e_min, edge_labels(0)
          fxbread, unit, e_max, edge_labels(1)
      edges = transpose( [[e_min],[e_max]])
       endelse
    end
    spectrum: begin     ;spectral data vs time

       if bfits then begin
         basetime=fxpar(header,'basetime',comment=comment)
         wtime=(where( strpos(types,'TIME') ne -1))(0)
         fxbread, unit, ut ,types(wtime)      
;We have assumed basetime is in TJD, times in seconds for the moment
         ut = ut + tjd2ymd( basetime)
         lo_chan  = fxpar(header, 'LO_CHAN')
         up_chan  = fxpar(header, 'UP_CHAN')
         channels = lo_chan + indgen(up_chan-lo_chan + 1)
       endif else begin 

;Here we follow one of the OGIP standards, the SPEX standard output format for binary tables
         fxbread, unit, ut, 'TIME'
         fxbread, unit, ut_del, 'TIMEDEL'
         ut = transpose( [[ut-ut_del/2],[ut+ut_del/2]])
     timeunit = (strlowcase(units(where(types eq 'TIME'))))(0)
         case timeunit of
       's':  ;no problem and no changes
           'd': ut = ut * 86400.0d0 ;change to seconds
          else: message,'Cannot resolve time system!'
         endcase
         tstart = fxpar( header, 'TSTART')
         tstop = fxpar( header, 'TSTOP')
         timezero = fxpar( header, 'TIMEZERO')
         timesys = strtrim(fxpar( header, 'TIMESYS'),2)
         timeunit = strtrim(fxpar( header, 'TIMEUNIT'),2)
         case timeunit of
       's': ut  =  ut +tjd2ymd(timezero) 
           'd': ut =   ut +tjd2ymd(timezero* 86400.0d0) ;change to seconds
          else: message,'Cannot resolve time system!'
         endcase
       endelse          ;closing OGIP time standard

;Deadtime corrections
         if fxpar( header, 'DEADAPP') or fxpar( header, 'LTIMECOR') then begin
           dead_index = where_are(types,'deadc', dcount)

         ;Here we use livetime as the variable name, but it means deadtime fraction until converted!!

           if dcount eq 1 then fxbread, unit, livetime, type_columns(dead_index(0)) else begin
             livetime = fxpar(header,'DEADC*' ,count=dcount)
             if dcount gt 1 then livetime = $
              rebin( reform(livetime,n_elements(livetime),1),n_elements(livetime),n_elements(ut(0,*)))
           endelse
           livetime = 1.-livetime       ;now convert deadtime to real livetime fraction
         endif

       checkvar, start_time, anytim(min(ut),/yohkoh)
       checkvar, end_time, anytim(max(ut),/yohkoh)
       rows = where( ut(0,*) ge anytim(start_time,/sec) and  $
                      ut(0,*) lt anytim(end_time,/sec) , nrows )
       if nrows eq 0 then begin
          call_procedure,printproc, ['No data found between start and end times', $
          anytim(start_time,/yohkoh)+' - '+anytim(end_time,/yohkoh), $
         'Start and Stop Time of data arrays: '+ atime(/yohkoh,limits(ut))]
         rate =-1
         erate=-1
       endif else begin 

;Yes, there should be some spectral data!
        data_list = ['count','rate','flux','spectrum']
    i = -1
        repeat begin
          i = i + 1
          wcount = where_are(types, data_list(i), exclude = 'err', ncount)
        endrep until ncount gt 0
        data_col = type_columns(wcount(0))
;Read the data between START and END times
        fxbread, unit, rate, data_col, 1+limits(rows)
        data_unit = units((where(unit_columns eq data_col))(0))

;Read the errors if they are present
        werr     = where_are(types, 'err', nerr )
        if nerr gt 0 then fxbread, unit, erate, type_columns(werr(0)), 1+limits(rows) $
           else erate = -1
        error = 0 
       endelse
   end
    else:   begin
    stop
    message,'Binary Table Format is unexpected here. Returning.'
    end
  endcase 


  fxbclose,unit

endfor


end