[Spice_discussion] Implementing True of Date in the SPICE Kernels
Nat Bachman
Nathaniel.Bachman at jpl.nasa.gov
Mon Jan 23 03:29:45 PST 2006
Hi Sam,
I've attached a better version of the test program sofa.f.
This one has some comments. :) Also, it explicitly uses
your frame kernel rather than just the low-level routines
that SPICE uses to implement your frame definition at run
time.
If you run this program, you'll see a very slight change
in the round-off errors---nothing significant. I've included
below the output I obtained.
-Nat
Program output
==============
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
Sofa - ECI_TOD nutation*precession diff at B1950
0.0000000000000000E+00 -0.1734723475976807E-17
-0.1734723475976807E-17
0.1734723475976807E-17 0.0000000000000000E+00
0.2767426045269250E-16
0.0000000000000000E+00 0.9874710099090633E-17
0.0000000000000000E+00
Nat Bachman wrote:
>
> 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
>
>
>
>------------------------------------------------------------------------
>
>
> 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
>
>
-------------- 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 )
DOUBLE PRECISION PXMAT ( 3, 3 )
INTEGER I
INTEGER J
C
C Load leapseconds and dynamic frame kernels. Dynamic frame
C kernel was provided by Sam Dupree.
C
CALL FURNSH ( 'lsk' )
CALL FURNSH ( 'tod.ker' )
C
C Set Julian date epochs compatible with the relevant SOFA
C subroutine interfaces.
C
JD1 = B1950()
WRITE (*, '(1X,A,D25.16)') 'JD1 = ', JD1
JD2 = 0.D0
C
C Convert to seconds past J2000 TDB for SPICE APIs.
C
ET = SPD() * ( (JD1 + JD2) - J2000() )
C
C Use SPICE private routines to find the nutation*precession
C matrix at B1950. This matrix gives the orientation of the
C earth true equator and equinox of date frame.
C
C Below, NPXF is a 6x6 state transformation matrix; NPM is
C the 3x3 rotation block of this matrix.
C
C PM is the 3x3 precession matrix.
C
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
C
C Get the SOFA precession matrix SPMAT. Compare to the SPICE
C matrix PM.
C
CALL IAU_PMAT76 ( JD1, JD2, SPMAT )
CALL VSUBG ( SPMAT, PM, 9, DIFF )
CALL M33 ( 'Sofa-Spice prec diff at B1950', DIFF )
C
C Get the SOFA nutation*precession matrix SNPMAT.
C
CALL IAU_PNM80 ( JD1, JD2, SNPMAT )
CALL M33 ( 'Sofa nutation*precession matrix is', SNPMAT )
C
C Compare to the SPICE version of this matrix NPM.
C
CALL VSUBG ( SNPMAT, NPM, 9, DIFF )
CALL M33 ( 'Sofa-Spice nutation*precession diff at B1950',
. DIFF )
C
C Repeat the above comparison using the high-level SPICE API
C routine PXFORM. Some additional round-off error may be
C introduced.
C
CALL PXFORM ( 'J2000', 'ECI_TOD', ET, PXMAT )
CALL VSUBG ( SNPMAT, PXMAT, 9, DIFF )
CALL M33 ( 'Sofa - ECI_TOD 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