; docformat = 'rst'
;+
;This function returns the spacecraft orientation quaternion at the given 
;input time.
;
; Notes
; * "Inertial" frame in a precise sense is the True Equator, Mean Equinox
;   (TEME) system used natively by SGP4. Origin is at Earth center of mass,
;   X axis points towards mean equinox in plane of true equator of date,
;   Z axis is along true axis of rotation normal to true equator, Y axis
;   completes a right handed set (points towards 6h right ascension).
; * ECEF frame is the above inertial frame rotated by the greenwich mean
;   siderial time angle. Origin is at Earth center of mass, X axis points
;   towards prime meridian, Z axis points towards north pole, Y axis
;   completes a right handed set (points at 90deg E)
;
; :Examples:
;   orient_quat = sc2global(gps_time, quat_table)
;-
;+
; :Params:
;   Time : in, required, type=double
;     Time in GPS microseconds (time unit of AIM telemetry).
;   quat : in, required, type=structure
;     Quaternion table. A structure of the following form: -sec: a pointer
;     to 1D array of doubles, representing GPS microseconds. Time must be 
;     in the range covered by sec. -q: a pointer to 2D array of doubles,
;     representing quaternion elements [-1,...,1] first index is 
;     quaternion element index (0=x, 1=y, 2=z, 3=scalar), second index is 
;     the time.
; quat_time : out, optional, double
;     The time tag of the output quaternion
;     
; :Keywords:
;   ecef : in, optional, type=boolean
;     Set to calculate the quaternion in ECEF frame. Default is J2000 \
;     inertial frame if flag is not set.
;   recurse : in, optional, type=int
;     Set to determine the minimum number of recursions through the quat
;     structure.
;   max_recurse : in, optional, type=int
;     Set to determine the maximum number of recursions through the quat
;     structure.
;     
; :Returns:
;   Returns a quaternion which can be used to transform a vector in the
;   spacecraft body frame to the global (ECI or J2000) frame selected 
;   above. With this quaternion q and a vector v_sc, in spacecraft frame,
;   do this: v_global = transformgrid(q, v_sc). 
;-
function sc2global,Time,quat,ecef=ecef,recurse=recurse,max_recurse=max_recurse, quat_time
  common orbit_data, sc2global_orbit_info

  if n_elements(max_recurse) eq 0 then max_recurse=3
  if n_elements(recurse) eq 0 then recurse=0
  ;Find correct time
  usec=(*(quat.usec))
  index = tpir(usec,Time)
  ;If we are past either end of the array, just move in so we look at, instead of past, the end.
  if index lt 1 then index=1
  if index gt n_elements(usec)-3 then index=n_elements(usec)-3
  if index lt 1 then index=1
  ;Get surrounding quaternions
  t0=usec[index-1]-usec[index];
  t2=usec[index+1]-usec[index];
  ;  t3=usec[index+2]-usec[index];  This is not used so I commented out - it was causing a problem when the # of usec was 3 instead of 4
  t=Time-usec[index];
  q0=(*(quat.q))[index-1,*]
  q1=(*(quat.q))[index+0,*]
  q2=(*(quat.q))[index+1,*]
  quat_time = usec[index]
  ;  q3=(*(quat.q))[index+2,*];;  This is not used so I commented out - it was causing a problem when the # of usec was 3 instead of 4

  ;interpolate and return result
  ;q=quat_normalize(quadterp(t0,q0,t1,q1,t2,q2,t)*0.5+quadterp(t1,q1,t2,q2,t3,q3,t)*0.5);
  if abs(t2) gt 60d6 then begin
    if recurse gt max_recurse then begin
      print, get_routine_name() + " Too long of a gap between quaternions. Asking to interpolate to "+ $
        usec2vms(t+usec[index])+", quat before "+usec2vms(usec[index])+" quat after "+usec2vms(t2+usec[index])
      return, -1
    endif else begin
      print, get_routine_name() + 'recursing, level = ' + strtrim(string(recurse), 2) + $
        ', time = ' + strtrim(string(long64(Time)), 2)
      recurse++
      if n_elements(sc2global_orbit_info) le 0 || $
        (Time lt sc2global_orbit_info.start_time || Time ge sc2global_orbit_info.stop_time) then begin
        sc2global_orbit_info = get_orbit_info(Time,Time)
      endif
      delta = sc2global_orbit_info.stop_time - sc2global_orbit_info.start_time
      return,sc2global(Time-delta,quat,ecef=ecef,recurse=recurse,max_recurse=max_recurse, quat_time)
    endelse
  endif

  q=slerp(q1,q2,linterp(0.0,0,t2,1,t))

  ;Converting J2000 to mean of date coordinates (close enough, within several arcsec of TEME)
  q=combine_rotation(q,precess(Time))

  if keyword_set(ecef) then begin
    return,combine_rotation(q,eci_to_ecef(Time))
  end else begin
    return,q
  end

end