Up to home page and other astronomy programs:
Get latest source code de118i-2.zip:
Get original version (1993-03-20) source code de118i.zip:
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.
Version de118i-2 can incorporate extra user-defined masses such as asteroids. (See the macro EXTRA_BODY and use aconst-aroid.h in place of aconst.h.)
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.
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.
When ssystem.exe is invoked, it asks several questions that require operator keyboard entries.
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.
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.
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.
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.
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.
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.
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.
Note: if your computer has a long double system math library, you should probably use it instead of these long double functions:
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.
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.
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 of DE118i: 3 March 1993
Last update of DE118i-1: 8 August 2000
Last update of DE118i-2: 13 May 2004
Up to home page:
Get original version (1993-03-20) de118i.zip:
Get previous version de118i-1.zip:
Get latest version de118i-2.zip:
The URL of this file is http://www.moshier.net/index.html.