[Spice_discussion] Implementing True of Date in the SPICE Kernels
Samuel H. Dupree, Jr.
sdupree at speakeasy.net
Sat Jan 21 01:28:56 PST 2006
Dear SPICE Kernel Users,
I'm working on a application that uses the coordinate transformation
facilities in the SPICE kernels. I need to be able to transform to true
equator and equinox of date (true of date or TOD). Now the frames.req
file gives a template for an inertial TOD frame definition. The frame
definition I came up with was
\begintext
Definition for an inertial Earth True Equator and Equinox of Date
frame.
The earth precession model is the 1976 IAU model.
The earth nutation model is the 1980 IAU model.
\begindata
FRAME_ECI_TOD = 1400010
FRAME_1400010_NAME = 'ECI_TOD'
FRAME_1400010_CLASS = 5
FRAME_1400010_CLASS_ID = 1400010
FRAME_1400010_CENTER = 399
FRAME_1400010_RELATIVE = 'J2000'
FRAME_1400010_DEF_STYLE = 'PARAMETERIZED'
FRAME_1400010_FAMILY = 'TRUE_EQUATOR_AND_EQUINOX_OF_DATE'
FRAME_1400010_PREC_MODEL = 'EARTH_IAU_1976'
FRAME_1400010_NUT_MODEL = 'EARTH_IAU_1980'
FRAME_1400010_ROTATION_STATE= 'INERTIAL'
\begintext
\begindata
BODY1400010_POLE_RA = ( 0. -0.641 0. )
BODY1400010_POLE_DEC = ( +90. -0.557 0. )
BODY1400010_PM = ( 190.16 +360.9856235 0. )
BODY1400010_LONG_AXIS = ( 0. )
\begintext
\begindata
BODY1400010_NUT_PREC_ANGLES = ( 125.045 -1935.5328
249.390 -3871.0656
196.694 -475263.3
176.630 +487269.6519
358.219 -35999.04 )
\begintext
Setting up the kernel file was not a problem. The problem I'm having is
validating my results.
For the validation I used the iau_PNM80 routine from SOFA (Standards of
Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). iau_PNM80 uses
the IAU 1976 Precession Model and the IAU 1980 Nutation Model. The
matrix returned rotates a vector in mean equator and equinox J2000 to
true equator and equinox of date. The rotation matrix returned by
iau_PNM80 for the true equator and equinox at B1950 (Julian Date
2433282.42345905) was:
9.999254945491462D-01 1.119498740037581D-02 4.865964215927001D-03
-1.119477515287043D-02 9.999373340309599D-01 -7.085431092707289D-05
-4.866452498682185D-03 1.637565659548431D-05 9.999881586158487D-01
On the other hand subroutine PXFORM from the SPICE kernels was called as
follows:
id = 1400010
et = 86400.d0 * ( B1950() - J2000() )
call PXFORM( 'J2000', 'ECI_TOD', et, Rot_SPICE(1,1) )
The rotation matrix returned was
9.999255120142722D-01 1.119367477663291D-02 4.865394947010554D-03
-1.119347855741818D-02 9.999373487741916D-01 -6.755904572630151D-05
-4.865846358039163D-03 1.309331937631033D-05 9.999881616138190D-01
where the difference between the two matrices ( SPICE minus SOFA ) was
1.746512601208394D-08 -1.312623742900920D-06 -5.692689164470133D-07
1.296595452247862D-06 1.474323185757243D-08 3.295265200771379D-06
6.061406430217520D-07 -3.282337219173985D-06 2.997970094398283D-09
I expected the differences to be much smaller. I expected to see
differences consistent with two floating point numbers x and y being
computationally equal to 15 digits. More on this later.
Not seeing smaller differences, I looked at the rotation matrix
generated by something as simple as a transformation from J2000 to
B1950. From SOFA I compared the rotation matrix returned from iau_PREC76
(the IAU 1976 Precession Model) with the rotation matrix returned from
subroutine IRFTRN (SPICE kernels). Subroutine iau_PREC76 returned
9.999257079523629D-01 1.117893812642771D-02 4.859003841454430D-03
-1.117893813777015D-02 9.999375133499889D-01 -2.715792625851089D-05
-4.859003815359272D-03 -2.716259471424703D-05 9.999881946023741D-01
while subroutine IRFTRN returned
9.999257079523629D-01 1.117893812642769D-02 4.859003841454428D-03
-1.117893813777014D-02 9.999375133499887D-01 -2.715792625851078D-05
-4.859003815359270D-03 -2.716259471424704D-05 9.999881946023741D-01
Computing the difference between the two subroutines (again, SPICE minus
SOFA) gave
0.000000000000000D+00 -1.734723475976807D-17 -1.734723475976807D-18
1.561251128379126D-17 -1.110223024625157D-16 1.118083490375676D-19
1.734723475976807D-18 -1.016439536705160D-20 0.000000000000000D+00
I went further and checked the (i,j) elements from these two matrices to
see if they were computationally equal to 15 digits. The following test
was applied
abs(x, y) .le. max( abs(x), abs(y) ) * eps, where eps = 1.d-15
The result was
0.000000000000000D+00 -1.734723475976807D-17 0.000000000000000D+00
1.561251128379126D-17 0.000000000000000D+00 1.118083490375676D-19
0.000000000000000D+00 0.000000000000000D+00 0.000000000000000D+00
In other words, the rotation matrix returned by iau_PREC76 for a J2000
to B1950 transformation is not computationally equal (to 15 digits) to
the rotation matrix returned by IRFTRN. Again, I expected the matrices
returned by the two routines to be computationally equal to 15 digits.
It is at this point I ask:
* Have I overlooked something in formulating the frame definition
for an inertial TOD transformation?
* Have I missed something in comparing the matrices returned by
iau_PREC76 and IRFTRN?
* Lastly, whose results are correct, SOFA or SPICE?
Any suggestions?
Sam Dupree.
More information about the Spice_discussion
mailing list