// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>
// SPDX-License-Identifier: MIT

//! Positions on the celestial sphere, and the conversions between the
//! ecliptic, equatorial and horizontal frames (Meeus, *Astronomical
//! Algorithms*, 2nd ed., ch. 13).
//!
//! Neither type says which equinox it is referred to; that is a property of
//! where the coordinates came from, and every function returning one says.
//! A conversion takes the obliquity belonging to that equinox: the true
//! obliquity for apparent coordinates of date, from `astro.earth`.

const std = @import("std");
const math = std.math;

const Angle = @import("Angle.zig");

/// A direction measured from the ecliptic and the vernal equinox.
pub const Ecliptic = struct {
    /// Celestial longitude λ, measured eastwards along the ecliptic from the
    /// equinox, in [0°, 360°).
    longitude: Angle,
    /// Celestial latitude β, positive north of the ecliptic, in
    /// [−90°, +90°].
    latitude: Angle,

    /// The same direction in equatorial coordinates (Meeus, eqs. 13.3 and
    /// 13.4).
    pub fn toEquatorial(self: Ecliptic, obliquity: Angle) Equatorial {
        const l = self.longitude.radians();
        const b = self.latitude.radians();
        const e = obliquity.radians();
        // atan2 rather than the book's tangent, which loses the quadrant.
        const ra = math.atan2(@sin(l) * @cos(e) - @tan(b) * @sin(e), @cos(l));
        const dec = math.asin(clamp(@sin(b) * @cos(e) + @cos(b) * @sin(e) * @sin(l)));
        return .{
            .right_ascension = Angle.fromRadians(ra).normalized(),
            .declination = Angle.fromRadians(dec),
        };
    }
};

/// A direction measured from the celestial equator and the vernal equinox.
pub const Equatorial = struct {
    /// Right ascension α, measured eastwards along the equator from the
    /// equinox, in [0°, 360°). `Angle.hours` reads it in hours.
    right_ascension: Angle,
    /// Declination δ, positive north of the equator, in [−90°, +90°].
    declination: Angle,

    /// The same direction in ecliptic coordinates (Meeus, eqs. 13.1 and
    /// 13.2).
    pub fn toEcliptic(self: Equatorial, obliquity: Angle) Ecliptic {
        const a = self.right_ascension.radians();
        const d = self.declination.radians();
        const e = obliquity.radians();
        const lon = math.atan2(@sin(a) * @cos(e) + @tan(d) * @sin(e), @cos(a));
        const lat = math.asin(clamp(@sin(d) * @cos(e) - @cos(d) * @sin(e) * @sin(a)));
        return .{
            .longitude = Angle.fromRadians(lon).normalized(),
            .latitude = Angle.fromRadians(lat),
        };
    }

    /// The same direction in an observer's horizontal coordinates, given
    /// its hour angle there (positive west of the meridian) and the
    /// observer's geodetic latitude (Meeus, eqs. 13.5 and 13.6, with the
    /// azimuth turned to count from the north; see `Horizontal`). Give it
    /// topocentric coordinates and the topocentric hour angle, or the
    /// result is out by the body's parallax. No refraction is applied: the
    /// altitude is the airless one.
    pub fn toHorizontal(self: Equatorial, hour_angle: Angle, latitude: Angle) Horizontal {
        const h = hour_angle.radians();
        const d = self.declination.radians();
        const phi = latitude.radians();
        // The direction's components towards the north point, the east
        // point and the zenith. Meeus's tan A = sin H / (cos H sin φ −
        // tan δ cos φ), from the south, is the first two multiplied
        // through by cos δ, so that a body at a pole needs no tangent, and
        // with both signs reversed to count from the north. His
        // sin h (eq. 13.6) is the third; the altitude is taken by atan2
        // rather than asin, which loses a millionth of a degree near the
        // zenith.
        const north = @sin(d) * @cos(phi) - @cos(d) * @cos(h) * @sin(phi);
        const east = -@cos(d) * @sin(h);
        const up = @sin(phi) * @sin(d) + @cos(phi) * @cos(d) * @cos(h);
        const azimuth = math.atan2(east, north);
        const altitude = math.atan2(up, @sqrt(north * north + east * east));
        return .{
            .azimuth = Angle.fromRadians(azimuth).normalized(),
            .altitude = Angle.fromRadians(altitude),
        };
    }
};

/// A direction as an observer on the Earth sees it, measured from the
/// horizon and the north point.
///
/// Azimuth counts from the north through the east: 0° north, 90° east,
/// 180° south, 270° west, as compasses, almanacs and JPL Horizons count it.
/// Meeus counts from the *south*, westwards (p. 91); add 180° to his.
pub const Horizontal = struct {
    /// Azimuth, from the north through the east, in [0°, 360°). At the
    /// zenith and the nadir it has no meaning, and comes out as whatever
    /// the rounding gives.
    azimuth: Angle,
    /// Altitude above the horizon, negative below it, in [−90°, +90°]. It
    /// is the airless altitude unless it says otherwise: the atmosphere
    /// raises what is seen by up to about 0.6°, at the horizon, and
    /// `astro.observer.observedAltitude` says by how much.
    altitude: Angle,
};

/// Rounding can carry a sine a hair past ±1 at the poles, where `asin`
/// would give NaN for a perfectly good direction. NaN stays NaN.
fn clamp(x: f64) f64 {
    if (x > 1) return 1;
    if (x < -1) return -1;
    return x;
}