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.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>// SPDX-License-Identifier: MIT
//! The Moon's libration and the position angle of its axis, as seen from//! the Earth's center or by an observer on its surface (Meeus,//! *Astronomical Algorithms*, 2nd ed., ch. 53).//!//! Selenographic longitude is positive towards Mare Crisium — selenographic//! east, the IAU's convention and Meeus's — and latitude towards Mare//! Serenitatis. A positive libration in longitude turns the west limb//! towards the Earth; a positive libration in latitude, the north limb.
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 moon_position = @import("position.zig");const Observer = @import("../observer.zig").Observer;const topocentricPosition = @import("../observer.zig").topocentric;const sun = @import("../sun.zig");
/// A point on the Moon's surface: its selenographic longitude, positive/// towards Mare Crisium, and latitude, positive towards Mare Serenitatis.pub const Selenographic = struct { /// In (−180°, +180°]. longitude: Angle, /// In [−90°, +90°]. latitude: Angle,};
/// The Moon's libration at an instant: the selenographic coordinates of the/// center of its disk as seen from the Earth's center, which is to say the/// selenographic position of the Earth, or as seen by an observer, which is/// the observer's.pub const Libration = struct { /// What the geometry of the orbit and the tilt of the Moon's axis alone /// would give, were it rotating uniformly (eq. 53.1). optical: Selenographic, /// The real irregularities in the Moon's rotation (eq. 53.2): a few /// hundredths of a degree. Meeus says never more than 0.04°, but his /// series as printed reaches 0.061° in latitude within 1900–2100, as /// the sizes of its coefficients allow. physical: Selenographic, /// The two together, which is what an observer sees. total: Selenographic, /// The direction of the Moon's north pole on the sky, measured from the /// north point of the disk towards the east, in [0°, 360°); always within /// about 25.4° of north, so in [0°, 90°] or [270°, 360°). position_angle_of_axis: Angle,};
/// The Moon's libration at an instant, in either time scale.////// It inherits the accuracy of the Moon's position, about 10″, and of the/// physical libration's series (Eckhardt), a few thousandths of a degree./// NaN if the instant is not finite.pub fn at(instant: Instant) Libration { const c = Common.at(instant); return c.libration(c.position.ecliptic, c.position.equatorial.right_ascension);}
/// The Moon's topocentric libration at an instant, in either time scale:/// the selenographic coordinates of the center of its disk as an observer/// sees it, and the position angle of its axis there. It differs from the/// geocentric libration by up to about a degree, enough to bring/// different features onto the limb for observers far apart.////// This is Meeus's direct calculation (p. 375): the libration is found as/// `at` finds it, from the Moon's topocentric ecliptic longitude and/// latitude and its topocentric right ascension in place of the geocentric/// ones, from `astro.observer.topocentric`. NaN if the instant is not/// finite or the observer is not valid.pub fn topocentric(instant: Instant, observer: Observer) Libration { const c = Common.at(instant); const t = topocentricPosition(instant, observer, c.position.equatorial, c.position.distance_km); return c.libration(t.ecliptic, t.equatorial.right_ascension);}
/// The subsolar point at an instant, in either time scale: where on the/// Moon the Sun is overhead, in selenographic coordinates, with physical/// libration included (Meeus, p. 376).////// It is found as the libration is, from the Moon's heliocentric position/// rather than its geocentric one. It inherits the accuracy of the Sun's/// position, about 0.01°. NaN if the instant is not finite.pub fn subsolarPoint(instant: Instant) Selenographic { const c = Common.at(instant); const s = sun.position(instant); const lambda = c.position.ecliptic.longitude.radians(); const beta = c.position.ecliptic.latitude.radians(); const lambda_sun = s.ecliptic.longitude.radians();
// The Moon's heliocentric coordinates, to sufficient accuracy (p. 376): // Δ/R is the ratio of the Moon's distance to the Sun's. const ratio = c.position.distance_km / (s.distance_au * km_per_au); const lambda_h = lambda_sun + math.pi + ratio * @cos(beta) * @sin(lambda_sun - lambda); const beta_h = ratio * beta; return c.locate(lambda_h, beta_h).total();}
/// The colongitude at an instant, in either time scale: the selenographic/// longitude of the morning terminator, where the Sun is rising, counted/// as 90° minus the subsolar point's longitude, in [0°, 360°). It is near/// 0° at first quarter, 90° at full Moon, 180° at last quarter and 270° at/// new Moon, and increases by about 12.19° a day. The Sun rises on a crater/// at selenographic longitude η when the colongitude is 360° − η (Meeus,/// p. 377). NaN if the instant is not finite.pub fn colongitude(instant: Instant) Angle { return Angle.fromDegrees(90).sub(subsolarPoint(instant).longitude).normalized();}
/// What the libration and the subsolar point share: the Moon's position,/// the nutation, the node, and Eckhardt's ρ, σ and τ, which do not depend/// on whose selenographic position is wanted.const Common = struct { position: moon_position.Position, f: f64, dpsi: f64, omega: f64, obliquity: f64, rho: f64, sigma: f64, tau: f64,
fn at(instant: Instant) Common { const t = instant.julianCenturies(); const a = moon_position.Arguments.at(t); const nutation = earth.nutation(instant); const omega = deg(polynomial(t, .{ 125.0445479, -1934.1362891, 0.0020754, 1.0 / 467441.0, -1.0 / 60616000.0 }));
// Eckhardt's series for the physical libration, in degrees (p. 373). const d = deg(a.d); const m = deg(a.m); const mp = deg(a.m_moon); const f = deg(a.f); const k1 = deg(119.75 + 131.849 * t); const k2 = deg(72.56 + 20.186 * t); return .{ .position = moon_position.position(instant), .f = f, .dpsi = nutation.longitude.radians(), .omega = omega, .obliquity = earth.meanObliquity(instant).add(nutation.obliquity).radians(), .rho = -0.02752 * @cos(mp) - 0.02245 * @sin(f) + 0.00684 * @cos(mp - 2 * f) - 0.00293 * @cos(2 * f) - 0.00085 * @cos(2 * f - 2 * d) - 0.00054 * @cos(mp - 2 * d) - 0.00020 * @sin(mp + f) - 0.00020 * @cos(mp + 2 * f) - 0.00020 * @cos(mp - f) + 0.00014 * @cos(mp + 2 * f - 2 * d), .sigma = -0.02816 * @sin(mp) + 0.02244 * @cos(f) - 0.00682 * @sin(mp - 2 * f) - 0.00279 * @sin(2 * f) - 0.00083 * @sin(2 * f - 2 * d) + 0.00069 * @sin(mp - 2 * d) + 0.00040 * @cos(mp + f) - 0.00025 * @sin(2 * mp) - 0.00023 * @sin(mp + 2 * f) + 0.00020 * @cos(mp - f) + 0.00019 * @sin(mp - f) + 0.00013 * @sin(mp + 2 * f - 2 * d) - 0.00010 * @cos(mp - 3 * f), .tau = 0.02520 * a.e * @sin(m) + 0.00473 * @sin(2 * mp - 2 * f) - 0.00467 * @sin(mp) + 0.00396 * @sin(k1) + 0.00276 * @sin(2 * mp - 2 * d) + 0.00196 * @sin(omega) - 0.00183 * @cos(mp - f) + 0.00115 * @sin(mp - 2 * d) - 0.00096 * @sin(mp - d) + 0.00046 * @sin(2 * f - 2 * d) - 0.00039 * @sin(mp - f) - 0.00032 * @sin(mp - m - d) + 0.00027 * @sin(2 * mp - m - 2 * d) + 0.00023 * @sin(k2) - 0.00014 * @sin(2 * d) + 0.00014 * @cos(2 * mp - 2 * f) - 0.00012 * @sin(mp - 2 * f) - 0.00012 * @sin(2 * mp) + 0.00011 * @sin(2 * mp - 2 * m - 2 * d), }; }
/// The libration and position angle of the axis seen from wherever the /// Moon has these apparent ecliptic coordinates and this right /// ascension: the Earth's center, or an observer (pp. 372–374). fn libration(c: Common, ecliptic: coord.Ecliptic, right_ascension: Angle) Libration { const point = c.locate(ecliptic.longitude.radians(), ecliptic.latitude.radians());
// The position angle of the axis (p. 374). const v = c.omega + c.dpsi + deg(c.sigma) / @sin(inclination); const i_rho = inclination + deg(c.rho); const x = @sin(i_rho) * @sin(v); const y = @sin(i_rho) * @cos(v) * @cos(c.obliquity) - @cos(i_rho) * @sin(c.obliquity); const small_omega = math.atan2(x, y); const alpha = right_ascension.radians(); const b_total = point.b_optical + point.b_physical; const p = math.asin(clamp(@sqrt(x * x + y * y) * @cos(alpha - small_omega) / @cos(b_total)));
return .{ .optical = .{ .longitude = Angle.fromRadians(point.l_optical).signed(), .latitude = Angle.fromRadians(point.b_optical), }, .physical = .{ .longitude = Angle.fromRadians(point.l_physical).signed(), .latitude = Angle.fromRadians(point.b_physical), }, .total = point.total(), .position_angle_of_axis = Angle.fromRadians(p).normalized(), }; }
/// The selenographic position of whatever is seen in the direction of /// ecliptic longitude λ and latitude β from the Moon's center — the /// Earth from the Moon's geocentric position, the Sun from its /// heliocentric one — in radians (eqs. 53.1 and 53.2). fn locate(self: Common, lambda: f64, beta: f64) Point { const w = lambda - self.dpsi - self.omega; const big_a = math.atan2(@sin(w) * @cos(beta) * @cos(inclination) - @sin(beta) * @sin(inclination), @cos(w) * @cos(beta)); const b_optical = math.asin(clamp(-@sin(w) * @cos(beta) * @sin(inclination) - @sin(beta) * @cos(inclination))); return .{ .l_optical = big_a - self.f, .b_optical = b_optical, .l_physical = deg(-self.tau + (self.rho * @cos(big_a) + self.sigma * @sin(big_a)) * @tan(b_optical)), .b_physical = deg(self.sigma * @cos(big_a) - self.rho * @sin(big_a)), }; }};
const Point = struct { l_optical: f64, b_optical: f64, l_physical: f64, b_physical: f64,
fn total(self: Point) Selenographic { return .{ .longitude = Angle.fromRadians(self.l_optical + self.l_physical).signed(), .latitude = Angle.fromRadians(self.b_optical + self.b_physical), }; }};
const km_per_au = 149597870.7;
/// The inclination of the mean lunar equator to the ecliptic, 1°32′32.7″,/// the IAU's value (p. 372).const inclination = 1.54242 * math.rad_per_deg;
fn deg(degrees: f64) f64 { return degrees * math.rad_per_deg;}
/// Rounding can carry a sine a hair past ±1; NaN stays NaN.fn clamp(x: f64) f64 { if (x > 1) return 1; if (x < -1) return -1; return x;}
/// 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;}