[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


     Definition for an inertial Earth True Equator and Equinox of Date
     The earth precession model is the 1976 IAU model.
     The earth nutation model is the 1980 IAU model.


        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_PREC_MODEL    = 'EARTH_IAU_1976'
        FRAME_1400010_NUT_MODEL     = 'EARTH_IAU_1980'



        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.                        )



        BODY1400010_NUT_PREC_ANGLES  = (  125.045    -1935.5328
                                          249.390    -3871.0656
                                          196.694  -475263.3
                                          176.630  +487269.6519
                                          358.219   -35999.04     )


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

      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