;docformat='rst'
;+
;   Get raw quaternion data from the database and calculate quaternion values
;
; :Author:
;    Chris Jeppesen, Bill Barrett
;
; :Params:
;    start_time : in, required, long64 or double
;      GPS seconds of start time for query
;    end_time : in, required, long64 or double
;      GPS seconds of end_time time for query.
; :Returns:
;    A structure as follows:
;      sec - Pointer to 1D array of quaternion timestamps, in GPS seconds
;      q   - Pointer to 2D array of quaternions. First index is quaternion component [0=x,1=y,2=z,3=scalar]
;            second corresponds to timestamp array above
;-
function loadquat, start_time=start_time,end_time=end_time,verbose=verbose

  ; The following three lines are correct for all real AIM data
  sid=1
  cycle=1
  external_element = 'ac'
  
  quaternion_columns = ['aqinr2bdyest0','aqinr2bdyest1','aqinr2bdyest2','aqinr2bdyest3']
  n_quaternions = n_elements(quaternion_columns)
  rows_from_interval = 2
  
  ; Note that for getting the time range (minmaxtimes) only the first quaternion column is considered.
  ; However, all four columns are retrieved in order to do the quaternion calculation
  min_and_max = minmaxtimes(external_element, quaternion_columns[0], $
    start_time, end_time, rows_from_interval, sid, cycle)

  ; 0 return for either of the min and max times means no data is present
  if min_and_max[0] gt 0 && min_and_max[1] gt 0 then begin
    one_second_time_tags = 10d6
    query_result = get_aim_hk_telemetry(SID=sid, sct_cycle=cycle, externalElement=external_element, $
    itemName=quaternion_columns, /gps, min_and_max[0], min_and_max[1] + one_second_time_tags, verbose=verbose)
  endif else query_result = 0
  if size(query_result,/type) ne 8 then begin
    print, get_routine_name() + 'No quaternions returned for time range: ' + usec2vms(start_time) + $
      ' to ' + usec2vms(end_time)
    return, -1
  end
  
  ; Get the timetags from the returned quaternion columns
  time_tags=(*(query_result.(0))).housekeeping[*].timetag

  ; n is the number of rows per quaternion
  n=lon64arr(n_quaternions)
  quats_from_db=ptrarr(n_quaternions)
  ; An IDL specific way of sequencing through the four items
  ;   in the data structure returned from the database via get_aim_hk_telemetry
  ; (*query_result.AC$AQINR2BDYEST0).housekeeping) through
  ;   (*query_result.AC$AQINR2BDYEST3).housekeeping)
  for i=0, n_quaternions - 1 do begin
    n[i]=n_elements((*(query_result.(i))).housekeeping[*].timetag)
    quats_from_db[i]=ptr_new((*(query_result.(i))).housekeeping)
  end

  str_free, query_result
  query_result=0

  ; This is checking the length of each of the time tag arrays.
  ; If they are all the same then 'nw' will be same as 'n_quaternions'
  w=where(n eq max(n),nw)
  ; If the number of rows per quaternion have different lengths, normalize everything to the
  ;   same length
  if nw ne n_quaternions then begin
    i=0L
    while i lt max(n) do begin
      time_tags=[(*quats_from_db[0])[i].timetag, (*quats_from_db[1])[i].timetag, $
        (*quats_from_db[2])[i].timetag, (*quats_from_db[3])[i].timetag]
      ;If a column is missing an entry here, it won't have a gap. Instead,
      ;the timestamp of this entry is greater than that of the others.
      ;Since we can't reconstruct the broken entry, we have to delete the others
      w=where(time_tags eq min(time_tags),ww) ;This gives the indices of the non-missing columns
      if(ww ne 4) then begin
        ;They aren't all non-missing, some of these things are not like the others.
        ;Drop the rows from the non-missing ones.
        for j=0,ww-1 do begin
          ;build a new array using the old elements up to but not including this one
          ;then the elements after but not including this one
          *quats_from_db[w[j]]=[(*quats_from_db[w[j]])[0:i-1],(*quats_from_db[w[j]])[i+1:n[w[j]]-1]]
          n[w[j]]--
          heap_gc
        end
      end else begin
        ;They're all the same at this point, so we can continue on
        i++
      end
    end
    ;Rebuild the time tag list. Grab the list from the zeroth component. It shouldn't matter
    ;which is used, the whole point of this exercise is to make them all the same.
    time_tags=(*quats_from_db[0]).timetag
  end
  
  transposed_columns = [transpose((*quats_from_db[0]).eu), transpose((*quats_from_db[1]).eu), $
    transpose((*quats_from_db[2]).eu), transpose((*quats_from_db[3]).eu)]
  str_free, quats_from_db

  ;Edit out cases where two consecutive lines have same timestamp
  unique_times=where((abs(time_tags[1:n_elements(time_tags)-1]-time_tags[0:n_elements(time_tags)-2])) ge 0.1)
  time_tags=time_tags[unique_times]
  transposed_columns = transposed_columns[*,unique_times]
  
  transposed_quats = transpose(transposed_columns)
  return, {quats, usec: ptr_new(time_tags), q: ptr_new(transposed_quats) }

end