// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>
// 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 },
};