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

//! Refraction in `astro.observer`, against Meeus's worked example and JPL
//! Horizons.

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

const astro = @import("astro");
const Angle = astro.Angle;
const observer = astro.observer;
const horizons = @import("horizons.zig");

fn arcmin(a: Angle) f64 {
    return a.degrees() * 60;
}

// Meeus, Astronomical Algorithms, 2nd ed., example 16.a (p. 107): the Sun's
// disk at the horizon, with its lower limb seen at exactly 0°30′. The book
// goes back to the airless altitude by Bennett's formula (eq. 16.3), and
// forward again by Sæmundsson's (eq. 16.4), both without the constants
// that make them zero at the zenith; it prints to 0.001′.
//
// `airlessAltitude` inverts Sæmundsson's formula instead, so its 28.812′
// is not Bennett's 28.754′: the two formulas are separate fits, which
// Meeus says agree to within 4″, 0.067′, and here are 0.058′ apart. The
// rest of the example starts from the book's own figures.
test "example 16.a: the flattened Sun at the horizon" {
    const lower = observer.airlessAltitude(Angle.fromDegrees(0.5), .{});
    try testing.expectApproxEqAbs(28.754, 30 - arcmin(lower), 0.059);
    // The upper limb is 32′ above the book's 0°01.246′, at 0°33.246′; by
    // Sæmundsson, R = 24.618′ there, so the upper limb is seen at 57.864′,
    // and the vertical diameter at 27.864′, 0.871 of the true one.
    const upper_airless = Angle.fromDegrees(33.246 / 60.0);
    const upper = observer.observedAltitude(upper_airless, .{});
    try testing.expectApproxEqAbs(24.618 + 0.0019279, arcmin(upper) - arcmin(upper_airless), 0.0005);
    const diameter = arcmin(upper) - 30;
    try testing.expectApproxEqAbs(27.864, diameter - 0.0019279, 0.0005);
    try testing.expectApproxEqAbs(0.871, diameter / 32, 0.0005);
}

test "refraction is never negative, is greatest at the horizon, and is held below −1°" {
    var previous: f64 = math.inf(f64);
    var h: f64 = -1;
    while (h <= 90) : (h += 0.05) {
        const r = arcmin(observer.observedAltitude(Angle.fromDegrees(h), .{})) - h * 60;
        try testing.expect(r >= 0);
        try testing.expect(r <= previous + 1e-9);
        previous = r;
    }
    // Zero at the zenith, to the formula's constant's rounding.
    try testing.expectApproxEqAbs(0, previous, 1e-6);
    // 38.80′ at −1°, and the same at any altitude below it.
    const held = arcmin(observer.observedAltitude(Angle.fromDegrees(-1), .{})) + 60;
    try testing.expectApproxEqAbs(38.80, held, 0.005);
    for ([_]f64{ -1.5, -5.11, -30, -90 }) |below| {
        const r = arcmin(observer.observedAltitude(Angle.fromDegrees(below), .{})) - below * 60;
        try testing.expectApproxEqAbs(held, r, 1e-9);
    }
}

test "the observed and airless altitudes both increase, without a jump" {
    // Steps of 0.001° in one move the other by more than nothing and, with
    // no jump anywhere, by no more than 0.0015°: near the horizon the
    // airless altitude moves 1.25 times as fast as the observed one,
    // because refraction falls off so quickly there.
    var h: f64 = -3;
    var last_observed = observer.observedAltitude(Angle.fromDegrees(h), .{}).degrees();
    var last_airless = observer.airlessAltitude(Angle.fromDegrees(h), .{}).degrees();
    while (h <= 90) {
        h += 0.001;
        const observed = observer.observedAltitude(Angle.fromDegrees(h), .{}).degrees();
        const back = observer.airlessAltitude(Angle.fromDegrees(h), .{}).degrees();
        try testing.expect(observed > last_observed and observed - last_observed < 0.0015);
        try testing.expect(back > last_airless and back - last_airless < 0.0015);
        last_observed = observed;
        last_airless = back;
    }
}

test "pressure and temperature scale the refraction" {
    const h = Angle.fromDegrees(10);
    const standard = arcmin(observer.observedAltitude(h, .{})) - 600;
    // Half the pressure, half the refraction; and at 10 °C colder, about
    // 3.6% more.
    const thin = arcmin(observer.observedAltitude(h, .{ .pressure_mb = 505 })) - 600;
    try testing.expectApproxEqAbs(standard / 2, thin, 1e-9);
    const cold = arcmin(observer.observedAltitude(h, .{ .temperature_c = 0 })) - 600;
    try testing.expectApproxEqAbs(standard * 283.0 / 273.0, cold, 1e-9);
}

// `airlessAltitude` inverts `observedAltitude` by iteration, in whatever
// air, from the nadir to the zenith, through the horizon and the altitude
// below which refraction is held.
test "there and back, in any air" {
    const airs = [_]observer.Atmosphere{
        .{},
        .{ .pressure_mb = 600, .temperature_c = -30 },
        .{ .pressure_mb = 1050, .temperature_c = -40 },
        .{ .pressure_mb = 0 },
        .{ .pressure_mb = 2020, .temperature_c = -100 },
    };
    var worst: f64 = 0;
    for (airs) |air| {
        var h: f64 = -90;
        while (h <= 90) : (h += 0.003) {
            const there = observer.observedAltitude(Angle.fromDegrees(h), air);
            const back = observer.airlessAltitude(there, air);
            worst = @max(worst, @abs(back.degrees() - h) * 3600);
        }
    }
    errdefer std.debug.print("worst: {e}″\n", .{worst});
    try testing.expect(worst < 1e-6);
}

const airless = horizons.parse(@embedFile("reference/refraction-airless.txt"));
const refracted = horizons.parse(@embedFile("reference/refraction-refracted.txt"));

// JPL Horizons: the Sun from Palomar every three minutes on 2026 June 21,
// airless and refracted, from 25° below the horizon to 80° above it (see
// tools/horizons.sh). Horizons refracts by Sæmundsson's formula with the
// zenith constant, in the default atmosphere whatever the site's height,
// held below −1° as here: `observedAltitude` reproduces its refracted
// elevations, and `airlessAltitude` its airless ones from them, to
// 0.00001″, the 0.000000001° it prints them to; held to 0.0001″.
test "refraction against JPL Horizons" {
    const jd = airless.column("Date");
    const airless_elevation = airless.column("Elevation");
    const refracted_jd = refracted.column("Date");
    const refracted_elevation = refracted.column("Elevation");

    var worst: [2]f64 = @splat(0);
    for (airless.rows[0..airless.row_count], refracted.rows[0..refracted.row_count]) |a, r| {
        try testing.expectEqual(a[jd], r[refracted_jd]);
        const observed = observer.observedAltitude(Angle.fromDegrees(a[airless_elevation]), .{}).degrees();
        worst[0] = @max(worst[0], @abs(observed - r[refracted_elevation]) * 3600);
        const back = observer.airlessAltitude(Angle.fromDegrees(r[refracted_elevation]), .{}).degrees();
        worst[1] = @max(worst[1], @abs(back - a[airless_elevation]) * 3600);
    }
    errdefer std.debug.print("worst: {d:.5}″ there, {d:.5}″ back\n", .{ worst[0], worst[1] });
    try testing.expect(worst[0] < 0.0001);
    try testing.expect(worst[1] < 0.0001);
}

test "refraction of a non-finite altitude is NaN" {
    for ([_]f64{ math.nan(f64), math.inf(f64), -math.inf(f64) }) |x| {
        try testing.expect(math.isNan(observer.observedAltitude(Angle.fromDegrees(x), .{}).radians()));
        try testing.expect(math.isNan(observer.airlessAltitude(Angle.fromDegrees(x), .{}).radians()));
    }
}