Simtel MSDOS - Coast to Coast/simteldosarchivecoasttocoast2.iso/astrnomy/de118i.zip/SSYSTEM.DOC

Text File | 1992-09-07 | 19KB | 425 lines

            DE118i.ARC

Summary:

This program performs an N-body numerical integration of the Sun,
Earth, Moon, and planets. The physics model includes
gravitational harmonics of the Earth and Moon, Earth tides, Lunar
librations, relativity theory corrections, and five asteroids.
The integrator is an Adams-Bashforth-Moulton predictor-corrector
method that is initialized by several Runge-Kutta steps.

This version "i" is extremely accurate due to use of long double
precision arithmetic throughout.  It has been used to compute
ephemerides extending from 9,000 B.C. to 13,000 A.D.

An effort has been made to achieve exactly the same physics
model as in the NASA-JPL integrations that form the basis of the
Astronomical Almanac.  The numerical results compare favorably
with JPL's DE118 and DE200 ephemerides.


The Model:

The principal forces that determine the motions are Newtonian
gravitational attractions between the solar system bodies modeled
as point masses or spheres.  The Newtonian forces are modified by
relativistic corrections as described in the first reference
(Newhall et al, 1983).  The relativistic center of mass of the
solar system is  constrained to stay at the origin of the
coordinate system by adjusting the barycentric location and
velocity of the Sun after each integration step.

The five asteroids Ceres, Pallas, Vesta, Iris, and Bamberga are
included in JPL's model because of their perturbations on the
motion of Mars.  Their motions are not integrated; the orbits are
calculated by Chebyshev expansions that approximate fixed
Keplerian ellipses.  Coefficients for these expansions are
supplied in the documentation records on the JPL ephemeris tapes.
The asteroids have a numerically significant effect (10^-12 au)
on the location of the barycenter of the solar system.  Therefore
it is important to include them in order to get an accurate
short-term reproduction of JPL's ephemerides.  The asteroid
orbits do not agree with the Russian Ephemerides of Minor
Planets.

Also not integrated is the Earth's orientation.  This is computed
by a precession matrix and a one-term nutation series.  The
orientation is required for oblateness calculations. In the
DE102, and apparently also the DE118, the precession was
calculated using Newcomb's precession constant and as if the
reference equatorial coordinate system were that of B1950. The
method used to calculate nutation is uncertain; however, the 1977
and 1980 IAU formulas give a good fit, at least in the tested
neighborhood of the initial epoch.  The 10,000-year precession
formulas of J. Laskar (1986) and the 1980 nutation formula are
the ones actually distributed with DE118i.ARC.

Lunar librations are computed by integrating Euler's equations
for the rotation of a rigid body, expressed in terms of the
Euler angles that relate the Moon's orientation to the fixed
reference equatorial coordinate system.  The rotation is
affected by external torques from the Earth and the Sun. Two of
the numerical parameters are ambiguous or not available in the
documentation.  In the case of the DE102 model, the mean radius
of the Moon is 1738.09 km.  In the DE118, the Lunar moment of
inertia coefficient C/MR^2 was not specified; its value has been
estimated to be 0.390689526131941 by fitting to the initial
second derivatives of the libration angles.

The DE200 is identical to the DE118 except for a rotation of the
computed states from B1950 to J2000.  Provision is made in this
program to perform the rotation prior to writing out the
results.  The default configuration produces DE200 outputs.


The Integrator:

Numerical integration of the system of ordinary differential
equations is carried out by an Adams-Bashforth-Moulton
predictor-corrector method.  This has to be initialized by
supplying a number of state variable points that have been
calculated by some other means (or else by the same method
but with very small initial steps).

In the integrator, if the Adams order is set to N, then the
first N+1 steps are performed by the Runge-Kutta integrator
subroutine runge.c.  These take about 7 times longer than the
subsequent steps done by the Adams-Bashforth-Moulton integrator
routine, adams4.c.

Neither the step size nor the approximation order has been made
self-adjusting.  Some experimentation was required to achieve a
balance between speed and accuracy.  The integrator settings
adopted for most of the tests were an Adams order of 11 and step
size of 3/16 day; these are near the optimum values for the
Moon's motion using IEEE 754 double precision floating point
arithmetic.  However, the Lunar librations are more accurate if
the step size is smaller. The numerical accuracy can be improved
by using IEEE 854 (80-bit) arithmetic and a step size of 1/8
day.  The software has been written so that the parameter
LDOUBLE in the file prec.h can select the extended precision
arithmetic, for language compilers that support such a feature.
Instability in the Moon's motion appears for step sizes of about
1/4 day or greater and Adams order >= 11.

JPL's integrator also uses an Adams method, but theirs has the
capability to modify the step size independently and adaptively
for all the integrated variables.  This ability makes it
potentially faster and more accurate than an integrator with
fixed step size. See the referenced book by Shampine and Gordon
and the paper by Krogh for more details.


Running the Program:


When ssystem.exe is invoked, it asks several questions that
require operator keyboard entries.

(1)
Do you want to restore a previously saved integrator state?

If starting a new ephemeris, answer by typing the letter n
and the Enter key.

To continue a computation that was interrupted, answer y.  In
that event a valid checkpoint file named de118.sav must be
available.  The program creates this file whenever it completes
a run or is interrupted gracefully.

To continue, but with a different tabular output interval (see
question (4) ), answer z.


(2)
Adams order ?

Asked only if the answer to (1) was n.  This is the integer
degree of the interpolation polynomial employed by the Adams -
Bashforth - Moulton integrator formula.  A small value gives less
accurate results than a larger value; but if too large, the
computation will be unstable and diverge to a meaningless result.
Maximum allowed value is 13. Suggested answer is the number 12
followed by Enter.

(3)
step size, days ?

Asked only if (1) was n.  This is the duration of each internal
integrator step.  A smaller step size gives a smaller error, down
to a minimum error limit imposed by the accuracy of the
arithmetic.  For accurate representation of the Lunar orbit, a
suggested answer is the number 0.125.

(4)
Interval between output samples, in steps ?

Asked if (1) was n or z.  This controls the tabulation interval
of the output ephemeris. For a suggested answer 3200, one
ephemeris record will be written for every 3200 integrator steps,
or 0.125 x 3200 = 400 days.

(5)
Output file name ?

Enter the name of a disc file.  The program will create this file
and write the ephemeris records into it.  Caution: if the file
already exists, its previous contents will be destroyed.

(6)
Number of output samples ?

Enter the desired number of records to be written to the
ephemeris file.  The program will stop, close files and exit
after it has computed the indicated number of records.
Be sure to rename the checkpoint file de118.sav that is
produced, otherwise it will be overwritten the next time
the program is run.


At any desired time, the user may interrupt the integration by
typing a capital S and the Enter key.  Computation will continue
until the next record has been written to the output file. The
file will then be closed and the complete internal state of the
integrator will be saved in a file named de118.sav.  Note,
the ability to detect that the user hit a key on the keyboard
is system dependent; this functionality has been included
only for IBM PC, using the library routine called kbhit().

It is important to retain all the .sav files so that it will be
possible to rerun interesting intervals later at a higher time
resolution. Copy each .sav file produced to a safe place,
renaming it by any convenient serial ordering scheme.  In the
event of equipment or power failure, the integration may be
continued from the last .sav file, but this file must be renamed
to have the name de118.sav. Note, the integrator does not merely
"restart" from this information; it carries on as if it had never
been interrupted at all.

The integrator output ephemeris files must be given different
names also, as the program will overwrite and destroy the data if
it is told a pre-existing file name for output.  Several output files
may be concatenated into one large file by the ordinary MSDOS
command of the form "copy /b name1+name2 name3". It is suggested
that the files be limited to 1.2 megabytes in size so they can be
copied onto floppy discs.

On startup, the program reads initial conditions and certain
constant parameters, such as the speed of light, from a text
file named aconst.h.  This file may be modified with a text
editor to explore the results of varying the initial conditions.


Numerical Comparisons:


In researching the model and debugging the program, the
discrepancy in the initial inertial acceleration of the Moon was
taken as an indication of the closeness of fit between this
program and JPL's original calculation. This discrepancy is less
than 10^-18 astronomical units (au)  per day^2 in each coordinate.

A complete set of starting conditions has been published for the
DE102. Using the DE102 model, the discrepancy in the geocentric
Moon location between this program and DE102 was 4x10^-14 au
after 400 days and 9x10^-12 au after 7200 days.  These results
are in agreement with the published estimate for the drift in
the DE102's Moon position, 10^-9 x days^1.7  kilometers.

The following table shows the result of estimating a set of
initial conditions for the DE200.  The initial libration state
was calculated by the analytical theory of Eckhardt, then
adjusted to obtain a least squares fit of the Moon's position
and velocity over a 6 year interval.  The table indicates the
maximum observed discrepancies in barycentric position vectors
between this program and DE200, measured in astronomical units,
over the indicated integration periods.  Comparisons were made
at 24-day intervals.

           408 days       7200 days
Mercury    1.2e-14         1.8e-11
Venus      7.0e-14         6.8e-12
EMB        3.5e-14         7.2e-13
Mars       2.0e-14         3.0e-13
Jupiter    1.0e-14         1.2e-12
Saturn     1.5e-14         2.9e-13
Uranus     5.0e-14         5.5e-13
Neptune    3.4e-14         5.4e-13
Pluto      6.2e-14         7.2e-13
Moon-Earth 2.8e-13         7.6e-12
Sun        9.6e-18         1.2e-15

The step size was 3/16 day and the Adams order was 11. The
7200-day calculation in IEEE double precision arithmetic took
about 6 hours on a 12MHz 68020/68882 computer.

Original 18-digit DE118 state vectors, including the librations,
were kindly supplied by JPL.  These yielded the best comparisons
to the DE200 Moon, so the DE118 model is the only one
distributed with this program.  The effect of arithmetic
precision can be seen in the following comparison of barycentric
states to the DE118 after integrating forward 400 days with a
1/8 day step size.  "Double" refers to IEEE 754 arithmetic and
"Long Double" is the 80-bit mode of an 80287 arithmetic
coprocessor.  The third and fourth columns of the table are the
results of 7200-day integrations in double arithmetic. The DE118
model was used, the output was rotated to the DE200 equinox, and
comparisons were made at 24-day intervals to the DE200
ephemerides; the maximum discrepancy, in au, is reported. The
last column is the same as the fourth but the integration was
done in long double arithmetic.

Arithmetic  Double      Long Double   Double       Double  Long Double
Interval     400 d        400 d       7200 d       7200 d    7200 d
Step Size    1/8 d        1/8 d       3/16 d        1/8 d     1/8 d

Mercury    1.2e-13      3.4e-15      1.9e-12      2.4e-12     1.9e-12
Venus      4.0e-14      3.9e-15      1.1e-12      6.9e-13     1.3e-12
EMB        9.4e-14      2.1e-16      1.5e-12      3.0e-12     1.0e-12
Mars       4.6e-14      1.7e-15      1.0e-12      1.3e-12     8.7e-13
Jupiter    1.0e-14      5.9e-15      3.9e-13      4.2e-13     2.8e-13
Saturn     1.5e-14      7.0e-15      1.6e-13      2.1e-13     3.3e-13
Uranus     5.1e-15      2.4e-14      3.3e-13      4.7e-13     3.4e-13
Neptune    9.4e-14      3.0e-14      3.7e-13      5.3e-13     4.2e-13
Pluto      1.6e-14      2.6e-14      6.3e-13      2.7e-13     5.5e-13
Moon-Earth 1.5e-14      1.6e-14      2.2e-12      6.4e-13     2.0e-13
Sun        1.5e-17      5.9e-18      4.0e-16      4.1e-16     2.4e-16
Librations 2.7e-12      2.6e-12         ?            ?           ?

Note that the variable step size in JPL's integrator was
adjusted for an accuracy of 2 x 10^-16 au/day for the orbits and
2 x 10^-14 rad/day for the librations.  The former figure applied
linearly with time would predict a magnitude error after 7200
days of 2.5 x 10^-12 au for the orbits.  See Moshier (1992) for
further comparison with the DE200.


Computer Progam Files:

aconst.h        Text file of initial conditions read by the program.
                These may be changed without recompiling the software.
de118.sav       Complete state of the integrator saved after each run.
                The program may carry on later by reading this file.
prec.h        Sets arithmetic to double or long double precision
int.h        Integrator settings
mconf.h        Sets hardware (680x0, IBM PC, DEC PDP-11/VAX)
ssystem.h    Software parameter settings
ini118d.h    Compiled initial conditions and test states from DE118
aconst.c    Compiled initial conditions for long double precision
adams4.c    Adams-Bashforth integrator
runge.c        Runge-Kutta integrator
ssystem.c    Main program and force calculations
jplmp.c        Asteroid orbit Chebyshev expansions
oblate.c    Earth-Moon figure model and tides
precess.c    Precession to or from basic epoch
nut1t.c        One-term nutation series
epsiln.c    Obliquity of Earth's equator to the ecliptic
findcent.c      Adjusts initial state to place barycenter at origin
asinl.c         Circular arcsine function in long double precision
atanl.c         Circular arctangent function in long double precision
sinl.c          Circular sine function in long double precision
tanl.c          Circular tangent function in long double precision
zatan2.c        Quadrant correct arctangent
polevll.c       Polynomial evaluator for long double precision
mtherr.c        Common error handler
const.h         References to some constants
rdnums.c        Reads initial conditions from aconst.h
ieee.c          Extra precise arithmetic used to read decimal
                constants at run time from aconst.h
econst.c        Extended precision numerical constants for ieee.c
etodec.c        DEC/IEEE floating point format conversion
ssysteml.mak    Microsoft makefile for long double precision
ssystemd.mak    Microsoft makefile for regular double precision
unix.mak        Unix makefile
ssystem.opt     DEC VAX/VMS makefile
descrip.mms     DEC VAX/VMS makefile

Program Output:

The data structure for the numbers that are written into the
output file from ssystem.c is given in the file ini118d.h. Each
output file record contains the Julian date followed by the
current velocity and position states. The state variables are
equatorial rectangular cartesian coordinates referred to a fixed
equator and ecliptic (B1950 or J2000).  The origin is at the
barycenter of the Solar system, and the unit of distance is the
astronomical unit. The unit of time is the day (86400 ephemeris
seconds).  All the numbers are floating point, either double or
long double precision.  There are six numbers for each object,
in the following order:

dx/dt, x, dy/dt, y, dz/dt, z

The vector (x, y, z) is the barycentric position vector of the
object; the vector (dx/dt, dy/dt, dz/dt) is the barycentric
velocity vector.

The order of the objects is Librations, Mercury, Venus, EMB,
Mars, Jupiter, Saturn, Uranus, Neptune, Pluto, MOON, Sun.

EMB is the arithmetic average of Earth and Moon, weighted by
their masses. The coordinates of the MOON vector are the the
solar system barycentric coordinates of the Moon minus those of
the Earth.

The Librations "object" is used to find the the selenocentric
direction of the Earth relative to the Lunar principal moment of
inertia axes. The coordinates integrated are the three Euler
angles, in radians, of the Moon's orientation relative to the
fixed equatorial coordinate system.  See the DE102 paper given in
the references for further explanation.


Acknowledgments:

  I thank J. G. Williams and E. M. Standish, of JPL, and R.
Babcock, of SAO, for their help in clarifying the details of the
model.  These people of course had nothing to do with writing
this program and bear no responsibility for my mistakes.


References:

Eckhardt, Donald H., "Theory of the Libration of the Moon,"
The Moon And the Planets 25, 3-49 (1981).
The introduction explains Cassini's laws and the libration
equations.

Goldstein, Herbert, _Classical Mechanics_, Addison-Wesley, 1950.
Contains a useful tutorial on rigid body rotation.

Laskar, J., "Secular terms of classical planetary theories
using the results of general theory," Astronomy and Astrophysics
157, 59070 (1986).

Moshier, S. L., "Comparison of a 7000-year Lunar ephemeris with
analytical theory,"  Astronomy and Astrophysics (1992), in press.
The output of this computer program is compared in detail
with the analytical Lunar theory of Chapront and Chapront.

Newhall, X. X., E. M. Standish, and J. G. Williams, "DE102: a
numerically integrated ephemeris of the Moon and planets
spanning forty-four centuries," Astronomy and Astrophysics
125, 150-167 (1983).
This is the primary summary reference for the physics model.

Standish, E. M., "Orientation of the JPL Ephemerides, DE200/LE200,
to the Dynamical Equinox of J2000," Astronomy and Astrophysics
114, 297-302 (1982)
Provides rotation matrices for conversions between DE102, DE118,
and DE200.

Shampine, L. F. and M. K. Gordon, _Computer Solution of
Ordinary Differential Equations_, W. H. Freeman, 1975.
This book gives a detailed explanation of Adams integrators
with either fixed or variable step size.

Thomas, Benku, "The Runge-Kutta Methods," BYTE magazine
11, #4, April 1986


-Steve Moshier
 2 April 1991
Last update: 16 August 1992