Numerical Ephemeris Calculation

Abstract

This numerical core computes apparent topocentric ephemerides for asteroids, comets and interstellar objects from osculating orbital elements. Elliptic, parabolic and hyperbolic initial states are handled without an external orbit service. Numerical integration includes solar gravity, the first relativistic correction, planets, Earth and Moon, Earth J2, optionally 16 massive asteroids, and Marsden A1/A2/A3. Planetary and asteroid states are obtained from compact tables by cubic Hermite interpolation.

The observed direction follows UTC-to-TDB conversion, light time, Earth rotation, nutation and precession. Six extreme objects and 20 frozen JPL Horizons rows form the common acceptance test for all 49 independent target-language programs.

1. Terms and symbols

Every astronomical and numerical term required by the calculation is defined here. Important concepts are also explained briefly at their first use in the main text.

Glossary of astronomical and numerical terms.

Solar and terrestrial gravitational parameters are kept unambiguous.

2. Conventions, gravitational parameters and base quantities

The dynamical state y=[x,y,z,vx,vy,vz] is heliocentric, right-handed and referred to the equatorial J2000 frame. Position is in AU, velocity in AU/day and acceleration in AU/day squared. Internal angles are radians; public orbital elements, geographic inputs and result angles are degrees.

The symbol mu denotes a gravitational parameter mu=G*M, the product of gravitational constant G and body mass M. To avoid ambiguity this article always writes mu_sun=G*M_sun for the Sun and mu_E=G*M_E for Earth. For a perturbing body b, eta_b=mu_b/mu_sun=M_b/M_sun.

3. Initial state from osculating elements

Inputs are epoch t0, eccentricity e, perihelion distance q, perihelion epoch Tp, ascending node Omega, argument of perihelion omega and inclination i. Osculating means that these elements reproduce the position and velocity of the perturbed orbit at their epoch. With mu_sun=G*M_sun and Delta=t0-Tp, the true anomaly nu is determined first.

3.1 Ellipse

Validity: e < 1 - 1e-8
a = q / (1 - e)
n = sqrt(mu_sun / a^3)
M = normalize(n Delta, -pi, +pi)
E0 = M                    when e < 0.8
E0 = sign(M) pi           otherwise
E(k+1) = Ek - (Ek - e sin Ek - M) / (1 - e cos Ek)
Stop: |correction| < 2e-14 or after 80 iterations
nu = 2 atan2(sqrt(1+e) sin(E/2), sqrt(1-e) cos(E/2))

3.2 Hyperbola

Validity: e > 1 + 1e-8
a = q / (e - 1)                    (positive magnitude convention)
M = sqrt(mu_sun / a^3) Delta
H0 = asinh(M / e)
H(k+1) = Hk - (e sinh Hk - Hk - M) / (e cosh Hk - 1)
Stop: |correction| < 2e-14 or after 80 iterations
nu = 2 atan2(sqrt(e+1) sinh(H/2), sqrt(e-1) cosh(H/2))

3.3 Parabola and the near-parabolic boundary

Validity: |e - 1| <= 1e-8
s = sqrt(2 q^3 / mu_sun)
B = Delta / s
D = 2 sinh(asinh(1.5 B) / 3)
nu = 2 atan D

Barker's form avoids the ill-conditioned division by 1-e or e-1 in the immediate parabolic boundary. The following equations then apply to all three conics:

p = q (1 + e)
r = p / (1 + e cos nu)
v0 = sqrt(mu_sun / p)
r_orb = [r cos nu, r sin nu, 0]
v_orb = [-v0 sin nu, v0 (e + cos nu), 0]

Rotation from the orbital plane uses Rz(Omega) Rx(i) Rz(omega). For a two-dimensional input [x,y], the expanded form is:

X = (cosOmega cosomega - sinOmega sinomega cosi) x
    + (-cosOmega sinomega - sinOmega cosomega cosi) y
Y = (sinOmega cosomega + cosOmega sinomega cosi) x
    + (-sinOmega sinomega + cosOmega cosomega cosi) y
Z = (sinomega sini) x + (cosomega sini) y

A fixed x-axis rotation by 84381.448 arcseconds then transforms J2000 ecliptic position and velocity to the equatorial frame.

4. Perturber tables and Hermite interpolation

Each table node stores position r0 and velocity v0. For two nodes separated by h and u=(t-t0)/h, the core uses the following complete cubic Hermite basis and exact numerical coefficients:

h00(u) =  2 u^3 - 3 u^2 + 1
h10(u) =    u^3 - 2 u^2 + u
h01(u) = -2 u^3 + 3 u^2
h11(u) =    u^3 -   u^2

r(u) = h00 r0 + h10 h v0 + h01 r1 + h11 h v1

dh00/dt = ( 6 u^2 - 6 u) / h
dh10/dt =   3 u^2 - 4 u + 1
dh01/dt = (-6 u^2 + 6 u) / h
dh11/dt =   3 u^2 - 2 u
v(u) = (dh00/dt) r0 + dh10 v0 + (dh01/dt) r1 + dh11 v1

5. Equation of motion

d/dt [r, v] = [v, a(r,v,t)]

The solar base acceleration is a_s=-mu_sun r/|r|^3. A perturbing body at heliocentric position r_b and with mass ratio eta_b contributes the direct and indirect terms together:

a_b = mu_sun eta_b [(r_b - r)/|r_b - r|^3 - r_b/|r_b|^3]

Earth and Moon are treated separately inside 0.2 AU and as the Earth-Moon barycentre outside it. Mercury, Venus, Mars, Jupiter, Saturn, Uranus, Neptune and Pluto are always active. N16 is added only when the switch is set and all data are present.

5.1 Relativistic solar correction

a_rel = mu_sun / (c^2 |r|^3)
        * [(4 mu_sun/|r| - |v|^2) r + 4 (r dot v) v]

5.2 Earth J2

Inside 0.02 AU from Earth, the usual J2 acceleration is added as a direct-minus-indirect term in this port's approximate Earth-axis frame:

mu_E = G*M_E = 8.887692445125148326e-10 AU^3/day^2
k = 1.5 J2 mu_E R_E^2 / rho^5
a_J2,x = k x (5 z^2/rho^2 - 1)
a_J2,y = k y (5 z^2/rho^2 - 1)
a_J2,z = k z (5 z^2/rho^2 - 3)

5.3 Non-gravitational acceleration

g(r) = alpha (r/r0)^(-m) [1 + (r/r0)^n]^(-k)
Defaults: r0=2.808, alpha=0.1112620426, m=2.15, n=5.093, k=4.6142
a_ng = g(r) [A1 e_r + A2 e_t + A3 e_n]

e_r is radial, e_n is normal to the orbital plane along r cross v, and e_t=e_n cross e_r is transverse. When A1=A2=A3=0 the exact zero vector is returned.

6. Dormand-Prince 5(4)

An embedded Runge-Kutta pair produces fifth- and fourth-order solutions from shared stages; their difference estimates the local error. A Butcher tableau is the tabular collection of stage abscissae c_i, intermediate weights a_ij and final weights b_i. It therefore defines the complete integration step. All coefficients used by the core are:

c2=1/5;       a21=1/5
c3=3/10;      a31=3/40; a32=9/40
c4=4/5;       a41=44/45; a42=-56/15; a43=32/9
c5=8/9;       a51=19372/6561; a52=-25360/2187;
              a53=64448/6561; a54=-212/729
c6=1;         a61=9017/3168; a62=-355/33; a63=46732/5247;
              a64=49/176; a65=-5103/18656
c7=1;         a71=35/384; a72=0; a73=500/1113;
              a74=125/192; a75=-2187/6784; a76=11/84

Order 5: b = [35/384, 0, 500/1113, 125/192,
              -2187/6784, 11/84, 0]
Order 4: b*= [5179/57600, 0, 7571/16695, 393/640,
              -92097/339200, 187/2100, 1/40]

For every state component, scale=2e-15+2e-13*max(|y_old|,|y_new|). The normalized RMS error is sqrt(sum((error/scale)^2)/6). A step is accepted for error <=1 or |h|<=1e-8 day. The next factor is clamp(0.9*error^(-0.2),0.1,5), or 5 for zero error. Step size stays between 1e-8 and 0.25 day, and one propagation may take at most 1,000,000 steps.

7. Time scales and observer

The public time scale is UTC. Dynamics use the approximation TDB=UTC+(TAI-UTC+32.184 s)/86400. The implemented leap-second steps are:

JD >= 2457754.5 : TAI-UTC = 37 s
JD >= 2457204.5 : TAI-UTC = 36 s
JD >= 2456109.5 : TAI-UTC = 35 s
JD >= 2454832.5 : TAI-UTC = 34 s
JD >= 2453736.5 : TAI-UTC = 33 s
otherwise       : TAI-UTC = 32 s

With T=(JD_UTC-2451545)/36525, GMST and the complete short nutation and precession series used by the port are:

GMST = 280.46061837 + 360.98564736629 (JD-2451545)
       + 0.000387933 T^2 - T^3/38710000

L_sun  = 280.4665  + 36000.7698 T              [degrees]
L_moon = 218.3165  + 481267.8813 T             [degrees]
Omega_moon = 125.04452 - 1934.136261 T         [degrees]

Delta_psi = -17.20 sin Omega_moon -1.32 sin(2 L_sun)
            -0.23 sin(2 L_moon) +0.21 sin(2 Omega_moon) [arcsec]
Delta_eps =  9.20 cos Omega_moon +0.57 cos(2 L_sun)
            +0.10 cos(2 L_moon) -0.09 cos(2 Omega_moon) [arcsec]

eps0 = 23 + 26/60 + 21.448/3600 degrees
       -(46.8150 T + 0.00059 T^2 - 0.001813 T^3) arcsec

zeta = 2306.2181 T + 0.30188 T^2 + 0.017998 T^3 arcsec
z = 2306.2181 T + 1.09468 T^2 + 0.018203 T^3 arcsec
theta = 2004.3109 T - 0.42665 T^2 - 0.041833 T^3 arcsec

The geocentric site vector is formed from MPC rho factors or geodetic longitude, latitude and height, rotated by apparent sidereal time, then returned to J2000 through nutation and precession.

8. Light time and apparent direction

Earth and observer states remain at reception time. Emission time is iterated at most five times:

t_emit,0 = t_TDB
rho_k = r_object(t_emit,k) - r_earth(t_TDB) - r_observer(t_UTC)
t_emit,k+1 = t_TDB - |rho_k| * 0.0057755183
Stop: |t_emit,k+1 - t_emit,k| < 1e-12 day

RA  = normalize(atan2(rho_y,rho_x), 0, 2 pi)
Dec = atan2(rho_z, sqrt(rho_x^2+rho_y^2))

9. Binary tables

Planet chunks cover 2000 through 2040 in nine overlapping intervals at 0.5-day spacing. N16 also uses nine chunks at exactly 1.9999817474948436 days. Overlap permits safe Hermite interpolation at boundaries. Every value is little-endian float64.

10. Extreme examples

The following tables list UTC, frozen JPL Horizons reference values, the C-port result and the approximately spherical angular separation. Station A68 and identical nominal elements were used. Angles are degrees.

Orbit form and perihelion distance of the six validation cases.

Perihelion epoch and orbital-plane orientation.

Non-gravitational parameters; zero omits that component.

Frozen JPL Horizons reference positions.

Calculated positions and spherical separation from the reference.

433 Eros checks an ordinary ellipse. C/2023 A3 is near-parabolic and contains A1/A2. 2P/Encke is a periodic comet with non-gravitational acceleration. 2I/Borisov enforces e=3.356475782676596 with A1/A2/A3. Apophis and Duende test rapid apparent motion during close Earth passages.

11. Reproducibility and limitations

• The one-arcsecond limit applies to the same frozen nominal elements, not to a later orbit solution.
• The short nutation, precession and TDB approximation targets the documented interval and accuracy.
• Leap seconds after the last built-in step must be added before later epochs are used.
• The tables cover 2000 through 2040; other intervals require a new dataset.
• Uncertainty ellipses and covariance propagation are outside this nominal core.
• Find_Orb is neither runtime component nor fallback; it is used only for private local comparison.

Appendix A. Comparison of programming languages

Every successful edition passed the same 20 reference rows. Accuracy is sorted alphabetically by language. The separate timing table is sorted by ascending median time, fastest first. Each language calculated six cases seven times in an alternating cross-language order. These numbers show only relative performance on this test computer.

A.1 Accuracy, alphabetical

Maximum spherical angular separation over all 20 reference rows.

A.2 Performance, fastest first

Factor 1.000 denotes the fastest median in the validation run.

Appendix B. Licence boundary

The project-written numerical core, generators and documentation covered by LICENSE-SCOPE.md are under Apache License 2.0. Data tables, external ephemerides, JPL kernels and third-party programs require separate assessment. Redistribution must include LICENSE and NOTICE, preserve copyright and licence notices, and identify modified files. The patent grant and termination rule are stated in section 3 of the Apache licence.

Project website: https://www.jostjahn.de/amrum-software/ephemeris/