Something went wrong. Try again.
A Zig library for astronomical calculations: the Sun and Moon from anywhere on Earth, the Moon's phases and libration, and time scales.
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153// 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())); }}