;+
; PROJECT:
;	SDAC
; NAME:
;	READ_DISCLA
; PURPOSE:
; 	Procedure to read BATSE FDB files and return time, DISCLA, and angle
; 	information.  The DISCLA data has been corrected for overflow and 
;	divided by livetime.
;
; 	You must supply the file name explicitly, the time you want, or the 
;	flare number in order for read_discla to find the correct input file.  
;	If you specify a file name or flare number, you may also specify 
;	time start and end to restrict the amount of data returned.  If 
;	you don't specify any times, fdbread will return all data in the 
; 	fdb, bdb, or the msfc discla file. 
; 	UTPLOT package times are set as follows:  base time is set to start 
;	of day of the current data, start and end times are set to 0 for 
;	autoscaling.
;	This procedure is based on FDBREAD but the data values aren't placed
;	into the common blocks, but read directly.
; RESTRICTIONS:
;	Doesn't fully support start and end times of all MSFC FITS formats.
; CALLING SEQUENCE:
;	read_discla, seconds, discla, cosines	
; CATEGORY:
;	BATSE
; OUTPUTS:
;   seconds - (output)   array of times at middle of accumulation interval in 
;                  seconds since start of day.  
;                  Base time for UTPLOT routines will be set to 0h 0m
;                  on day of data. (Function GETUTBASE(0) will return base 
;                  time in seconds since 79/1/1.) 
;   discla  - (output)   (i,j,n) floating point array of discla data where
;                  i=6 - channels 1-4 (energy bins 25-50, 50-100, 100-300, and
;                        >300 keV), $
;			 NaI-(channels 1-4 +CPD, effectively 15-25 keV), and CPD
;                  j=8 - LADs 1 through 8
;                  n - number of data intervals 
;   cosines - (output) floating point array of cosine of angles between 8 LADs
;                  and Sun (cosines(j) corresponds to y(*,j,*)
;
; OPTIONAL KEYWORD INPUTS:
;   filename -  string FDB, BDB, or DISCLA file name
;   starttime - start time of data to accumulate in r*8 seconds 
;               since 79/1/1,0 or as ASCII string 'YY/MM//DD,HHMM:SS'.
;   endtime -   end time of data to accumulate. See starttime above for format.
;   flare -     BATSE flare number to accumulate data for
;   verbose -   0/1 means don't/do print accumulation times message
;
; OPTIONAL KEYWORD OUTPUTS:
;   livetime -  (8,n) floating point array of live time in seconds for 
;               each LAD for each data interval 
;   error -     (keyword output) =0/1 means no error / error
;   cos_spec -  (output) floating point array of cosine of angles between 8 SPECs
;                  and Sun (cos_spec(j) corresponds to y(*,j,*)
;		The cosines are selected from the start of the time interval
;		which could be in error if there is a pointing change
;		which comes about once every two weeks.
;   USEBDB   -  If set, will use .bdb file if user has privilege
;   REORDER  -  If set, put low energy channel (index 4) into index 0,
;		channel indices 0-3 are moved to 1-3.  This properly
;		orders the channel indices by energy response.  The original
;		channel 4 is a residual channel just below the DISCLA1 (Msfc nomenclature)
;		energy.
; MODIFICATION HISTORY:
;
;	richard.schwartz@gsfc.nasa.gov, 26-sep-1997.
;	Version 2, fixed overflow time selection, richard.schwartz@gsfc.nasa.gov, 24-feb-1998.
;-
;*******************************************************************
;
pro read_discla, seconds, discla, cosines, livetime=live_time, filename=filename, $
             starttime=st_input, endtime=en_input, cos_spec=cos_spec,  $
             flare=flare, reorder=reorder, verbose=verbose, error=error, usebdb=usebdb
;
;*******************************************************************
!quiet=1
;

checkvar, verbose, 1   ; set verbose to true if not passed as keyword.

error = 0

;print,'st_input,en_input = ',st_input,en_input
if keyword_set(st_input) then $
   sec_st_input = anytim(st_input,/utime)


if keyword_set(en_input) then $
   sec_en_input = anytim(en_input,/utime)

usefdb_old=fcheck(usefdb,0)
; find and open requested input file
if keyword_set(usebdb) then begin
	usefdb=0
	dd_type=0
endif else begin
	usefdb=1
	dd_type=1
endelse
if strpos(strupcase(fcheck(filename,'')),'FITS') ne -1 then begin
	st_input = 0.0d0
	read_discla_fits, filename, seconds, discla, start_time=st_input, end_time=en_input
	sec_en = max(seconds)
	
; See if any overflow intervals overlap with requested time interval.  If so,
; expand time interval to accumulate to include the 1 minute before the overflow
; started.
	new_sec_st = seconds(0)
	check_overflow, seconds(0), sec_en, nfound, ovr_sec_st, ovr_sec_en
	if nfound gt 0 then begin
   	print, 'Correcting ', nfound, ' overflow intervals.'
   	if st_input gt 0.0 then begin
	new_sec_st = min ([seconds(0), ovr_sec_st-60.])
	read_discla_fits, filename, seconds, discla, start_time=new_sec_st, end_time=en_input
	sec_en = max(seconds)
	endif
	
	endif
	gro_point,anytim(seconds(0),/yoh), solephut(seconds(0)), cos8, cos_spec
	discla = float(discla)	
endif else begin
fs_open, file=filename, flare=flare, time=st_input, dd_type=dd_type, $
   dd_open=dd_open, verbose=verbose, error=error
if error then begin
   print, 'Can''t find input file.'
   print, 'Must provide file name, flare number, or a time ' + $
          'interval in calling arguments.'
   goto, errorexit
endif

; If user supplied start or end time, set sec_st and sec_en in common fscom
; to those values.  Otherwise use start and end of file.
if keyword_set(st_input) then sec_st=sec_st_input else sec_st=dd_open.startsec
if keyword_set(en_input) then sec_en=sec_en_input else sec_en=dd_open.endsec
sec_st = sec_st > dd_open.startsec
sec_en = sec_en < dd_open.endsec
; See if any overflow intervals overlap with requested time interval.  If so,
; expand time interval to accumulate to include the 1 minute before the overflow
; started.
new_sec_st = sec_st

check_overflow, sec_st, sec_en, nfound, ovr_sec_st, ovr_sec_en
if nfound gt 0 then begin
   print, 'Correcting ', nfound, ' overflow intervals.'
   new_sec_st = min ([sec_st, ovr_sec_st-60.])
   endif


; ras, 27-mar-1996, not needed
; set up arrays to read data into, 
;setup_arrays

; read all requested data into the above arrays
;print,atime(sec_en)
read_dd, seconds, discla, hsk, dd_open=dd_open, startt=new_sec_st, endt=sec_en, $
   error=error
free_lun, abs(dd_open.lun)
maxindex = n_elements(seconds)-1
if error then goto,errorexit
                                                       
; initialized the time array to 0, then filled in the time and data arrays 
; using time as an index into the array).
gaps = where (seconds eq 0.d0, count)
if count ne 0 then seconds(gaps) = seconds(0) + 1.024*gaps

pseconds = seconds(lindgen(n_elements(seconds)/2)*2)
maxposindex = n_elements(hsk(5,*))-1

xra = hsk(5,*)
xdec= hsk(6,*)
zra=  hsk(3,*)
zdec = hsk(4,*)
radec_fill, pseconds, xra, xdec, zra, zdec, maxposindex, error=error
if error then begin
	print, 'Error, no good pointing data found. Select a longer interval.'
	goto, errorexit
endif

det_point, pseconds(0:maxposindex), xra(0:maxposindex), $
   xdec(0:maxposindex), p_ndx, cos_arr, cos_spec=cos_spec, $
   zra(0:maxposindex), zdec(0:maxposindex), bad_data=bad_data

if bad_data then begin
   print,'No valid pointing data in requested time interval.'
   goto,errorexit
endif

; Save just one set of cosines for 8 detectors
cos8 = cos_arr (*,0)
cos_spec = cos_spec(*,0)

endelse
; If there were any overflows in any of the detectors, correct them
; Pass only non-zero elements to reset routine.  Reset will replace
; the overflowed valued with good values in the array temp. 

restore_overflow, discla, rotate(sort(cos8),2), maxindex=maxindex
;
;DISCLA5 (index 4) is now known to contain an independent
;signal below the threshold of DISCLA1.  Therefore we
;now subtract the others from it to obtain this signal for display.
;We have yet to determine the implication of this signal for
;livetime considerations.
discla(4,*,*) = discla(4,*,*) - total( discla([indgen(4),5],*,*),1)
wnzero = where( discla(5,0,*) ne 0, nz)
for i=0,7 do begin
        temp= (discla(4,i,*))(wnzero)
        jumper, temp, reset
        discla(4,i,wnzero) = reset
endfor
                                                                                                              

; Compute live time arrays for data.
livet = fltarr(8, maxindex+1)
for i=0,7 do livet(i,0:maxindex) = $
	1.024*livetime( reform(discla(0:3,i,0:maxindex)), det_id=i )
; Calculate rate.

; livetime will only be 0 where count_raw is 0, so set it to 1 wherever
; it is 0 just so we can divide by it.

bad_data = where (livet(*,0:maxindex) le 0, count_bad)
if count_bad ne 0 then livet (bad_data) = 1.
; Don't divide by livetime for channel 6 (CPD detector), just 1-5 (0-4)
for ich=0,4 do begin
   discla(ich,0,0) = discla(ich,*,0:maxindex) / livet(*,0:maxindex)
   if count_bad ne 0 then discla (6*bad_data+ich) = 0.
endfor

; just return the data that was requested (we might have started at
; an earlier time to include a complete overflow interval.
early = where (seconds(0:maxindex) lt fcheck(sec_st,seconds(0)), kearly)
minindex = 0
if kearly gt 0 then minindex = early(kearly-1) + 1
seconds = seconds(kearly:*)
discla = discla(*,*,kearly:*)
cosines = cos8
live_time = livet(*,kearly:*)

;   REORDER  -  If set, put low energy channel (index 4) into index 0,
;		channel indices 0-3 are moved to 1-3.  This properly
;		orders the channel indices by energy response.  The original
;		channel 4 is a residual channel just below the DISCLA1 (Msfc nomenclature)
;		energy.
if keyword_set(reorder) then discla= discla([4,indgen(4),5],*,*)

; Set base time for utplot package to start of day.  Subtract start of day
; from time array so it will be seconds relative to start of day.
; Also set start and end time for utplot package to 0 - find_dbfile sets
; these to start and end of the FDB file, which may not be correct here.
setutbase, anytim(seconds(0),/sec,/date)
setutstart, 0
setutend, 0
seconds = seconds - getutbase(0)

goto,getout

errorexit:
error = 1
goto,getout


getout:
usefdb=fcheck(usefdb_old,usefdb)

end
