// 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;
}