// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie // SPDX-License-Identifier: MIT //! The Moon's phase: how much of it is lit and from which side, at an //! instant (Meeus, *Astronomical Algorithms*, 2nd ed., ch. 48). const std = @import("std"); const math = std.math; const Angle = @import("../Angle.zig"); const Instant = @import("../time.zig").Instant; const sun = @import("../sun.zig"); const moon_position = @import("position.zig"); /// The conventional name for the Moon's appearance, by which 45° sector of /// elongation in longitude it is in. Each of the four principal phases is /// named for the sector centered on it, and the four between are centered /// on 45°, 135°, 225° and 315°. pub const NamedPhase = enum { new_moon, waxing_crescent, first_quarter, waxing_gibbous, full_moon, waning_gibbous, last_quarter, waning_crescent, }; /// The Moon's phase at an instant, as seen from the Earth's center. pub const Phase = struct { /// The angle between the Moon and the Sun, from 0° to 180°. elongation: Angle, /// The Moon's apparent ecliptic longitude minus the Sun's, in /// [0°, 360°): 0° at new Moon, 90° at first quarter, 180° at full Moon /// and 270° at last quarter. elongation_in_longitude: Angle, /// The angle between the Sun and the Earth as seen from the Moon, from /// 0° (full) to 180° (new). phase_angle: Angle, /// The fraction of the Moon's disk that is lit, from 0 to 1. illuminated_fraction: f64, /// The position angle of the midpoint of the bright limb, in /// [0°, 360°), measured from the north point of the disk towards the /// east — from the direction of the celestial pole, not of the Moon's /// axis. Near 270° at first quarter and near 90° after full Moon. bright_limb: Angle, /// True while the Moon is waxing (elongation in longitude under 180°), /// false while it is waning; null exactly at new or full Moon, or for a /// non-finite instant. pub fn waxing(self: Phase) ?bool { const l = self.elongation_in_longitude.degrees(); if (l > 0 and l < 180) return true; if (l > 180 and l < 360) return false; return null; } /// The named phase. Each sector includes its lower boundary and not its /// upper: an elongation in longitude of exactly 22.5° is a waxing /// crescent. Null for a non-finite instant. pub fn named(self: Phase) ?NamedPhase { const l = self.elongation_in_longitude.degrees(); if (!(l >= 0 and l < 360)) return null; const sector: u3 = @intFromFloat(@mod(@floor((l + 22.5) / 45), 8)); return @fromBackingInt(sector); } }; /// The Moon's phase at an instant, in either time scale, from the apparent /// positions of the Moon and the Sun (Meeus, eqs. 48.1–48.3 and 48.5). /// /// The Sun's position is good to about 0.01°, and that is what limits the /// elongation and the phase angle; Meeus notes that 1′ is ample for the /// illuminated fraction, which is good to about 0.0001. /// /// NaN in every field if the instant is not finite. pub fn at(instant: Instant) Phase { const m = moon_position.position(instant); const s = sun.position(instant); const a = m.equatorial.right_ascension.radians(); const d = m.equatorial.declination.radians(); const a0 = s.equatorial.right_ascension.radians(); const d0 = s.equatorial.declination.radians(); // The elongation (eq. 48.2), by atan2 of the sine and cosine rather than // the arccosine, which loses precision near new and full Moon. const cos_psi = @sin(d0) * @sin(d) + @cos(d0) * @cos(d) * @cos(a0 - a); const sin_psi = sinOfSeparation(a, d, a0, d0); const psi = math.atan2(sin_psi, cos_psi); // The phase angle (eq. 48.3), with both distances in kilometers. const r = s.distance_au * km_per_au; const i = math.atan2(r * sin_psi, m.distance_km - r * cos_psi); // The bright limb (eq. 48.5). const chi = math.atan2( @cos(d0) * @sin(a0 - a), @sin(d0) * @cos(d) - @cos(d0) * @sin(d) * @cos(a0 - a), ); return .{ .elongation = Angle.fromRadians(psi), .elongation_in_longitude = m.ecliptic.longitude.sub(s.ecliptic.longitude).normalized(), .phase_angle = Angle.fromRadians(i), .illuminated_fraction = (1 + @cos(i)) / 2, .bright_limb = Angle.fromRadians(chi).normalized(), }; } const km_per_au = 149597870.7; /// The sine of the angle between two directions, from the length of the /// cross product of their unit vectors: always positive, and exact where /// the angle is small. fn sinOfSeparation(a: f64, d: f64, a0: f64, d0: f64) f64 { const x = @cos(d) * @cos(a); const y = @cos(d) * @sin(a); const z = @sin(d); const x0 = @cos(d0) * @cos(a0); const y0 = @cos(d0) * @sin(a0); const z0 = @sin(d0); const cx = y * z0 - z * y0; const cy = z * x0 - x * z0; const cz = x * y0 - y * x0; return @sqrt(cx * cx + cy * cy + cz * cz); } /// One of the four instants in each lunation when the Moon's apparent /// geocentric longitude exceeds the Sun's by exactly 0°, 90°, 180° or 270°. pub const PrincipalPhase = enum(u2) { new_moon, first_quarter, full_moon, last_quarter, /// The fraction of a lunation after new Moon at which this phase falls /// in Meeus's numbering. fn offset(self: PrincipalPhase) f64 { return @as(f64, @floatFromInt(@backingInt(self))) / 4; } }; /// A principal phase and when it happens. pub const Event = struct { phase: PrincipalPhase, /// The instant, in TT. instant: Instant, }; /// The first principal phase of `kind` after `instant`, in TT (Meeus, ch. /// 49). /// /// Meeus compared the method with the complete ELP-2000/82 and VSOP87 /// theories over 1980–2020 and found a mean error of 3.7 s and a largest of /// 17.4 s. Like the theory it comes from, it grows less accurate with /// distance from the present. /// /// NaN if the instant is not finite, or more than about 100,000 years from /// 2000. Beyond that the series' factor for the Earth's eccentricity, which /// is a polynomial in time, grows until its corrections exceed the week /// between phases, and it stops producing phases in order at all. pub fn next(instant: Instant, kind: PrincipalPhase) Instant { return .{ .julian_day = search(instant, kind, .after), .scale = .tt }; } /// The last principal phase of `kind` before `instant`, in TT. NaN where /// `next` would be. pub fn previous(instant: Instant, kind: PrincipalPhase) Instant { return .{ .julian_day = search(instant, kind, .before), .scale = .tt }; } /// The time since the most recent new Moon at or before `instant`, in days: /// zero at the instant of new Moon itself, and never more than the longest /// lunation, about 29.83 days. NaN where `next` would be. pub fn lunarAgeDays(instant: Instant) f64 { const jd = instant.to(.tt, .model).julian_day; return jd - search(instant, .new_moon, .at_or_before); } /// Every principal phase from `from` (inclusive) to `until` (exclusive), in /// order. An infinite `until` iterates to the last phase `next` will compute, /// some 100,000 years on; a non-finite `from` yields nothing. pub fn between(from: Instant, until: Instant) Iterator { const start = from.to(.tt, .model).julian_day; // `Instant.to` makes every non-finite Julian Day NaN, so +∞ is kept // aside; NaN and −∞ yield nothing. const end = if (math.isPositiveInf(until.julian_day)) math.inf(f64) else until.to(.tt, .model).julian_day; // The quarter-lunation before `from`, so that the first phase at or // after it is not missed; `Iterator.next` skips any that are earlier. const quarters = @floor((start - mean_epoch) / mean_lunation * 4) - 2; return .{ .quarter = if (math.isFinite(quarters) and @abs(quarters) < max_quarters) quarters else math.nan(f64), .start = start, .end = end, }; } pub const Iterator = struct { /// The next lunation number to try, in quarters: k × 4. quarter: f64, start: f64, end: f64, pub fn next(self: *Iterator) ?Event { while (math.isFinite(self.quarter)) { // Past the last lunation computed, there is nothing more to // find, and nothing to skip forward through. if (@abs(self.quarter) > max_quarters) { self.quarter = math.nan(f64); return null; } const k = self.quarter / 4; self.quarter += 1; const jde = trueJde(k); if (!(jde >= self.start)) continue; if (!(jde < self.end)) { self.quarter = math.nan(f64); return null; } const kind: PrincipalPhase = @fromBackingInt(@as(u2, @intFromFloat(@mod(self.quarter - 1, 4)))); return .{ .phase = kind, .instant = .{ .julian_day = jde, .scale = .tt } }; } return null; } }; /// The mean new Moon of 2000 January 6, k = 0 (eq. 49.1). const mean_epoch = 2451550.09766; const mean_lunation = 29.530588861; /// The furthest lunation from 2000 for which a phase is computed: 1,000 /// centuries, where T = k / 1236.85 is 1000. There E² is about 80, the /// periodic corrections reach a few days, and the phases still come in /// order; some way beyond, they do not. const max_lunations = 1236.85 * 1000; const max_quarters = max_lunations * 4; const Direction = enum { after, before, at_or_before }; /// The principal phase of `kind` nearest `instant` in `direction`. The true /// phase is never more than about 0.6 days from the mean one and phases of a /// kind are 29.5 days apart, so it is one of the four lunations around the /// mean estimate: no loop, and nothing to run away on. fn search(instant: Instant, kind: PrincipalPhase, direction: Direction) f64 { const jd = instant.to(.tt, .model).julian_day; const estimate = @floor((jd - mean_epoch) / mean_lunation - kind.offset()); if (!math.isFinite(estimate) or @abs(estimate) > max_lunations) return math.nan(f64); var best = math.nan(f64); for ([_]f64{ -1, 0, 1, 2 }) |step| { const jde = trueJde(estimate + step + kind.offset()); const ok = switch (direction) { .after => jde > jd, .before => jde < jd, .at_or_before => jde <= jd, }; if (!ok) continue; if (math.isNan(best)) { best = jde; } else if (direction == .after) { best = @min(best, jde); } else { best = @max(best, jde); } } return best; } /// The instant of the true principal phase with lunation number `k`, as a /// Julian Ephemeris Day (Meeus, ch. 49). `k` is an integer for new Moon and /// an integer plus 0.25, 0.5 or 0.75 for the others. fn trueJde(k: f64) f64 { const t = k / 1236.85; const t2 = t * t; const t3 = t2 * t; const t4 = t3 * t; const mean = mean_epoch + mean_lunation * k + 0.00015437 * t2 - 0.000000150 * t3 + 0.00000000073 * t4; const e = 1 - 0.002516 * t - 0.0000074 * t2; const m = deg(2.5534 + 29.10535670 * k - 0.0000014 * t2 - 0.00000011 * t3); const mp = deg(201.5643 + 385.81693528 * k + 0.0107582 * t2 + 0.00001238 * t3 - 0.000000058 * t4); const f = deg(160.7108 + 390.67050284 * k - 0.0016118 * t2 - 0.00000227 * t3 + 0.000000011 * t4); const omega = deg(124.7746 - 1.56375588 * k + 0.0020672 * t2 + 0.00000215 * t3); const fraction = k - @floor(k); var correction: f64 = 0; if (fraction < 0.125 or fraction > 0.875) { correction = periodic(&new_moon_terms, e, m, mp, f, omega); } else if (fraction > 0.375 and fraction < 0.625) { correction = periodic(&full_moon_terms, e, m, mp, f, omega); } else { correction = periodic(&quarter_terms, e, m, mp, f, omega); const w = 0.00306 - 0.00038 * e * @cos(m) + 0.00026 * @cos(mp) - 0.00002 * @cos(mp - m) + 0.00002 * @cos(mp + m) + 0.00002 * @cos(2 * f); correction += if (fraction < 0.5) w else -w; } // The planetary arguments and their corrections, common to all phases. const a = [14]f64{ deg(299.77 + 0.107408 * k - 0.009173 * t2), deg(251.88 + 0.016321 * k), deg(251.83 + 26.651886 * k), deg(349.42 + 36.412478 * k), deg(84.66 + 18.206239 * k), deg(141.74 + 53.303771 * k), deg(207.14 + 2.453732 * k), deg(154.84 + 7.306860 * k), deg(34.52 + 27.261239 * k), deg(207.19 + 0.121824 * k), deg(291.34 + 1.844379 * k), deg(161.72 + 24.198154 * k), deg(239.56 + 25.513099 * k), deg(331.55 + 3.592518 * k), }; const planetary = [14]f64{ 325, 165, 164, 126, 110, 62, 60, 56, 47, 42, 40, 37, 35, 23 }; for (a, planetary) |angle, coefficient| correction += coefficient * 1e-6 * @sin(angle); return mean + correction; } fn deg(degrees: f64) f64 { return degrees * math.rad_per_deg; } /// A periodic correction to the mean phase, in days: the coefficient, the /// power of E it is multiplied by, and the multiples of M, M′, F and Ω in /// its argument. const Term = struct { c: f64, e: u2 = 0, m: i8 = 0, mp: i8 = 0, f: i8 = 0, omega: i8 = 0 }; fn periodic(comptime terms: []const Term, e: f64, m: f64, mp: f64, f: f64, omega: f64) f64 { var sum: f64 = 0; inline for (terms) |term| { const argument = @as(f64, @floatFromInt(term.m)) * m + @as(f64, @floatFromInt(term.mp)) * mp + @as(f64, @floatFromInt(term.f)) * f + @as(f64, @floatFromInt(term.omega)) * omega; const factor = switch (term.e) { 0 => 1, 1 => e, else => e * e, }; sum += term.c * factor * @sin(argument); } return sum; } /// Meeus, p. 351, the New Moon column. const new_moon_terms = [_]Term{ .{ .c = -0.40720, .mp = 1 }, .{ .c = 0.17241, .e = 1, .m = 1 }, .{ .c = 0.01608, .mp = 2 }, .{ .c = 0.01039, .f = 2 }, .{ .c = 0.00739, .e = 1, .mp = 1, .m = -1 }, .{ .c = -0.00514, .e = 1, .mp = 1, .m = 1 }, .{ .c = 0.00208, .e = 2, .m = 2 }, .{ .c = -0.00111, .mp = 1, .f = -2 }, .{ .c = -0.00057, .mp = 1, .f = 2 }, .{ .c = 0.00056, .e = 1, .mp = 2, .m = 1 }, .{ .c = -0.00042, .mp = 3 }, .{ .c = 0.00042, .e = 1, .m = 1, .f = 2 }, .{ .c = 0.00038, .e = 1, .m = 1, .f = -2 }, .{ .c = -0.00024, .e = 1, .mp = 2, .m = -1 }, .{ .c = -0.00017, .omega = 1 }, .{ .c = -0.00007, .mp = 1, .m = 2 }, .{ .c = 0.00004, .mp = 2, .f = -2 }, .{ .c = 0.00004, .m = 3 }, .{ .c = 0.00003, .mp = 1, .m = 1, .f = -2 }, .{ .c = 0.00003, .mp = 2, .f = 2 }, .{ .c = -0.00003, .mp = 1, .m = 1, .f = 2 }, .{ .c = 0.00003, .mp = 1, .m = -1, .f = 2 }, .{ .c = -0.00002, .mp = 1, .m = -1, .f = -2 }, .{ .c = -0.00002, .mp = 3, .m = 1 }, .{ .c = 0.00002, .mp = 4 }, }; /// Meeus, p. 351, the Full Moon column: the same arguments as New Moon. const full_moon_terms = [_]Term{ .{ .c = -0.40614, .mp = 1 }, .{ .c = 0.17302, .e = 1, .m = 1 }, .{ .c = 0.01614, .mp = 2 }, .{ .c = 0.01043, .f = 2 }, .{ .c = 0.00734, .e = 1, .mp = 1, .m = -1 }, .{ .c = -0.00515, .e = 1, .mp = 1, .m = 1 }, .{ .c = 0.00209, .e = 2, .m = 2 }, .{ .c = -0.00111, .mp = 1, .f = -2 }, .{ .c = -0.00057, .mp = 1, .f = 2 }, .{ .c = 0.00056, .e = 1, .mp = 2, .m = 1 }, .{ .c = -0.00042, .mp = 3 }, .{ .c = 0.00042, .e = 1, .m = 1, .f = 2 }, .{ .c = 0.00038, .e = 1, .m = 1, .f = -2 }, .{ .c = -0.00024, .e = 1, .mp = 2, .m = -1 }, .{ .c = -0.00017, .omega = 1 }, .{ .c = -0.00007, .mp = 1, .m = 2 }, .{ .c = 0.00004, .mp = 2, .f = -2 }, .{ .c = 0.00004, .m = 3 }, .{ .c = 0.00003, .mp = 1, .m = 1, .f = -2 }, .{ .c = 0.00003, .mp = 2, .f = 2 }, .{ .c = -0.00003, .mp = 1, .m = 1, .f = 2 }, .{ .c = 0.00003, .mp = 1, .m = -1, .f = 2 }, .{ .c = -0.00002, .mp = 1, .m = -1, .f = -2 }, .{ .c = -0.00002, .mp = 3, .m = 1 }, .{ .c = 0.00002, .mp = 4 }, }; /// Meeus, p. 352, First and Last Quarters. const quarter_terms = [_]Term{ .{ .c = -0.62801, .mp = 1 }, .{ .c = 0.17172, .e = 1, .m = 1 }, .{ .c = -0.01183, .e = 1, .mp = 1, .m = 1 }, .{ .c = 0.00862, .mp = 2 }, .{ .c = 0.00804, .f = 2 }, .{ .c = 0.00454, .e = 1, .mp = 1, .m = -1 }, .{ .c = 0.00204, .e = 2, .m = 2 }, .{ .c = -0.00180, .mp = 1, .f = -2 }, .{ .c = -0.00070, .mp = 1, .f = 2 }, .{ .c = -0.00040, .mp = 3 }, .{ .c = -0.00034, .e = 1, .mp = 2, .m = -1 }, .{ .c = 0.00032, .e = 1, .m = 1, .f = 2 }, .{ .c = 0.00032, .e = 1, .m = 1, .f = -2 }, .{ .c = -0.00028, .e = 2, .mp = 1, .m = 2 }, .{ .c = 0.00027, .e = 1, .mp = 2, .m = 1 }, .{ .c = -0.00017, .omega = 1 }, .{ .c = -0.00005, .mp = 1, .m = -1, .f = -2 }, .{ .c = 0.00004, .mp = 2, .f = 2 }, .{ .c = -0.00004, .mp = 1, .m = 1, .f = 2 }, .{ .c = 0.00004, .mp = 1, .m = -2 }, .{ .c = 0.00003, .mp = 1, .m = 1, .f = -2 }, .{ .c = 0.00003, .m = 3 }, .{ .c = 0.00002, .mp = 2, .f = -2 }, .{ .c = 0.00002, .mp = 1, .m = -1, .f = 2 }, .{ .c = -0.00002, .mp = 3, .m = 1 }, };