; docformat = 'rst'
;+
;The function precess converts a GPS microsecond time to rotational 
;matrices of quaternions and then combines the rotational matrices.
;
;NOTES: The constants, used by zetaAPoly, thetaAPoly, and zAPoly in 
;precess, are computed from (s / 3600) * !const.dtor where s has the values of 
;0, 2306.2181, 2004.3109, 0.30188, -0.42665, 1.09468, 0.017998, and 
;-0.041833, 0.018203 respectively. These numbers are all used to calculate 
;the corrections (in arcseconds) needed in order to correct for both 
;luni-polar precession as well as planetary precession which is also known 
;as general precesssion. General precession results from the gravitational 
;attractions of the moon, sun, and planets causing the Earth's rotational 
;axis to slowly gyrate about the vertical axis.
;
;The funcion storadians takes a latitude / longitude in seconds and 
;converts it to degrees (s / 3600d) and maps it to radians by multiplying 
;the result by the IDL constant !const.dtor.
;
; :Examples:
;   comb_rot_mats = precess(gps_time
;   radians = storadians(sector)
;-
;+
; :Params:
;   s : in, required, type=double
;     The sector to be converted into radians.
;     
; :Returns:
;   Returns the radian conversion of the sector.
;-
function sToRadians,s

  if ~real_numeric_type(s) then begin
    doc_library, 'precess'
    return, -1
  endif
  return, (s / 3600d) * !const.dtor
end

;+
; :Params:
;   usec : in, required, type=double
;     The GPS microsecond used for the calulation.
;     
; :Returns:
;   Returns the combined rotational matrices of quaternions.
;-
function precess,usec

  if ~real_numeric_type(usec) then begin
    doc_library, 'precess'
    return, -1
  endif
  
  zetaAPoly= [0, 0.011180860378820828d,         0.0000014635555389981416378d,  0.0000000872567685924788776d]
  thetaAPoly=[0, 0.0097171734433392966540088d, -0.0000020684574947870930276d, -0.0000002028120978660714530d]
  zAPoly=    [0, 0.011180860378820828d,         0.0000053071581289374870959d,  0.0000000882506317821132837d]
  JD=usec2jd(usec)
  
  ;input time normalized to number of Julian days since 01-Jan-2000 
  ;12:00:00.00 divided by 365.25 (days) multiplied by 1000.
  TT=(JD-2451545d)/36525d; 
  zeta=poly(TT,zetaAPoly);
  z=poly(TT,zAPoly);
  theta=poly(TT,thetaAPoly);
;  print,'zetaA (\"):        ',zeta*180d/!dpi*3600d
;  print,'zA (\"):           ',z*180d/!dpi*3600d
;  print,'thetaA (\"):       ',theta*180d/!dpi*3600d;
  return,combine_rotation(combine_rotation(quatx(3,zeta),quatx(2,-theta)),quatx(3,z));
end