Something went wrong. Try again.
A Zig library for astronomical calculations: the Sun and Moon from anywhere on Earth, the Moon's phases and libration, and time scales.
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>// SPDX-License-Identifier: MIT
//! The Moon's apparent position (Meeus, *Astronomical Algorithms*, 2nd ed.,//! ch. 47): the principal terms of Chapront's ELP-2000/82 lunar theory, with//! his later mean arguments.
const std = @import("std");const math = std.math;
const Angle = @import("../Angle.zig");const coord = @import("../coord.zig");const earth = @import("../earth.zig");const Instant = @import("../time.zig").Instant;const observer_ = @import("../observer.zig");const Observer = observer_.Observer;const Topocentric = observer_.Topocentric;
/// Where the Moon appears from the Earth's center at an instant.pub const Position = struct { /// Apparent ecliptic coordinates, referred to the true equinox of date. ecliptic: coord.Ecliptic, /// Apparent equatorial coordinates, referred to the true equator and /// equinox of date. equatorial: coord.Equatorial, /// The distance between the centers of the Earth and the Moon, in /// kilometers. distance_km: f64,};
/// The Moon's apparent position at an instant, in either time scale.////// Meeus keeps the 60 largest terms of the theory in longitude and distance/// and the 60 largest in latitude, which he gives as good to about 10″ in/// longitude and 4″ in latitude. The theory was fitted to the centuries/// around 1900–2100, and its mean arguments are polynomials in time, so its/// error grows with distance from them: some tens of seconds of arc within a/// thousand years, and degrees beyond a few thousand.////// The mean longitude already includes the constant part of the light-time/// correction, and so the apparent longitude is the geometric one plus the/// nutation in longitude, as Meeus has it (p. 343). The equatorial/// coordinates use the true obliquity.////// NaN if the instant is not finite.pub fn position(instant: Instant) Position { const t = instant.julianCenturies(); const a = Arguments.at(t);
var sum_l: f64 = 0; var sum_r: f64 = 0; for (table_47a) |term| { const argument = a.combine(term.d, term.m, term.m_moon, term.f); const e = a.eccentricityFactor(term.m); sum_l += term.l * e * @sin(argument); sum_r += term.r * e * @cos(argument); } var sum_b: f64 = 0; for (table_47b) |term| { const argument = a.combine(term.d, term.m, term.m_moon, term.f); sum_b += term.b * a.eccentricityFactor(term.m) * @sin(argument); }
// The additive terms: Venus (A₁), Jupiter (A₂), and the flattening of // the Earth (L′) (p. 342). const l_moon = a.l_moon * math.rad_per_deg; const m_moon = a.m_moon * math.rad_per_deg; const f = a.f * math.rad_per_deg; const a1 = (119.75 + 131.849 * t) * math.rad_per_deg; const a2 = (53.09 + 479264.290 * t) * math.rad_per_deg; const a3 = (313.45 + 481266.484 * t) * math.rad_per_deg; sum_l += 3958 * @sin(a1) + 1962 * @sin(l_moon - f) + 318 * @sin(a2); sum_b += -2235 * @sin(l_moon) + 382 * @sin(a3) + 175 * @sin(a1 - f) + 175 * @sin(a1 + f) + 127 * @sin(l_moon - m_moon) - 115 * @sin(l_moon + m_moon);
const nutation = earth.nutation(instant); const ecliptic: coord.Ecliptic = .{ .longitude = Angle.fromDegrees(a.l_moon + sum_l / 1e6).add(nutation.longitude).normalized(), .latitude = Angle.fromDegrees(sum_b / 1e6), }; const obliquity = earth.meanObliquity(instant).add(nutation.obliquity); return .{ .ecliptic = ecliptic, .equatorial = ecliptic.toEquatorial(obliquity), .distance_km = 385000.56 + sum_r / 1000, };}
/// The Moon's apparent topocentric position at an instant, in either time/// scale: its apparent position passed through `astro.observer.topocentric`./// The parallax reaches about 1°, so an observer sees the Moon up to that/// far from where `position` puts it, and lower in the sky. NaN if the/// instant is not finite or the observer is not valid.pub fn topocentric(instant: Instant, observer: Observer) Topocentric { const p = position(instant); return observer_.topocentric(instant, observer, p.equatorial, p.distance_km);}
/// The theory's mean arguments at T, in degrees, unreduced (eqs./// 47.1–47.6). Shared with the calculations of the Moon's phase and/// libration that are built on the same theory.pub const Arguments = struct { t: f64, /// The Moon's mean longitude, L′, including the constant term of the /// light-time (−0.70″). l_moon: f64, /// The mean elongation of the Moon, D. d: f64, /// The Sun's mean anomaly, M. m: f64, /// The Moon's mean anomaly, M′. m_moon: f64, /// The Moon's argument of latitude, F. f: f64, /// The factor for the decreasing eccentricity of the Earth's orbit, E. e: f64,
pub fn at(t: f64) Arguments { return .{ .t = t, .l_moon = polynomial(t, .{ 218.3164477, 481267.88123421, -0.0015786, 1.0 / 538841.0, -1.0 / 65194000.0 }), .d = polynomial(t, .{ 297.8501921, 445267.1114034, -0.0018819, 1.0 / 545868.0, -1.0 / 113065000.0 }), .m = polynomial(t, .{ 357.5291092, 35999.0502909, -0.0001536, 1.0 / 24490000.0 }), .m_moon = polynomial(t, .{ 134.9633964, 477198.8675055, 0.0087414, 1.0 / 69699.0, -1.0 / 14712000.0 }), .f = polynomial(t, .{ 93.2720950, 483202.0175233, -0.0036539, -1.0 / 3526000.0, 1.0 / 863310000.0 }), .e = polynomial(t, .{ 1, -0.002516, -0.0000074 }), }; }
/// The argument d·D + m·M + m′·M′ + f·F, in radians. fn combine(self: Arguments, d: i8, m: i8, m_moon: i8, f: i8) f64 { return (@as(f64, @floatFromInt(d)) * self.d + @as(f64, @floatFromInt(m)) * self.m + @as(f64, @floatFromInt(m_moon)) * self.m_moon + @as(f64, @floatFromInt(f)) * self.f) * math.rad_per_deg; }
/// E for a term whose argument contains ±M, E² for ±2M (p. 338). fn eccentricityFactor(self: Arguments, m: i8) f64 { return switch (m) { 0 => 1, -1, 1 => self.e, else => self.e * self.e, }; }};
/// c[0] + c[1]x + c[2]x² + …, by Horner's method.fn polynomial(x: f64, comptime c: anytype) f64 { var sum: f64 = 0; comptime var i = c.len; inline while (i > 0) { i -= 1; sum = sum * x + c[i]; } return sum;}
/// A row of Table 47.A: the multiples of D, M, M′ and F, and the/// coefficients of the sine in Σl (0.000001°) and of the cosine in Σr/// (0.001 km).const LongitudeTerm = struct { d: i8, m: i8, m_moon: i8, f: i8, l: f64, r: f64 };
/// A row of Table 47.B: the multiples of D, M, M′ and F, and the coefficient/// of the sine in Σb (0.000001°).const LatitudeTerm = struct { d: i8, m: i8, m_moon: i8, f: i8, b: f64 };
fn lr(d: i8, m: i8, m_moon: i8, f: i8, l: f64, r: f64) LongitudeTerm { return .{ .d = d, .m = m, .m_moon = m_moon, .f = f, .l = l, .r = r };}
fn lb(d: i8, m: i8, m_moon: i8, f: i8, b: f64) LatitudeTerm { return .{ .d = d, .m = m, .m_moon = m_moon, .f = f, .b = b };}
/// Meeus, Table 47.A (pp. 339–340), in the book's order.const table_47a = [_]LongitudeTerm{ lr(0, 0, 1, 0, 6288774, -20905355), lr(2, 0, -1, 0, 1274027, -3699111), lr(2, 0, 0, 0, 658314, -2955968), lr(0, 0, 2, 0, 213618, -569925), lr(0, 1, 0, 0, -185116, 48888), lr(0, 0, 0, 2, -114332, -3149), lr(2, 0, -2, 0, 58793, 246158), lr(2, -1, -1, 0, 57066, -152138), lr(2, 0, 1, 0, 53322, -170733), lr(2, -1, 0, 0, 45758, -204586), lr(0, 1, -1, 0, -40923, -129620), lr(1, 0, 0, 0, -34720, 108743), lr(0, 1, 1, 0, -30383, 104755), lr(2, 0, 0, -2, 15327, 10321), lr(0, 0, 1, 2, -12528, 0), lr(0, 0, 1, -2, 10980, 79661), lr(4, 0, -1, 0, 10675, -34782), lr(0, 0, 3, 0, 10034, -23210), lr(4, 0, -2, 0, 8548, -21636), lr(2, 1, -1, 0, -7888, 24208), lr(2, 1, 0, 0, -6766, 30824), lr(1, 0, -1, 0, -5163, -8379), lr(1, 1, 0, 0, 4987, -16675), lr(2, -1, 1, 0, 4036, -12831), lr(2, 0, 2, 0, 3994, -10445), lr(4, 0, 0, 0, 3861, -11650), lr(2, 0, -3, 0, 3665, 14403), lr(0, 1, -2, 0, -2689, -7003), lr(2, 0, -1, 2, -2602, 0), lr(2, -1, -2, 0, 2390, 10056), lr(1, 0, 1, 0, -2348, 6322), lr(2, -2, 0, 0, 2236, -9884), lr(0, 1, 2, 0, -2120, 5751), lr(0, 2, 0, 0, -2069, 0), lr(2, -2, -1, 0, 2048, -4950), lr(2, 0, 1, -2, -1773, 4130), lr(2, 0, 0, 2, -1595, 0), lr(4, -1, -1, 0, 1215, -3958), lr(0, 0, 2, 2, -1110, 0), lr(3, 0, -1, 0, -892, 3258), lr(2, 1, 1, 0, -810, 2616), lr(4, -1, -2, 0, 759, -1897), lr(0, 2, -1, 0, -713, -2117), lr(2, 2, -1, 0, -700, 2354), lr(2, 1, -2, 0, 691, 0), lr(2, -1, 0, -2, 596, 0), lr(4, 0, 1, 0, 549, -1423), lr(0, 0, 4, 0, 537, -1117), lr(4, -1, 0, 0, 520, -1571), lr(1, 0, -2, 0, -487, -1739), lr(2, 1, 0, -2, -399, 0), lr(0, 0, 2, -2, -381, -4421), lr(1, 1, 1, 0, 351, 0), lr(3, 0, -2, 0, -340, 0), lr(4, 0, -3, 0, 330, 0), lr(2, -1, 2, 0, 327, 0), lr(0, 2, 1, 0, -323, 1165), lr(1, 1, -1, 0, 299, 0), lr(2, 0, 3, 0, 294, 0), lr(2, 0, -1, -2, 0, 8752),};
/// Meeus, Table 47.B (p. 341): the left half of the page, then the right.const table_47b = [_]LatitudeTerm{ lb(0, 0, 0, 1, 5128122), lb(0, 0, 1, 1, 280602), lb(0, 0, 1, -1, 277693), lb(2, 0, 0, -1, 173237), lb(2, 0, -1, 1, 55413), lb(2, 0, -1, -1, 46271), lb(2, 0, 0, 1, 32573), lb(0, 0, 2, 1, 17198), lb(2, 0, 1, -1, 9266), lb(0, 0, 2, -1, 8822), lb(2, -1, 0, -1, 8216), lb(2, 0, -2, -1, 4324), lb(2, 0, 1, 1, 4200), lb(2, 1, 0, -1, -3359), lb(2, -1, -1, 1, 2463), lb(2, -1, 0, 1, 2211), lb(2, -1, -1, -1, 2065), lb(0, 1, -1, -1, -1870), lb(4, 0, -1, -1, 1828), lb(0, 1, 0, 1, -1794), lb(0, 0, 0, 3, -1749), lb(0, 1, -1, 1, -1565), lb(1, 0, 0, 1, -1491), lb(0, 1, 1, 1, -1475), lb(0, 1, 1, -1, -1410), lb(0, 1, 0, -1, -1344), lb(1, 0, 0, -1, -1335), lb(0, 0, 3, 1, 1107), lb(4, 0, 0, -1, 1021), lb(4, 0, -1, 1, 833), lb(0, 0, 1, -3, 777), lb(4, 0, -2, 1, 671), lb(2, 0, 0, -3, 607), lb(2, 0, 2, -1, 596), lb(2, -1, 1, -1, 491), lb(2, 0, -2, 1, -451), lb(0, 0, 3, -1, 439), lb(2, 0, 2, 1, 422), lb(2, 0, -3, -1, 421), lb(2, 1, -1, 1, -366), lb(2, 1, 0, 1, -351), lb(4, 0, 0, 1, 331), lb(2, -1, 1, 1, 315), lb(2, -2, 0, -1, 302), lb(0, 0, 1, 3, -283), lb(2, 1, 1, -1, -229), lb(1, 1, 0, -1, 223), lb(1, 1, 0, 1, 223), lb(0, 1, -2, -1, -220), lb(2, 1, -1, -1, -220), lb(1, 0, 1, 1, -185), lb(2, -1, -2, -1, 181), lb(0, 1, 2, 1, -177), lb(4, 0, -2, -1, 176), lb(4, -1, -1, -1, 166), lb(1, 0, 1, -1, -164), lb(4, 0, 1, -1, 132), lb(1, 0, -1, -1, -119), lb(4, -1, 0, -1, 115), lb(2, -2, 0, 1, 107),};
comptime { std.debug.assert(table_47a.len == 60); std.debug.assert(table_47b.len == 60);}