[Spice_discussion] Implementing True of Date in the SPICE Kernels
Nat Bachman
Nathaniel.Bachman at jpl.nasa.gov
Mon Jan 23 02:49:07 PST 2006
Hi Sam,
A few comments regarding your questions:
1) Your ECI_TOD frame definition is correct; none of the
errors you found are due to that definition.
Note that the inertial/rotating status of the frame is
a matter of choice; I'm not saying that using the
"rotating" state specification would be incorrect.
2) There appears to be something amiss with your
comparison of the ECI_TOD frame the B1950 epoch
with that obtained from the SOFA software. When
I attempted to duplicate your results using a small
test program "sofa.f" (which is attached), I found
round-off level differences between results from
SPICE, using your frame definition, and SOFA (third matrix below).
The error magnitude you found in this comparison is indeed much too
large: it looks a lot like what you get if
you truncate the JD inputs to the SOFA iau_PNM80 routine to
the single precision level.
3) The differences between outputs from IRFTRN and those from
iau_PREC76 are *slightly* larger than one might expect.
This is due to some loss of precision in the implementation
of IRFTRN. However, the error is many orders of magnitude
smaller than the precision level of the constants used to
define the IAU 1976 precession model, so I don't see that
the difference is computationally significant. We're
talking about angular errors on the order of 1e-16 radians.
The polynomial coefficients used to define the Euler angles
for the IAU 1976 precession model are specified to at most eight
digits; the implied precision level is coarser than 1.e-11
radians/century.
A note on SPICE implementation: the precession implementation used
to support your dynamic frame definition is not taken from IRFTRN,
but rather from the private SPICE routine ZZEPRC76. As
you can see in the attached output from the test program
sofa.f, the agreement between that private routine and
the corresponding SOFA routine at the B1950 epoch is (in this test,
at any rate) perfect. This result is of course platform-dependent.
4) The differences in the SPICE and SOFA implementations of the
IAU 1976 precession and IAU 1980 nutation models are at noise
level.
I'm glad to see that you're carefully testing your frame definition
kernel. This is an extremely important step, particularly if the
kernel is intended for widespread distribution.
Best regards,
Nat Bachman (NAIF/JPL)
Nathaniel.Bachman at jpl.nasa.gov
Output generated by sofa.f
==========================
JD1 = 0.2433282423459050E+07
Sofa-Spice prec diff at B1950
0.0000000000000000E+00 0.0000000000000000E+00
0.0000000000000000E+00
0.0000000000000000E+00 0.0000000000000000E+00
0.0000000000000000E+00
0.0000000000000000E+00 0.0000000000000000E+00
0.0000000000000000E+00
Sofa nutation*precession matrix is
0.9999255120142722E+00 0.1119367477663290E-01
0.4865394947010550E-02
-0.1119347855741818E-01 0.9999373487741918E+00
-0.6755904572624971E-04
-0.4865846358039160E-02 0.1309331937628749E-04
0.9999881616138188E+00
Sofa-Spice nutation*precession diff at B1950
0.0000000000000000E+00 -0.1734723475976807E-17
-0.8673617379884035E-18
0.0000000000000000E+00 0.0000000000000000E+00
0.2767426045269250E-16
0.0000000000000000E+00 0.9873016033196125E-17
0.0000000000000000E+00
Samuel H. Dupree, Jr. wrote:
> 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.
>
>
>
> _______________________________________________
> Spice_discussion mailing list
> Spice_discussion at naif.jpl.nasa.gov
> http://naif.jpl.nasa.gov/mailman/listinfo/spice_discussion
-------------- next part --------------
PROGRAM SOFA
IMPLICIT NONE
DOUBLE PRECISION SPD
DOUBLE PRECISION J2000
DOUBLE PRECISION B1950
DOUBLE PRECISION ET
DOUBLE PRECISION JD1
DOUBLE PRECISION JD2
DOUBLE PRECISION DIFF ( 3, 3 )
DOUBLE PRECISION SPMAT ( 3, 3 )
DOUBLE PRECISION SNPMAT ( 3, 3 )
DOUBLE PRECISION PRECXF ( 6, 6 )
DOUBLE PRECISION NUTXF ( 6, 6 )
DOUBLE PRECISION NPXF ( 6, 6 )
DOUBLE PRECISION PM ( 3, 3 )
DOUBLE PRECISION NPM ( 3, 3 )
INTEGER I
INTEGER J
JD1 = B1950()
WRITE (*, '(1X,A,D25.16)') 'JD1 = ', JD1
JD2 = 0.D0
ET = SPD() * ( (JD1 + JD2) - J2000() )
ET = SPD() * ( (B1950() + JD2) - J2000() )
CALL ZZEPRC76 ( ET, PRECXF )
CALL ZZENUT80 ( ET, NUTXF )
CALL MXMG ( NUTXF, PRECXF, 6, 6, 6, NPXF )
DO I = 1, 3
DO J = 1, 3
NPM(I,J) = NPXF (I,J)
PM (I,J) = PRECXF(I,J)
END DO
END DO
CALL IAU_PMAT76 ( JD1, JD2, SPMAT )
CALL VSUBG ( SPMAT, PM, 9, DIFF )
CALL M33 ( 'Sofa-Spice prec diff at B1950', DIFF )
CALL IAU_PNM80 ( JD1, JD2, SNPMAT )
CALL M33 ( 'Sofa nutation*precession matrix is', SNPMAT )
CALL VSUBG ( SNPMAT, NPM, 9, DIFF )
CALL M33 ( 'Sofa-Spice nutation*precession diff at B1950', DIFF )
END
SUBROUTINE M33 ( TITLE, XFORM )
IMPLICIT NONE
CHARACTER*(*) TITLE
DOUBLE PRECISION XFORM ( 3, 3 )
INTEGER I
INTEGER J
CALL TOSTDO ( ' ' )
CALL TOSTDO ( TITLE )
CALL TOSTDO ( ' ' )
DO I = 1, 3
WRITE (*, '( 3(1X,E25.16) )' ) ( XFORM(I,J), J=1,3 )
END DO
CALL TOSTDO ( ' ' )
END
More information about the Spice_discussion
mailing list