Repository navigation
Algorithm Reference
This page documents the algorithms, formulae, and references used by the Ephemeris library. Each section maps to a namespace in the core library and cites the primary sources.
Primary reference: Jean Meeus, Astronomical Algorithms, 2nd ed. (Willmann-Bell, 1998). Cited as "Meeus Ch. N".
- Timekeeping (Chronology)
- Solar Ephemeris (Heliology)
- Lunar Ephemeris (Selenography)
- Topocentric Parallax
- Planetary Positions (Planetology)
- Coordinate Transforms (Geometry)
- Nutation & Precession (Geodesy)
- Observable Phenomena (Phenomenology)
- Fixed Stars (Stellarography)
- SPICE/BSP Import (Import)
Class: TimeUtils, TimeZoneUtils
Source: Meeus Ch. 7 (Julian Day), Ch. 10 (ΔT), Ch. 12 (GMST)
JD = 365.25 × (Y + 4716) + 30.6001 × (M + 1) + D + B − 1524.5
where B is the Gregorian calendar correction. Dates before 15 Oct 1582 use the Julian calendar (B = 0).
Julian Century:
T = (JD − 2451545.0) / 36525.0
J2000.0 epoch = JD 2451545.0 = 2000 January 1.5 TT.
Polynomial approximations by era from Morrison & Stephenson (2004) and the IERS. Five time-range branches are used; the post-2005 branch is:
ΔT ≈ 62.92 + 0.32217(y − 2000) + 0.005589(y − 2000)² (seconds)
Meeus Eq. 12.4 (degrees, normalized to [0, 360)):
GMST = 280.46061837 + 360.98564736629 × (JD − J2000)
+ 0.000387933 × T² − T³ / 38710000
Class: SunEphemeris
Source: Meeus Ch. 25 (low-precision solar coordinates), Ch. 22 (aberration), Ch. 22 (nutation)
Accuracy: ~0.01°
L₀ = 280.46646 + 36000.76983 T + 0.0003032 T² (degrees)
M = 357.52911 + 35999.05029 T − 0.0001537 T² (degrees)
C = (1.914602 − 0.004817T − 0.000014T²) sin M
+ (0.019993 − 0.000101T) sin 2M
+ 0.000289 sin 3M
Sun's true longitude: Θ = L₀ + C
Apparent longitude with aberration: λ = Θ − 0.00569 − 0.00478 sin Ω
where Ω = Moon's ascending node longitude (Meeus Eq. 25.9).
ε₀ = 23° 26′ 21.448″ − 4680.93″T − 1.55″T² + 1999.25″T³ − …
ε = ε₀ + 0.00256 cos Ω (apparent obliquity including nutation)
Class: MoonEphemeris
Source: Meeus Ch. 47 (ELP-2000/82 truncated series)
Accuracy: geocentric ~0.1°; topocentric after parallax correction ~0.01° additional error
| Symbol | Meaning | Meeus Eq. |
|---|---|---|
| L′ | Moon's mean longitude | 47.1 |
| D | Moon's mean elongation | 47.2 |
| M | Sun's mean anomaly | 47.3 |
| M′ | Moon's mean anomaly | 47.4 |
| F | Moon's argument of latitude | 47.5 |
Longitude correction Σl: 60-term series summing A_i sin(arg_i), where each argument is a linear combination of D, M, M′, F.
Latitude correction Σb: 60-term series summing B_i sin(arg_i).
Distance correction Σr: 25-term series summing C_i cos(arg_i).
Eccentricity factor e = 1 − 0.002516T − 0.0000074T² modifies terms involving M.
Class: TopocentricParallax
Source: Meeus Ch. 40
Accuracy: ~0.01° additional correction to geocentric position
For the Moon: sin π = 6378.14 km / Δ where Δ is geocentric distance.
For the Sun: π☉ ≈ 8.794″.
For planets: π = 8.794″ × (1 AU / Δ).
Observer's reduced latitude:
ρ sin φ′ = 0.99664719 sin φ + (h/6378140) sin φ
ρ cos φ′ = cos φ + (h/6378140) cos φ
ΔRA (Meeus Eq. 40.6):
ΔRA = atan[ −ρ cos φ′ sin π sin H / (cos δ − ρ cos φ′ sin π cos H) ]
Topocentric Dec (Meeus Eq. 40.7):
δ′ = atan[ (sin δ − ρ sin φ′ sin π) cos ΔRA / (cos δ − ρ cos φ′ sin π cos H) ]
Class: PlanetEphemeris, PlanetPhysicalEphemeris
Source: Meeus Ch. 33 (simplified orbital elements), Ch. 41 (magnitudes), Ch. 26 (physical ephemeris)
Accuracy: 0.5°–5° depending on planet and epoch
Each planet is represented by six osculating elements at J2000.0 with linear drift in T:
| Element | Symbol |
|---|---|
| Longitude of ascending node | Ω |
| Inclination | i |
| Argument of perihelion | ω |
| Semi-major axis | a (AU) |
| Eccentricity | e |
| Mean anomaly | M |
The eccentric anomaly E satisfies E − e sin E = M. Solved iteratively:
E₀ = M
E_{n+1} = E_n + (M − E_n + e sin E_n) / (1 − e cos E_n)
Convergence in 10–15 iterations for e < 0.9.
xh = r [cos Ω cos(ω+ν) − sin Ω sin(ω+ν) cos i]
yh = r [sin Ω cos(ω+ν) + cos Ω sin(ω+ν) cos i]
zh = r sin(ω+ν) sin i
Subtract Earth's heliocentric position (from SunEphemeris.HeliocentricLongitude) to obtain geocentric ecliptic XYZ, then convert to RA/Dec.
Classes: ObserverGeometry, CoordinateConverter
Source: Meeus Ch. 13 (equatorial ↔ horizontal), Ch. 93 (atmospheric refraction)
H = GMST + λ − RA (local hour angle, degrees)
Alt = arcsin(sin φ sin δ + cos φ cos δ cos H)
Az = atan2(sin H, cos H sin φ − tan δ cos φ)
Az = Az + 180° if sin H > 0 (quadrant correction)
For apparent altitude h_app in degrees:
R = 1.02 / tan(h_app + 10.3 / (h_app + 5.11)) (arcminutes)
Cutoff: not applied below h_app = −1°.
Inverse (Saemundsson): h_true = h_app − R(h_app)
sin δ = sin ε sin λ + cos ε cos λ sin β … (full transform via rotation by ε)
Classes: NutationCalculator, PrecessionCalculator
Source: Meeus Ch. 22; IAU 1980 nutation theory
48 of the 106 standard terms retained (covering > 99.9% of amplitude).
Nutation in longitude Δψ (arcseconds):
Δψ = Σ (S_i + S′_i T) sin(arg_i)
Nutation in obliquity Δε (arcseconds):
Δε = Σ (C_i + C′_i T) cos(arg_i)
Each argument is a linear combination of: D, M☉, M☽, F, Ω.
Source: IAU 2006 precession model; Meeus Ch. 21
Accumulated precession angles from J2000.0:
ψ_A = 5038.481507″T − 1.0790069″T² − 0.00114045″T³ + …
ω_A = ε₀ − 0.025754″T + 0.0512623″T² − 0.00772503″T³ + …
χ_A = 10.556403″T − 2.3814292″T² − 0.00121197″T³ + …
Class: RiseSetCalculator
-
Compute approximate HA at rise/set:
cos H₀ = (sin h₀ − sin φ sin δ) / (cos φ cos δ)where h₀ is the standard altitude (−0.8333° Sun, −0.5667° stars).
-
Estimate fractional day:
m_transit = (RA − λ − θ₀) / 360 -
Three-iteration correction using Meeus Eq. 15.1–15.3 with three-point interpolation for RA/Dec and ΔT correction.
Class: EclipseCalculator
Lunation index k such that k integer = new moon, k + 0.5 = full moon. Julian centuries T = k / 1236.85.
Quick filter: if |sin F| > 0.36 (Moon far from node), no eclipse possible.
Eclipse parameters:
- gamma: shadow axis distance from Earth's centre in Earth radii
- u: penumbral cone parameter
Classification thresholds (Meeus Table 54.a):
| Condition | Type |
|---|---|
| gamma | |
| gamma | |
| gamma | |
| 0.9972 < | gamma |
| sin F₁ | |
| 0.9972 < | sin F₁ |
Class: SeasonCalculator
JDE of mean March equinox:
JDE₀ = 2451623.80984 + 365242.37404T + 0.05169T² − 0.00411T³ − 0.00057T⁴
(analogous polynomials for June solstice, September equinox, December solstice)
Corrections for solar perturbations via 24-coefficient series in W = 2π(JDE₀ − 2451545) / 365.25.
Class: PlanetaryEventCalculator, InnerPlanetEventCalculator
Signed elongation ε ∈ (−180°, +180°] (positive = east of Sun):
ε = normalize(λ_planet − λ_sun, −180, +180)
Events detected by scanning in 0.5-day steps and detecting the target sign change:
| Event | Detection |
|---|---|
| Opposition | ε wraps from ≈ −180 to ≈ +180 |
| Conjunction | ε crosses 0 (pos → neg) |
| East quadrature | ε decreases through +90° |
| West quadrature | ε decreases through −90° |
Sub-step interpolation gives ~hours precision. Greatest elongation for inner planets is found by a golden-section maximisation of |ε| over the bracketing interval.
Classes: StarEphemeris, BrightStarCatalog, StarCatalog
Source: Meeus Ch. 21 (precession), Ch. 21 (proper motion)
Position at epoch JD from J2000.0 catalogue position:
RA(JD) = RA₀ + μ_α cos δ × Δt (μ in arcsec/yr, Δt in Julian years)
Dec(JD) = Dec₀ + μ_δ × Δt
Rigorous precession matrix using IAU 2006 angles ψ_A, ω_A, χ_A applied to J2000.0 ICRS unit vector.
Classes: SpkReader, SpiceKernelDatabase, BspImporter
Reference: NAIF DAF/SPK format
Binary SPK kernels use the Double Array File (DAF) format:
- 1024-byte file record (ASCII + binary header)
- Linked list of summary records, each containing segment descriptors
- Each descriptor identifies: target body, centre body, frame, data type, start/end ET
Segment data is divided into equal-length records. Each record contains:
- Start epoch (ET seconds), interval length, N Chebyshev coefficients per axis (X, Y, Z)
Evaluation at time t:
t_norm = 2(t − t_mid) / interval ∈ [−1, 1]
pos[axis] = Σ c_k T_k(t_norm) (Clenshaw recurrence)
ET = (JD_UTC − J2000_JD) × 86400 + ΔAT + 32.184
where ΔAT is the number of leap seconds since 1972. The leap second table is embedded in SpkLeapSeconds.
See also: SPK-BSP-Format, SE1-Ephemeris-Format, SEFStars-Catalog-Format, Yale-BSC5-Format