\documentclass[10pt,a4paper]{article}
\usepackage[english]{babel}
\usepackage{fontspec}
\setmainfont{Source Serif Pro}
\setsansfont{Source Sans Pro}
\setmonofont{Source Code Pro}
\usepackage{microtype}
\usepackage[a4paper,inner=23mm,outer=19mm,top=21mm,bottom=24mm]{geometry}
\usepackage{amsmath,amssymb}
\usepackage{booktabs,longtable,array}
\usepackage[table]{xcolor}
\usepackage{hyperref}
\usepackage{fancyhdr}
\usepackage{enumitem}
\usepackage{fvextra}
\definecolor{ephblue}{HTML}{1F4D78}
\definecolor{ephstripe}{HTML}{F4F8FB}
\hypersetup{
  colorlinks=true,
  linkcolor=ephblue,
  urlcolor=ephblue,
  pdftitle={Numerical Ephemeris Calculation},
  pdfauthor={Jost Jahn},
  pdfsubject={Scientific and technical description of the independent Ephemeris numerical core}
}
\pagestyle{fancy}
\fancyhf{}
\fancyhead[LE,RO]{\small\color{ephblue} Numerical Ephemeris Calculation}
\fancyhead[RE,LO]{\small\color{gray} Amrum Software}
\fancyfoot[C]{\thepage}
\setlength{\headheight}{14pt}
\setlist{nosep,leftmargin=5mm}
\renewcommand{\arraystretch}{1.18}
\setlength{\tabcolsep}{4pt}
\arrayrulecolor{ephblue!35}
\newcommand{\thead}[1]{\color{white}\bfseries #1}
\title{\color{ephblue}\bfseries Numerical Ephemeris Calculation\\[2mm]
\large Mathematical method, constants, series and cross-language validation}
\author{Jost Jahn}
\date{July 25, 2026 \quad Version 0.5.3}

\begin{document}
\maketitle
\tableofcontents
\newpage

\section*{Abstract}
\addcontentsline{toc}{section}{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.

\section{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.


\noindent\textit{Glossary of astronomical and numerical terms.}\par\smallskip
\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.2667\linewidth}p{0.6933\linewidth}@{}}
\rowcolor{ephblue}\thead{Term} & \thead{Meaning in this article}\\
\endfirsthead
\rowcolor{ephblue}\thead{Term} & \thead{Meaning in this article}\\
\endhead
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. \\
\end{longtable}


\noindent\textit{Solar and terrestrial gravitational parameters are kept unambiguous.}\par\smallskip
\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.2051\linewidth}p{0.4677\linewidth}p{0.2872\linewidth}@{}}
\rowcolor{ephblue}\thead{Symbol} & \thead{Meaning} & \thead{Value or unit}\\
\endfirsthead
\rowcolor{ephblue}\thead{Symbol} & \thead{Meaning} & \thead{Value or unit}\\
\endhead
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 \\
\end{longtable}

\section{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.



\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.2564\linewidth}p{0.4410\linewidth}p{0.2626\linewidth}@{}}
\rowcolor{ephblue}\thead{Quantity} & \thead{Exact value} & \thead{Unit}\\
\endfirsthead
\rowcolor{ephblue}\thead{Quantity} & \thead{Exact value} & \thead{Unit}\\
\endhead
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 \\
\end{longtable}

\section{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.

\subsection{Ellipse}
\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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))
\end{Verbatim}
\subsection{Hyperbola}
\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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))
\end{Verbatim}
\subsection{Parabola and the near-parabolic boundary}
\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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
\end{Verbatim}
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:

\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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]
\end{Verbatim}
Rotation from the orbital plane uses Rz(Omega) Rx(i) Rz(omega). For a two-dimensional input [x,y], the expanded form is:

\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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
\end{Verbatim}
A fixed x-axis rotation by 84381.448 arcseconds then transforms J2000 ecliptic position and velocity to the equatorial frame.

\section{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:

\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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
\end{Verbatim}


\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.4410\linewidth}p{0.5190\linewidth}@{}}
\rowcolor{ephblue}\thead{Perturbing body} & \thead{Mass/Sun mass}\\
\endfirsthead
\rowcolor{ephblue}\thead{Perturbing body} & \thead{Mass/Sun mass}\\
\endhead
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 \\
\end{longtable}



\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.4410\linewidth}p{0.5190\linewidth}@{}}
\rowcolor{ephblue}\thead{N16 body} & \thead{Mass/Sun mass}\\
\endfirsthead
\rowcolor{ephblue}\thead{N16 body} & \thead{Mass/Sun mass}\\
\endhead
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 \\
\end{longtable}

\section{Equation of motion}
\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
d/dt [r, v] = [v, a(r,v,t)]
\end{Verbatim}
The solar base acceleration is a\_s=-mu\_sun r/|r|\textasciicircum{}3. A perturbing body at heliocentric position r\_b and with mass ratio eta\_b contributes the direct and indirect terms together:

\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
a_b = mu_sun eta_b [(r_b - r)/|r_b - r|^3 - r_b/|r_b|^3]
\end{Verbatim}
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.

\subsection{Relativistic solar correction}
\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
a_rel = mu_sun / (c^2 |r|^3)
        * [(4 mu_sun/|r| - |v|^2) r + 4 (r dot v) v]
\end{Verbatim}
\subsection{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:

\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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)
\end{Verbatim}
\subsection{Non-gravitational acceleration}
\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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]
\end{Verbatim}
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.

\section{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:

\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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]
\end{Verbatim}
For every state component, scale=2e-15+2e-13*max(|y\_old|,|y\_new|). The normalized RMS error is sqrt(sum((error/scale)\textasciicircum{}2)/6). A step is accepted for error <=1 or |h|<=1e-8 day. The next factor is clamp(0.9*error\textasciicircum{}(-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.

\section{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:

\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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
\end{Verbatim}
With T=(JD\_UTC-2451545)/36525, GMST and the complete short nutation and precession series used by the port are:

\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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
\end{Verbatim}
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.

\section{Light time and apparent direction}
Earth and observer states remain at reception time. Emission time is iterated at most five times:

\begin{Verbatim}[breaklines=true,breakanywhere=true,fontsize=\small]
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))
\end{Verbatim}
\section{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.



\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.1744\linewidth}p{0.1846\linewidth}p{0.2256\linewidth}p{0.3754\linewidth}@{}}
\rowcolor{ephblue}\thead{Dataset} & \thead{Interval} & \thead{Step} & \thead{Content}\\
\endfirsthead
\rowcolor{ephblue}\thead{Dataset} & \thead{Interval} & \thead{Step} & \thead{Content}\\
\endhead
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 \\
\end{longtable}

\section{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.


\noindent\textit{Orbit form and perihelion distance of the six validation cases.}\par\smallskip
\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.1846\linewidth}p{0.2256\linewidth}p{0.2933\linewidth}p{0.2564\linewidth}@{}}
\rowcolor{ephblue}\thead{Object} & \thead{Epoch JD} & \thead{e} & \thead{q AU}\\
\endfirsthead
\rowcolor{ephblue}\thead{Object} & \thead{Epoch JD} & \thead{e} & \thead{q AU}\\
\endhead
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 \\
\end{longtable}


\noindent\textit{Perihelion epoch and orbital-plane orientation.}\par\smallskip
\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.1846\linewidth}p{0.2523\linewidth}p{0.1744\linewidth}p{0.1744\linewidth}p{0.1744\linewidth}@{}}
\rowcolor{ephblue}\thead{Object} & \thead{Tp TDB} & \thead{Omega deg} & \thead{omega deg} & \thead{i deg}\\
\endfirsthead
\rowcolor{ephblue}\thead{Object} & \thead{Tp TDB} & \thead{Omega deg} & \thead{omega deg} & \thead{i deg}\\
\endhead
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 \\
\end{longtable}


\noindent\textit{Non-gravitational parameters; zero omits that component.}\par\smallskip
\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.2051\linewidth}p{0.7549\linewidth}@{}}
\rowcolor{ephblue}\thead{Object} & \thead{A1 / A2 / A3}\\
\endfirsthead
\rowcolor{ephblue}\thead{Object} & \thead{A1 / A2 / A3}\\
\endhead
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 \\
\end{longtable}


\noindent\textit{Frozen JPL Horizons reference positions.}\par\smallskip
\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.1949\linewidth}p{0.2564\linewidth}p{0.2544\linewidth}p{0.2544\linewidth}@{}}
\rowcolor{ephblue}\thead{Object} & \thead{UTC} & \thead{JPL RA deg} & \thead{JPL Dec deg}\\
\endfirsthead
\rowcolor{ephblue}\thead{Object} & \thead{UTC} & \thead{JPL RA deg} & \thead{JPL Dec deg}\\
\endhead
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 \\
\end{longtable}


\noindent\textit{Calculated positions and spherical separation from the reference.}\par\smallskip
\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.1744\linewidth}p{0.2359\linewidth}p{0.1949\linewidth}p{0.1949\linewidth}p{0.1600\linewidth}@{}}
\rowcolor{ephblue}\thead{Object} & \thead{UTC} & \thead{Core RA deg} & \thead{Core Dec deg} & \thead{Delta arcsec}\\
\endfirsthead
\rowcolor{ephblue}\thead{Object} & \thead{UTC} & \thead{Core RA deg} & \thead{Core Dec deg} & \thead{Delta arcsec}\\
\endhead
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 \\
\end{longtable}

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.

\section{Reproducibility and limitations}
\begin{itemize}
\item The one-arcsecond limit applies to the same frozen nominal elements, not to a later orbit solution.
\item The short nutation, precession and TDB approximation targets the documented interval and accuracy.
\item Leap seconds after the last built-in step must be added before later epochs are used.
\item The tables cover 1850-01-01 through 2150-01-01; other intervals require a new dataset.
\item Uncertainty ellipses and covariance propagation are outside this nominal core.
\item Find\_Orb is neither runtime component nor fallback; it is used only for private local comparison.
\end{itemize}
\appendix
\section{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.

\subsection{Accuracy, alphabetical}

\noindent\textit{Maximum spherical angular separation over all 20 reference rows.}\par\smallskip
\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.1538\linewidth}p{0.1231\linewidth}p{0.0974\linewidth}p{0.1231\linewidth}p{0.1231\linewidth}p{0.3395\linewidth}@{}}
\rowcolor{ephblue}\thead{Language} & \thead{Max. Delta arcsec} & \thead{Result} & \thead{Source B} & \thead{Package B} & \thead{Toolchain}\\
\endfirsthead
\rowcolor{ephblue}\thead{Language} & \thead{Max. Delta arcsec} & \thead{Result} & \thead{Source B} & \thead{Package B} & \thead{Toolchain}\\
\endhead
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 \\
\end{longtable}

\subsection{Performance, fastest first}

\noindent\textit{Factor 1.000 denotes the fastest median in the validation run.}\par\smallskip
\rowcolors{2}{ephstripe}{white}
\begin{longtable}{@{}p{0.0718\linewidth}p{0.2010\linewidth}p{0.1641\linewidth}p{0.1641\linewidth}p{0.1641\linewidth}p{0.1949\linewidth}@{}}
\rowcolor{ephblue}\thead{Rank} & \thead{Language} & \thead{Minimum ms} & \thead{Median ms} & \thead{Maximum ms} & \thead{Factor}\\
\endfirsthead
\rowcolor{ephblue}\thead{Rank} & \thead{Language} & \thead{Minimum ms} & \thead{Median ms} & \thead{Maximum ms} & \thead{Factor}\\
\endhead
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 \\
\end{longtable}

\section{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: \url{https://www.jostjahn.de/amrum-software/ephemeris/}


\end{document}
