TECHNICAL DOCUMENTATION
Numerical Ephemeris Calculation
Mathematical method, constants, series and cross-language validation
Version 0.5.3
July 25, 2026

Contents
 •	Abstract
 •	1. Terms and symbols
 •	2. Conventions, gravitational parameters and base quantities
 •	3. Initial state from osculating elements
 •	3.1 Ellipse
 •	3.2 Hyperbola
 •	3.3 Parabola and the near-parabolic boundary
 •	4. Perturber tables and Hermite interpolation
 •	5. Equation of motion
 •	5.1 Relativistic solar correction
 •	5.2 Earth J2
 •	5.3 Non-gravitational acceleration
 •	6. Dormand-Prince 5(4)
 •	7. Time scales and observer
 •	8. Light time and apparent direction
 •	9. Binary tables
 •	10. Extreme examples
 •	11. Reproducibility and limitations
 •	Appendix A. Comparison of programming languages
 •	A.1 Accuracy, alphabetical
 •	A.2 Performance, fastest first
 •	Appendix B. Licence boundary

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.
Term	Meaning in this article	
Ephemeris	Predicted sequence of sky positions and related quantities for specified epochs and an observer.	
apparent	Direction after light time, Earth rotation, nutation and precession; atmospheric refraction is outside this core.	
topocentric	Referred to the actual observer on or near Earth's surface, rather than Earth's centre.	
geocentric / heliocentric	Referred respectively to the centre of Earth or the centre of the Sun.	
barycentre	Common centre of mass of several bodies; EMB denotes the Earth-Moon barycentre.	
osculating orbital elements	Six elements of the instantaneous Kepler orbit matching position and velocity at an epoch.	
epoch	Instant for which orbital elements or a state vector are valid.	
conic	Ellipse, parabola or hyperbola; eccentricity e distinguishes the three orbit forms.	
eccentricity e	Dimensionless shape parameter: e<1 ellipse, e=1 parabola, e>1 hyperbola.	
q, Tp, i, Omega, omega	Perihelion distance, perihelion epoch, inclination, ascending-node longitude and perihelion argument.	
M, E, H, nu	Mean, eccentric, hyperbolic and true anomaly; angle or auxiliary value locating the body on its orbit.	
J2000	The fixed equatorial reference orientation at standard epoch JD 2451545.0 used by this core.	
AU	Astronomical unit, exactly 149597870.7 kilometres.	
JD	Julian date, the continuous astronomical day count.	
UTC / TAI / TDB	Coordinated Universal Time, International Atomic Time and barycentric dynamical time.	
leap second	Integer correction tying UTC to irregular Earth rotation while TAI remains uniform.	
gravitational parameter mu	Product G*M of gravitational constant and body mass; mu_sun is the solar value.	
mass ratio eta_b	Ratio M_b/M_sun of a perturbing body to the Sun, equal to mu_b/mu_sun.	
perturbing body	Body whose gravity adds acceleration beyond the pure two-body orbit.	
direct and indirect term	The body acts directly on the target; the simultaneous acceleration of the heliocentric origin is subtracted.	
N16	Optional set of the 16 most massive asteroids integrated as additional perturbers.	
cubic Hermite interpolation	Third-degree interpolation reproducing position and velocity at both table nodes.	
little-endian float64	Binary format storing the least significant byte first and IEEE-754 64-bit floating-point values.	
Runge-Kutta method	One-step method combining several derivative evaluations to advance the state.	
embedded method	Two Runge-Kutta solutions of different order share stages; their difference estimates local error.	
Dormand-Prince 5(4)	Embedded fifth- and fourth-order Runge-Kutta pair with seven stages.	
Butcher tableau	Tabular arrangement of Runge-Kutta coefficients c_i, a_ij and b_i defining stages and final weights.	
local error	Estimated error of one integration step, here the difference between fifth- and fourth-order solutions.	
adaptive step size	Automatic reduction or enlargement of the time step based on the normalized error estimate.	
RMS	Root mean square; here the square root of the mean of six squared, scaled component errors.	
clamp	Restriction of a number to a closed minimum-maximum interval.	
light time	Time required for light to travel from target to observer, shifting the evaluated emission epoch.	
precession / nutation	Slow secular and shorter-period changes in Earth's axis orientation.	
GMST	Greenwich mean sidereal time, the Earth rotation angle relative to the mean celestial frame.	
Earth J2	Leading coefficient of Earth's gravity-field oblateness relative to an ideal point mass.	
Marsden A1/A2/A3	Empirical comet acceleration in radial, transverse and normal directions.	
RA / declination	Right ascension and declination, the equatorial coordinate pair of the apparent direction.	
MPC rho factors	Dimensionless Minor Planet Center site factors for distance from Earth's axis and equatorial plane.	
arcsecond	1/3600 degree, denoted by the double-prime symbol.	
reference ephemeris / residual	JPL Horizons comparison position and the remaining spherical angular separation.	
covariance propagation	Propagation of statistical orbit uncertainties and correlations; outside this nominal core.	
acceptance test	Fixed comparison whose limit a port must meet before release.	

Solar and terrestrial gravitational parameters are kept unambiguous.
Symbol	Meaning	Value or unit	
mu_sun = G*M_sun	Heliocentric gravitational parameter of the Sun	0.0002959122082855911 AU cubed/day squared	
mu_E = G*M_E	Geocentric gravitational parameter of Earth	8.887692445125148326e-10 AU cubed/day squared	
mu_Moon = G*M_Moon	Gravitational parameter of the Moon	1.09318946230584718e-11 AU cubed/day squared	
eta_b = mu_b/mu_sun	Mass ratio of perturbing body b to the Sun	dimensionless	
r, v, a	Position, velocity and acceleration vector	AU; AU/day; AU/day squared	
t0, Tp, Delta	Element epoch, perihelion epoch and t0-Tp	JD or days	
q, a, p, e	Perihelion distance, semimajor axis, semilatus rectum and eccentricity	AU or dimensionless	
M, E, H, nu	Mean, eccentric, hyperbolic and true anomaly	radians	
Omega, omega, i	Ascending node, argument of perihelion and inclination	input degrees; internal radians	
c	Speed of light	173.144632674240 AU/day	
rho	Topocentric target vector, or Earth distance in the J2 section	AU	
J2, R_E	Earth oblateness coefficient and equatorial radius	dimensionless and AU	
A1, A2, A3	Marsden non-gravitational parameters	AU/day squared	
T	Julian centuries since J2000	dimensionless	

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.
Quantity	Exact value	Unit	
pi	3.141592653589793238462643383279502884	dimensionless	
AU	149597870.7	km	
mu_sun = G*M_sun	0.0002959122082855911	AU cubed/day squared	
mu_E = G*M_E	8.887692445125148326e-10	AU cubed/day squared	
mu_Moon = G*M_Moon	1.09318946230584718e-11	AU cubed/day squared	
speed of light	173.144632674240	AU/day	
light time	0.0057755183	day/AU	
Earth/Sun mass	3.0034896149156e-6	dimensionless	
Moon/Sun mass	3.6943033497651e-8	dimensionless	
Earth radius	6378.1363 / 149597870.7	AU	
Earth J2	1.08262668e-3	dimensionless	
J2000 obliquity	84381.448	arcseconds	

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
Perturbing body	Mass/Sun mass	
Mercury	1.660114153054348e-7	
Venus	2.447838287784771e-6	
Mars	3.227156037554996e-7	
Jupiter	9.547919101886966e-4	
Saturn	2.858856727222416e-4	
Uranus	4.366249662744965e-5	
Neptune	5.151383772628673e-5	
Pluto	7.350487833457740e-9	

N16 body	Mass/Sun mass	
Ceres	4.719154706681575e-10	
Pallas	1.0297360324108826e-10	
Juno	1.4471670474939703e-11	
Vesta	1.3026875547328177e-10	
Iris	8.589038999387399e-12	
Hygiea	4.238598614896659e-11	
Eunomia	1.5243642468487352e-11	
Psyche	1.1978215785019767e-11	
Euphrosyne	8.133159614634648e-12	
Europa	2.0216913527078825e-11	
Cybele	7.06871003275745e-12	
Sylvia	1.6337820878172953e-11	
Thisbe	8.965306556356798e-12	
Camilla	1.0878696848144257e-11	
Davida	2.9345275747768696e-11	
Interamnia	2.1327387533901138e-11	

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 1850-01-01 through 2150-01-01 in 60 overlapping intervals at 0.5-day spacing. N16 also uses 60 chunks at exactly 1.9999817474948436 days. Overlap permits safe Hermite interpolation at boundaries. Every value is little-endian float64.
Dataset	Interval	Step	Content	
Planets	1850-01-01–2150-01-01	0.5 day	11 bodies: EMB, Earth, Moon, Mercury through Pluto	
N16	1850-01-01–2150-01-01	1.9999817474948436 days	16 massive asteroids	

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.
Object	Epoch JD	e	q AU	
433 Eros	2453311.5	0.2228078944584026	1.133355399799004	
C/2023 A3	2460448.5	1.000095309222603	0.3914300307809727	
2P Encke	2459847.5	0.8477496967533629	0.3379482792219925	
2I Borisov	2458853.5	3.356475782676596	2.006520878500843	
99942 Apophis	2459215.5	0.19152168935010219	0.74582704784665232	
367943 Duende	2456158.5	0.10812843436824709	0.893491631652941	

Perihelion epoch and orbital-plane orientation.
Object	Tp TDB	Omega deg	omega deg	i deg	
433 Eros	2453371.5859943051	304.40102733795362	178.665326776373	10.8291838260782	
C/2023 A3	2460581.24084083	21.55947863619833	308.4917712569641	139.1121095087364	
2P Encke	2460239.6495185494	334.1935846036774	187.1342463695676	11.41227811179314	
2I Borisov	2458826.052845906	308.1477292269942	209.1236864378081	44.05264247909138	
99942 Apophis	2459101.0394224627	204.0389272089208	126.65205183685531	3.3367513200667558	
367943 Duende	2456261.554400227	147.2625010881184	271.08746874977783	10.33728752192523	

Non-gravitational parameters; zero omits that component.
Object	A1 / A2 / A3	
433 Eros	0 / 0 / 0	
C/2023 A3	1.4218744044776211e-8 / 1.3486232284695711e-9 / 0	
2P Encke	2.34758107815004e-10 / -6.540690743894436e-13 / 0	
2I Borisov	4.8595882640274848e-8 / 1.8798986909513609e-8 / 2.8355032201284751e-9	
99942 Apophis	5e-13 / -2.9017666371531647e-14 / 0	
367943 Duende	0 / 0 / 0	

Frozen JPL Horizons reference positions.
Object	UTC	JPL RA deg	JPL Dec deg	
433 Eros	2026-07-21 22:00:00	188.352878680	-15.485731924	
433 Eros	2026-07-22 22:00:00	189.033073621	-15.673602566	
433 Eros	2026-07-23 22:00:00	189.712493431	-15.860306533	
C/2023 A3	2026-07-21 22:00:00	270.131469713	20.765253967	
C/2023 A3	2026-07-22 22:00:00	269.997011490	20.710979861	
C/2023 A3	2026-07-23 22:00:00	269.864375552	20.655473422	
2P Encke	2026-07-21 22:00:00	22.086326480	16.151079557	
2P Encke	2026-07-22 22:00:00	22.204227929	16.258348832	
2P Encke	2026-07-23 22:00:00	22.318311808	16.365096989	
2I Borisov	2019-12-08 00:00:00	171.852537543	-17.895950115	
2I Borisov	2019-12-09 00:00:00	172.269466247	-18.665188729	
2I Borisov	2019-12-10 00:00:00	172.685173642	-19.435471809	
99942 Apophis	2029-04-13 20:30:00	162.187610650	-1.853259832	
99942 Apophis	2029-04-13 20:40:00	158.279644424	0.774049246	
99942 Apophis	2029-04-13 20:50:00	153.882081064	3.710242552	
99942 Apophis	2029-04-13 21:00:00	148.942949743	6.944475706	
367943 Duende	2013-02-15 19:15:00	180.339825182	-21.468034850	
367943 Duende	2013-02-15 19:20:00	180.569814701	-17.847154119	
367943 Duende	2013-02-15 19:25:00	180.799653390	-14.126076156	
367943 Duende	2013-02-15 19:30:00	181.029779486	-10.326525675	

Calculated positions and spherical separation from the reference.
Object	UTC	Core RA deg	Core Dec deg	Delta arcsec	
433 Eros	2026-07-21 22:00:00	188.352878496	-15.485732361	0.001696	
433 Eros	2026-07-22 22:00:00	189.033073460	-15.673603008	0.001687	
433 Eros	2026-07-23 22:00:00	189.712493295	-15.860306981	0.001680	
C/2023 A3	2026-07-21 22:00:00	270.131471956	20.765254400	0.007710	
C/2023 A3	2026-07-22 22:00:00	269.997013729	20.710980297	0.007701	
C/2023 A3	2026-07-23 22:00:00	269.864377787	20.655473860	0.007693	
2P Encke	2026-07-21 22:00:00	22.086326254	16.151079161	0.001626	
2P Encke	2026-07-22 22:00:00	22.204227701	16.258348433	0.001638	
2P Encke	2026-07-23 22:00:00	22.318311578	16.365096587	0.001650	
2I Borisov	2019-12-08 00:00:00	171.852535053	-17.895966055	0.058015	
2I Borisov	2019-12-09 00:00:00	172.269463755	-18.665203001	0.052076	
2I Borisov	2019-12-10 00:00:00	172.685171171	-19.435484512	0.046495	
99942 Apophis	2029-04-13 20:30:00	162.187667572	-1.853286529	0.226242	
99942 Apophis	2029-04-13 20:40:00	158.279702838	0.774021995	0.232031	
99942 Apophis	2029-04-13 20:50:00	153.882138225	3.710216547	0.225683	
99942 Apophis	2029-04-13 21:00:00	148.943005706	6.944451473	0.218189	
367943 Duende	2013-02-15 19:15:00	180.339837945	-21.468014339	0.085326	
367943 Duende	2013-02-15 19:20:00	180.569829577	-17.847141796	0.067577	
367943 Duende	2013-02-15 19:25:00	180.799668504	-14.126072752	0.054168	
367943 Duende	2013-02-15 19:30:00	181.029797039	-10.326531491	0.065598	

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 1850-01-01 through 2150-01-01; 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.
Language	Max. Delta arcsec	Result	Source B	Package B	Toolchain	
Ada	0.223931	PASS	66,911	288,626	GNAT/GCC 15.2.0 (MSYS2)	
B4X	0.232382	PASS	65,725	299,027	B4J 10.50	
Ballerina	0.238317	PASS	55,513	298,119	Ballerina Swan Lake 2201.13.4	
BlitzMax	1.203843	PASS	58,926	295,334	BlitzMax NG 0.154.3.58	
C	0.232031	PASS	56,797	315,498	GCC 16.1.0	
C#	0.232031	PASS	50,914	301,338	.NET SDK 8.0.423	
C++	0.238317	PASS	57,599	305,519	Microsoft Visual C++ Community	
Clojure	0.235614	PASS	45,664	266,025	Clojure 1.12.5 / OpenJDK	
Common Lisp	0.223203	PASS	45,788	271,099	SBCL 2.6.6	
D	0.250215	PASS	38,779	283,920	DMD 2.112.0	
Dart	0.228185	PASS	43,237	277,530	Dart 3.12.2	
Delphi/Object Pascal	0.232316	PASS	49,102	286,420	Free Pascal 3.2.2	
Elixir	0.243787	PASS	53,836	314,208	Elixir 1.20.2 / Erlang OTP 29	
Erlang	0.243787	PASS	53,038	318,577	Erlang/OTP 29.0.3	
F#	0.216606	PASS	122,496	409,990	.NET SDK 8.0.423	
Factor	51.049773	PASS	43,875	269,748	Factor 0.101	
Fantom	0.238317	PASS	55,416	300,877	Fantom 1.0.83	
Fortran	0.226948	PASS	51,600	272,368	GNU Fortran 16.1.0	
FreeBASIC	0.252308	PASS	68,277	316,719	FreeBASIC 1.10.1	
GNU Octave	38.034360	PASS	36,553	276,734	GNU Octave 11.3.0	
Go	0.239381	PASS	40,236	289,032	Go 1.26.5	
Groovy	2.810988	PASS	61,373	302,461	Apache Groovy 5.0.2 / OpenJDK	
Haskell	0.259614	PASS	45,149	267,224	GHC 9.10.3	
Haxe	0.228185	PASS	47,482	276,893	Haxe 4.3.7	
Java	0.238317	PASS	47,394	300,055	Microsoft OpenJDK 25.0.3	
JavaScript	0.236219	PASS	108,905	324,762	Node.js 24.15.0 / Browser ES Modules	
Julia	0.229764	PASS	31,521	259,975	Julia 1.12.6	
Kotlin	0.238317	PASS	40,456	264,701	Kotlin 2.4.10 / OpenJDK	
Lua	0.216606	PASS	42,288	285,772	Lua 5.4	
Nim	0.228185	PASS	38,587	285,635	Nim 2.2.10	
OCaml	0.232382	PASS	114,349	391,270	OCaml 5.4.0	
Perl	0.236719	PASS	43,670	256,197	Perl 5	
PHP	0.232031	PASS	33,586	281,611	PHP 8.3.30 CLI	
PowerShell	37.997961	PASS	50,536	261,593	PowerShell 7.6.4	
Prolog	38.032226	PASS	56,296	340,512	SWI-Prolog 10.0.2	
Python	0.082129	PASS	219,557	435,161	Python 3.11.9	
R	0.231526	PASS	34,543	275,740	R 4.6.1	
Racket	0.218560	PASS	45,215	264,552	Racket 9.2	
Ruby	0.238317	PASS	37,728	278,057	Ruby 3.4	
Rust	0.228185	PASS	50,204	297,171	rustc 1.94.1	
Scala	0.238317	PASS	39,182	253,170	Scala CLI 1.15.0 / OpenJDK	
Smalltalk	42.792737	PASS	35,322	228,707	Pharo 12	
Squirrel	0.228185	PASS	45,488	289,555	Squirrel 3 with float64 build	
Tcl	0.232382	PASS	48,183	256,273	Tcl 8.6	
TypeScript + WebAssembly	0.239511	PASS	57,437	310,839	TypeScript 7.0.2, Rust WASI, Node.js 24.15.0	
V	0.236584	PASS	40,910	278,165	V 0.4	
Vala/Genie	0.228185	PASS	62,755	315,396	Vala 0.56	
Visual Basic .NET	0.232031	PASS	105,633	350,275	.NET SDK 8.0.423	
Zig	0.239824	PASS	45,124	286,945	Zig 0.16.0	

A.2 Performance, fastest first
Factor 1.000 denotes the fastest median in the validation run.
Rank	Language	Minimum ms	Median ms	Maximum ms	Factor	
1	TypeScript + WebAssembly	597.377	700.999	1047.449	1.000	
2	Fortran	596.688	750.360	1158.344	1.070	
3	Rust	523.705	823.898	1294.582	1.175	
4	B4X	775.000	971.000	2194.000	1.385	
5	Haxe	39.747	1077.578	2013.806	1.537	
6	Delphi/Object Pascal	937.000	1078.000	1672.000	1.538	
7	Zig	1035.454	1160.447	2320.055	1.655	
8	C++	1223.592	1441.870	2671.024	2.057	
9	Java	1284.711	1458.187	3353.160	2.080	
10	C	1343.938	1486.065	3056.110	2.120	
11	Kotlin	1176.612	1593.443	5530.074	2.273	
12	Ada	1460.014	1640.974	3057.029	2.341	
13	V	1575.736	1878.741	2571.347	2.680	
14	D	1546.526	1890.757	4220.870	2.697	
15	Go	1468.523	2141.171	9256.961	3.054	
16	Visual Basic .NET	1696.939	2216.958	3742.810	3.163	
17	F#	1783.725	2231.113	3564.241	3.183	
18	Nim	2119.599	2326.564	5060.611	3.319	
19	C#	2374.419	2552.401	6731.288	3.641	
20	Dart	2190.595	2733.704	5027.791	3.900	
21	Scala	2073.900	3242.849	7436.041	4.626	
22	OCaml	2669.212	3467.746	8892.401	4.947	
23	FreeBASIC	2767.911	3570.477	5810.864	5.093	
24	Julia	3434.719	4304.078	8335.834	6.140	
25	Fantom	4241.671	5794.304	8545.390	8.266	
26	Vala/Genie	4412.042	5820.308	8666.353	8.303	
27	Haskell	6047.352	8690.853	14693.710	12.398	
28	BlitzMax	7899.000	8742.000	12931.000	12.471	
29	JavaScript	7724.900	10157.720	23217.210	14.490	
30	Common Lisp	11056.269	14392.095	22933.360	20.531	
31	Ballerina	12341.616	16894.464	33663.612	24.101	
32	Racket	16694.584	21229.729	44060.672	30.285	
33	Erlang	19354.829	24802.406	45274.726	35.381	
34	Elixir	23434.035	29699.584	53299.712	42.367	
35	PHP	41451.985	49798.621	104949.488	71.039	
36	Prolog	49043.000	69317.000	255244.000	98.883	
37	Lua	53249.000	69977.000	130427.000	99.825	
38	Squirrel	55843.000	71891.000	119351.000	102.555	
39	Clojure	55241.707	72206.563	171065.047	103.005	
40	Ruby	66049.112	90263.925	168102.051	128.765	
41	Factor	70004.415	94462.332	151220.285	134.754	
42	Smalltalk	65113.000	95089.000	157550.000	135.648	
43	Groovy	68129.578	95282.139	143061.528	135.923	
44	Python	74058.281	138890.584	308314.450	198.132	
45	Perl	108115.457	154930.755	276448.215	221.014	
46	PowerShell	206315.845	250942.054	717706.460	357.978	
47	Tcl	202078.602	265143.768	494701.078	378.237	
48	R	168780.000	343730.000	910730.000	490.343	
49	GNU Octave	315324.607	419255.178	922894.954	598.082	

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/