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

//! An observer on the Earth's surface, and what a body looks like from
//! there: its topocentric position, shifted from the geocentric one by its
//! parallax (Meeus, *Astronomical Algorithms*, 2nd ed., chs. 11 and 40),
//! and the refraction that lifts it from its airless altitude to where it
//! is seen (ch. 16).
//!
//! Longitude counts positive east of Greenwich, as the IAU, WGS 84 and JPL
//! Horizons have it. Meeus counts it positive *west*; his longitudes change
//! sign here (see ADR 0003).
//!
//! Left out, each well inside the accuracy of the Sun's and the Moon's
//! positions: diurnal aberration, at most 0.32″; polar motion, a few tenths
//! of a second of arc; the difference between UTC and UT1, under 0.9 s of
//! the Earth's rotation; and the deflection of the vertical, by which the
//! local plumb line departs from the ellipsoid's normal by up to about a
//! minute of arc in mountains.

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;

/// A point on the Earth's surface.
pub const Observer = struct {
    /// Geodetic latitude, the latitude on maps: the angle between the
    /// equator and the normal to the WGS 84 ellipsoid, positive north. It
    /// must be within [−90°, +90°]; beyond that, every result is NaN.
    latitude: Angle,
    /// Longitude, positive east of Greenwich. Any finite value is taken
    /// modulo 360°.
    longitude: Angle,
    /// Height above the WGS 84 ellipsoid, in meters. Height above sea level
    /// differs from it by the geoid, at most about 100 m, which moves the
    /// Moon by about 0.003″. Any finite height is accepted, including
    /// negative ones.
    height_m: f64 = 0,

    /// The observer's place in the Earth's meridian plane: ρ sin φ′ and
    /// ρ cos φ′, where ρ is the distance from the Earth's center in units of
    /// the equatorial radius and φ′ the geocentric latitude (Meeus, p. 82,
    /// with the WGS 84 ellipsoid in place of the IAU 1976 one). NaN if the
    /// latitude is outside [−90°, +90°] or anything is not finite.
    pub fn geocentric(self: Observer) Geocentric {
        if (!self.isValid()) return .{ .rho_sin_phi = math.nan(f64), .rho_cos_phi = math.nan(f64) };
        const phi = self.latitude.radians();
        const u = math.atan2(polar_ratio * @sin(phi), @cos(phi));
        const h = self.height_m / (equatorial_radius_km * 1000);
        return .{
            .rho_sin_phi = polar_ratio * @sin(u) + h * @sin(phi),
            .rho_cos_phi = @cos(u) + h * @cos(phi),
        };
    }

    fn isValid(self: Observer) bool {
        return @abs(self.latitude.radians()) <= math.pi / 2.0 and
            math.isFinite(self.longitude.radians()) and
            math.isFinite(self.height_m);
    }
};

/// An observer's place in the Earth's meridian plane; see
/// `Observer.geocentric`.
pub const Geocentric = struct {
    /// ρ sin φ′: positive in the northern hemisphere.
    rho_sin_phi: f64,
    /// ρ cos φ′: never negative.
    rho_cos_phi: f64,
};

/// Where a body appears to an observer.
pub const Topocentric = struct {
    /// Topocentric right ascension and declination, referred to the true
    /// equator and equinox of date.
    equatorial: coord.Equatorial,
    /// Topocentric ecliptic longitude and latitude, referred to the true
    /// equinox of date.
    ecliptic: coord.Ecliptic,
    /// The topocentric hour angle: how far west of the observer's meridian
    /// the body is, in (−180°, +180°]; negative while it is still rising
    /// towards the meridian.
    hour_angle: Angle,
    /// Azimuth, from the north through the east, and the airless altitude:
    /// no refraction is applied. `observedAltitude` gives where the body is
    /// seen.
    horizontal: coord.Horizontal,
    /// The distance from the observer to the body's center, in kilometers.
    distance_km: f64,
};

/// The equatorial radius of the WGS 84 ellipsoid, in kilometers.
pub const equatorial_radius_km = 6378.137;

/// b/a = 1 − f for the WGS 84 ellipsoid, f = 1/298.257223563.
const polar_ratio = 1 - 1.0 / 298.257223563;

/// A body's topocentric position, from its apparent geocentric equatorial
/// coordinates and distance at an instant, which may be in either time
/// scale. This is what every body's own `topocentric` calls, and what a
/// body without one can be passed through.
///
/// The parallax is applied rigorously, as the difference of the body's and
/// the observer's geocentric vectors. That is Meeus's eqs. 40.6 and 40.7,
/// with sin π the ratio of the Earth's equatorial radius to the distance,
/// and it gives the topocentric distance as well. The hour angle comes from
/// the apparent sidereal time, the ecliptic coordinates from the
/// topocentric equatorial ones with the true obliquity, and the horizontal
/// ones from the topocentric hour angle and declination with the
/// observer's geodetic latitude (Meeus, eqs. 13.5 and 13.6).
///
/// NaN if the instant, the coordinates or the distance is not finite, or the
/// observer is not valid (see `Observer`).
pub fn topocentric(
    instant: Instant,
    observer: Observer,
    equatorial: coord.Equatorial,
    distance_km: f64,
) Topocentric {
    if (!observer.isValid()) return nan_topocentric;
    const place = observer.geocentric();
    // Reduced before anything is subtracted from it: a longitude of 10³⁰⁰°
    // is finite, and would swallow any right ascension.
    const local_sidereal_time = earth.greenwichApparentSiderealTime(instant).add(observer.longitude).normalized().radians();
    const alpha = equatorial.right_ascension.radians();
    const delta = equatorial.declination.radians();
    const hour_angle = local_sidereal_time - alpha;

    // The body as seen from the Earth's center, in kilometers, in a frame
    // whose x axis points to the observer's meridian on the equator, y to
    // the west point and z to the north celestial pole; then the observer's
    // own position taken away.
    const x = distance_km * @cos(delta) * @cos(hour_angle) - equatorial_radius_km * place.rho_cos_phi;
    const y = distance_km * @cos(delta) * @sin(hour_angle);
    const z = distance_km * @sin(delta) - equatorial_radius_km * place.rho_sin_phi;

    const topocentric_hour_angle = math.atan2(y, x);
    const topocentric_equatorial: coord.Equatorial = .{
        .right_ascension = Angle.fromRadians(local_sidereal_time - topocentric_hour_angle).normalized(),
        .declination = Angle.fromRadians(math.atan2(z, @sqrt(x * x + y * y))),
    };
    const h = Angle.fromRadians(topocentric_hour_angle).signed();
    return .{
        .equatorial = topocentric_equatorial,
        .ecliptic = topocentric_equatorial.toEcliptic(earth.trueObliquity(instant)),
        .hour_angle = h,
        .horizontal = topocentric_equatorial.toHorizontal(h, observer.latitude),
        .distance_km = @sqrt(x * x + y * y + z * z),
    };
}

/// The air an observer looks through, which sets how strongly it
/// refracts. The defaults are the conditions the refraction formulas were
/// written for: sea level, 1010 mb and 10 °C.
pub const Atmosphere = struct {
    /// Pressure at the observer, in millibars (hectopascals).
    pressure_mb: f64 = 1010,
    /// Air temperature at the observer, in degrees Celsius.
    temperature_c: f64 = 10,

    /// How much stronger or weaker refraction is than under the default
    /// conditions: P/1010 × 283/(273 + T) (Meeus, p. 107). Meeus calls it
    /// only approximately correct, and refraction near the horizon
    /// wanders by tenths of a degree with the weather whatever it is
    /// given.
    fn factor(self: Atmosphere) f64 {
        return self.pressure_mb / 1010 * 283 / (273 + self.temperature_c);
    }
};

/// Where a body at a given airless altitude is seen: the altitude plus
/// the refraction, by Sæmundsson's formula (Meeus, eq. 16.4), with
/// Meeus's constant that makes it zero at the zenith, scaled for the
/// atmosphere.
///
/// Refraction lifts a body by about 29′ at the horizon, 5.4′ at 10° and
/// 1′ at 45°. Below an airless altitude of −1° it is held at its value
/// there, 38.80′ in the default atmosphere, as JPL Horizons holds it: the
/// formula means nothing so far below the horizon, and holding it keeps
/// the observed altitude increasing with the airless one, so that
/// `airlessAltitude` has an answer for every observed altitude. Refraction
/// a few degrees from the horizon is uncertain by much more than the
/// formula's accuracy; see Meeus, p. 107.
///
/// NaN if the altitude is not finite.
pub fn observedAltitude(airless: Angle, atmosphere: Atmosphere) Angle {
    if (!math.isFinite(airless.radians())) return nan;
    const h = airless.degrees();
    return Angle.fromDegrees(h + atmosphere.factor() * refraction(h) / 60);
}

/// Where a body seen at a given observed altitude is without the
/// atmosphere: the exact inverse of `observedAltitude`, in any atmosphere,
/// found by Newton's method.
///
/// Meeus gives Bennett's formula (eq. 16.3) for this direction, but it and
/// Sæmundsson's are separate fits that agree only to about 4″ in the
/// default atmosphere, and once the refraction is scaled for any other
/// they are evaluated at different altitudes and part by minutes of arc
/// near the horizon. Inverting Sæmundsson's keeps the two directions
/// consistent, and consistent with JPL Horizons.
///
/// NaN if the altitude is not finite.
pub fn airlessAltitude(observed: Angle, atmosphere: Atmosphere) Angle {
    if (!math.isFinite(observed.radians())) return nan;
    const h0 = observed.degrees();
    const f = atmosphere.factor();
    // Refraction changes by at most about 0.17° a degree of altitude in
    // the default atmosphere, so h + f R(h) climbs steadily and Newton's
    // method, started from h₀ less the refraction there, converges in a
    // few steps; below −1°, where refraction is held, in one.
    var h = h0 - f * refraction(h0) / 60;
    for (0..8) |_| {
        const step = (h + f * refraction(h) / 60 - h0) / (1 + f * refractionSlope(h) / 60);
        h -= step;
        if (@abs(step) < 1e-12) break;
    }
    return Angle.fromDegrees(h);
}

/// Sæmundsson's refraction at an airless altitude h in degrees, in the
/// default atmosphere, in minutes of arc, held below −1°.
fn refraction(h: f64) f64 {
    const a = @max(h, hold_below_degrees);
    const r = 1.02 / @tan(saemundssonArgument(a) * math.rad_per_deg) + 0.0019279;
    return @max(r, 0);
}

/// The derivative of `refraction`, in minutes of arc a degree.
fn refractionSlope(h: f64) f64 {
    if (h < hold_below_degrees) return 0;
    const x = saemundssonArgument(h) * math.rad_per_deg;
    const s = @sin(x);
    const dx_dh = 1 - 10.3 / ((h + 5.11) * (h + 5.11));
    return -1.02 * math.rad_per_deg * dx_dh / (s * s);
}

/// h + 10.3/(h + 5.11), in degrees.
fn saemundssonArgument(h: f64) f64 {
    return h + 10.3 / (h + 5.11);
}

/// The airless altitude in degrees below which refraction is held.
const hold_below_degrees: f64 = -1.0;

const nan = Angle.fromRadians(math.nan(f64));

const nan_topocentric: Topocentric = .{
    .equatorial = .{ .right_ascension = nan, .declination = nan },
    .ecliptic = .{ .longitude = nan, .latitude = nan },
    .hour_angle = nan,
    .horizontal = .{ .azimuth = nan, .altitude = nan },
    .distance_km = math.nan(f64),
};