;+
;
; NAME:
;	READ_HXRS_FITS
;
; PURPOSE:
;       Reads hxrs fits files into idl, then
;       computes counts/sec and live times for shielded and non-shielded
;       detectors and effective total livetime for sum of both s. and ns.
;
; CATEGORY:
;       HXRS
;
; CALLING SEQUENCE:
;
; CALLS:
;       CHECKVAR, EDGE_HXRS, F_DIV, ATIME, UTPLOT, PRINTX
;
; INPUTS:
;	Through keywords so you don't have to remember the input order.
;
; OUTPUTS:
;	Through keywords so you don't have to remember the output order.
;
; KEYWORDS:
;        file:	hxrs fits file name to read (input)
;        ut:	time in seconds since 1-jan-1979, middle of datapoint,
;		dblarr(2,n) (start and end of each interval) (output)
;        syarray, nyarray, totyarray: shielded, nonsh and total raw counts/interval
;        esyarray, enyarray,etotyarray: uncertainties, i.e. sqrt(yarray(10,N)
;        deltat:	accumulation time interval, scalar, in seconds (output)
;        slivetime, nlivetime, totlivetime: livetimes per interval for detectors and their total
;        edges:	output energy loss bands in keV  fltarr(2,10) (output)
;        h0:	primary header (string array) (output)
;        h1:	header from first (time) extension (string array) (output)
;        error:	if =1, then error reading file.  Otherwise 0. (output)
;
; PROCEDURE:
;	User may explicitly pass the name of a HXRS FITS file;
; 	Sample call:
;	read_hxrs_fits, file='hxrs_fits:hxrs_23.fits', ut=ut, syarray=syarray
;
; MODIFICATION HISTORY:
;        Adapted from the get_hxrbs_fits routine
;	Written by Kim Tolbert (NASA/GSFC), 9-sep-1994.
;        Adapted for HXRS by Alan Kiplinger May 2000.
;   5-June-2000 Paul Bilodeau, general re-organization, error handling
;-


PRO READ_HXRS_FITS, file=file, ut=ut, syarray=syarray, nyarray=nyarray, $
         totyarray=totyarray, esyarray=esyarray, enyarray=enyarray,     $
         etotyarray=etotyarray, deltat=deltat, slivetime=slivetime,     $
         nlivetime=nlivetime, totlivetime=totlivetime, edges=edges,     $
         h0=h0, h1=h1, error=error

; Assume an error
error = 1

; handle errors via printed messages and return
catch, error_status
if error_status ne 0 then begin
	printx,'Error in read_hxrs_fits:'
	printx, !err_string
	error=1
	return
endif

; if no file name was given, prompt the user for one
if (size(file))[1] eq 0 then begin
   file=' '
   read,' Enter input filename: ',file
endif

; read first extension and its header
datain=mrdfits(file,1,h1)

; relative seconds since first interval
times = double(datain[*].time - datain[0].time)

; millisecond corrections to start times
finetimes = double(datain[*].finetime)

; real startime of 1st interval
date_obs = STRTRIM(fxpar(h1, 'DATE_OBS'), 2)

; get time of first interval since Jan. 1, 1979
t0 = anytim(date_obs)

 ; start times of each interval corrected for t0 and finetimes
times = times + finetimes/1000.d0 + t0

; calculate deltat, setting the final interval length to either 1s or 0.2s
deltat = double(shift(times,-1) - times)
if datain[n_elements(deltat)-1L].blocktype eq 8 then $
	deltat[n_elements(deltat)-1L]=0.2 $
else deltat[n_elements(deltat)-1L]=1.

; use data with blocktype values of 0 or 8 and a positive time interval length
use=where( (datain.blocktype eq 0b or datain.blocktype eq 8b) and deltat gt 0.d0, nbin)
if nbin gt 0L then begin
	counts_s = datain[use].counts_s
	counts_n = datain[use].counts_n
	times = times[use]
	deltat = deltat[use]
	blocktype = datain[use].blocktype
	endif else begin
	printx,'Error in read_hxrs_fits:'
	printx,'No flare mode or normal mode data.'
	error=1
	return
	endelse

; catch deltat values that are too long
flare = where(blocktype eq 8, n_flare)
if n_flare gt 0L then deltat[flare]=deltat[flare] < 0.2
normal = where(blocktype eq 0, n_normal)
if n_normal gt 0L then deltat[normal]=deltat[normal] < 1.

; ut array
ut = dblarr(2, nbin)
ut[0,*] = times
ut[1,*] = times + deltat

; shielded data
syarray = float(counts_s + 32768)
esyarray = sqrt(syarray)
slivetime = deltat - ( total(syarray(0:7,*),1)*4e-6 +  syarray(8,*)*4e-5 )

; non-shielded data
nyarray=float(counts_n + 32768)
enyarray=sqrt(nyarray)
nlivetime = deltat - ( total(nyarray(0:7,*),1)*4e-6 +  nyarray(8,*)*4e-5 )

; total data
totyarray = syarray + nyarray
etotyarray = sqrt(totyarray)
totlivetime = make_array(nbin, /float, value=1.)
q = where( (slivetime gt 0.) and (nlivetime gt 0.) and (totyarray gt 0.), nq )
if nq gt 0 then totlivetime[q] = totyarray[q] / ( (syarray[q]/slivetime[q]) + (nyarray[q]/nlivetime[q]) )

; compute energy channel widths
edges_hxrs, edges

; hxrs central detector area
area = 4.909;cm^2

error=0

return

end