;; File Version: 8.21
;; Last Modified: February 6, 2010
;; By: JEC

FUNCTION HSP_BINNER, data, bin, desc, time, aux_time, sec = sec

; Bin two dimensional HSP data arrays returned by HSP_READER and
; return as a one dimensional vector.
; Input: data - raw data read from DAPS file by HSP_READER
;        bin - number of points to bin the data by.
;        desc - description structure returned by HSP_READER
;        aux_time - structure of time values from the auxilliary time file
;        sec (keyword) - set to make value of bin correspond to number of 
;        seconds rather than number of points.
; Output: time - time vector, ET, corresponding to returned binned data vector.
; time values correspond to the middle of each integration period.
;
; November 1, 2007: Complete change to how times are calculated: get
; the full set of times for the marker points, then use interpol to
; get the full set of times.
; November 5, 2007: Agrees with Sremcevic times, except in first
; record where the disagreement is ~ 0.01 integration periods.
; January 1, 2009: Allow for aux_time = -1000 flagging a dummy
; auxiliary time file.
; January 16, 2009: Fix bug in dummy aux time file treatment.
; February 2, 2009: compute ET1 from header info if there is no valid
; pointer.
; Make "row" type LONG to avoid error when there are more than 32000 rows.
; August 26, 2009: Define packetsize for two different values of
; compression. Previously was hard-wired for compressed data
; (packesize=945).
; February 6, 2010: Prevent array out of bounds error at end
; calculating fulltime array for binned data.

ncol = LONG(N_ELEMENTS(data(*,0)))
nrow = LONG(N_ELEMENTS(data(0,*)))

CASE desc.compression OF
   0: packetsize = 532L
   1: packetsize = 945L
   ELSE: STOP,'Unknown compression.'
ENDCASE

IF KEYWORD_SET(sec) THEN binsize = bin*1000D0/desc.integration $
    ELSE binsize = bin

IF binsize LT 1. THEN BEGIN
    PRINT,'Binsize = ',binsize,' is invalid.'
    PRINT,'Integration period = ',desc.integration,' msec.'
    STOP,'Increase bin size.'
ENDIF
; Check the last row to see how many zero elements are there
; for buffering that need to be removed:
lastrow = data[*,nrow-1]
index = ncol
REPEAT index-- UNTIL (lastrow[index] GT 0) OR (index EQ 0)
nzero = ncol-index-1>0
print,nzero,' zeros removed.'

data_copy = data[0:N_ELEMENTS(data)-nzero-1]
npts = N_ELEMENTS(data_copy)

; Now bin it:
nbins = LONG(npts/binsize)

;;Leftover fractional bin at end:
leftover = FLOAT(npts)/FLOAT(binsize)-FLOAT(nbins)

bindata = ULONARR(nbins)
FOR ibin = 0L, nbins-1 DO bindata(ibin)=TOTAL(data_copy(ibin*binsize:$
    (ibin+1)*binsize-1))

;; Find the row of the aux_time data array which contains the timing
;; for the first record in this file:
IF aux_time.ticks[0] NE -1000. THEN BEGIN
   trow = -1L
   nticks = N_ELEMENTS(aux_time.ticks)
   REPEAT trow++ UNTIL (aux_time.ticks[trow] GE desc.sctime_sec_start) OR $
      (trow EQ nticks-1)
   IF aux_time.ticks[trow] NE desc.sctime_sec_start THEN BEGIN
      PRINT,'ERROR: The auxilliary time file does not match the input data.'
      STOP,'Stopping in HSP_BINNER.'
   ENDIF
ENDIF ELSE BEGIN
;; We don't have a valid auxiliary time file, so these are dummy values.
trow = 0
nticks = nrow+1
ENDELSE


;; Build unbinned time vector
fulltime = DBLARR(N_ELEMENTS(data_copy))

;; Approach in this version of HSP_BINNER is to get the times for the
;; individual points in each record that have valid times, record the 
;; the locations of those points, and then interpolate amongst all
;; those points to get the full time vector.

timepointvec = LONARR(nrow)
timevec = DBLARR(nrow)
row = 0L
half_int_period = 0.5*desc.integration/1000D0
REPEAT BEGIN
   aux_time_row = trow+row
   IF aux_time.ticks[0] EQ -1000. THEN BEGIN
      CSPICE_SCS2E, -82L, STRING(desc.sctime_sec_start), et1
      timevec[row] = et1
      timepointvec[row] = row*packetsize
   ENDIF ELSE IF aux_time.f1_valid[aux_time_row] AND NOT $
      (desc.integration EQ 1 AND aux_time.f1_ptr[aux_time_row] EQ 0) THEN BEGIN
      ;; Get the SCLK value at the first point in this row for which the S/C
      ;; RTI0 occurred. That is one tick after the tick corresponding to the
      ;; start of the packet:
      ticks1 = STRTRIM(STRING(aux_time.ticks[aux_time_row]+1L),2)
      ;; Get the ephemeris time of that point:
      CSPICE_SCS2E, -82L, ticks1, et1
      timevec[row] = et1
      timepointvec[row] = aux_time.f1_ptr[aux_time_row]+row*packetsize
   ENDIF ELSE BEGIN
      ;; no valid pointer for this row.
      CSPICE_SCS2E, -82L, STRING(desc.sctime_sec_start), et1
      timevec[row] = -1D
      timepointvec[row] = -1L
   ENDELSE
   row++
ENDREP UNTIL row EQ nrow OR aux_time_row EQ nticks-1


fullpts = LINDGEN(npts)
usethesepoints = WHERE(timevec GT 0.)
IF N_ELEMENTS(usethesepoints) GE 2 THEN $
   fulltime = INTERPOL(timevec[usethesepoints],timepointvec[usethesepoints], $
                       fullpts)-half_int_period ELSE $
                          fulltime = et1+DINDGEN(nbins)*half_int_period*2D0 -$
                                     half_int_period

npoints_to_use = nbins*binsize < N_ELEMENTS(fulltime)
time = REBIN(fulltime[0:npoints_to_use-1], nbins)

RETURN, bindata
END