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.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>// SPDX-License-Identifier: MIT
//! What the calculations must do with any instant at all.//!//! Every target is a property, not an example: whatever `f64` arrives — NaN,//! ±∞, a subnormal, ±1e308 — nothing panics, a non-finite instant gives NaN,//! and a finite one gives results in their documented ranges.//!//! Each target is an ordinary test as well as a fuzz target. Without `--fuzz`//! it runs the seeds beside it, so `zig build test` checks the same//! properties on inputs that have been interesting once. With `--fuzz` they//! are the targets of Zig's coverage-guided fuzzer: `zig build fuzz --fuzz`.//!//! A seed is the bytes `Smith.value` reads: eight little-endian bytes for each//! `f64`, then one byte for each enum.
const std = @import("std");const math = std.math;const testing = std.testing;const Smith = testing.Smith;
const astro = @import("astro");const time = astro.time;
/// The eight bytes `Smith.value(f64)` reads as `x`.fn seed(comptime x: f64) [8]u8 { return @bitCast(std.mem.nativeToLittle(u64, @as(u64, @bitCast(x))));}
const instant_seeds = [_][]const u8{ &seed(2451545.0), &seed(0.0), &seed(-0.0), &seed(2299160.5), &seed(-1e6), &seed(7.8e11), &seed(-7.8e11), &seed(7.9e11), &seed(1e308), &seed(-1e308), &seed(math.floatTrueMin(f64)), &seed(math.nan(f64)), &seed(math.inf(f64)), &seed(-math.inf(f64)), // Found by the fuzzer: the year 10530, where the terms in T of the // nutation take Δψ past 20″. &seed(5566261.295362834),};
fn fuzzInstant(_: void, smith: *Smith) anyerror!void { const jd = smith.value(f64); const scale = smith.value(time.TimeScale); try instantProperty(jd, scale);}
/// The model, the conversions and the calendars, from any Julian Day.fn instantProperty(jd: f64, scale: time.TimeScale) !void { const instant = time.Instant.fromJulianDay(jd, scale);
const delta_t = instant.deltaTSeconds(); const tt = instant.to(.tt, .model); const ut = instant.to(.ut, .model); const unix = instant.toUnixSeconds(.model); const gregorian = instant.toCalendar(.gregorian); const julian = instant.toCalendar(.julian);
if (!math.isFinite(jd)) { try testing.expect(math.isNan(delta_t)); try testing.expect(math.isNan(tt.julian_day)); try testing.expect(math.isNan(ut.julian_day)); try testing.expect(math.isNan(unix)); try testing.expect(gregorian == null); try testing.expect(julian == null); return; }
try testing.expectEqual(time.TimeScale.tt, tt.scale); try testing.expectEqual(time.TimeScale.ut, ut.scale);
// Within ±2 billion years or so every result is a real number, and ΔT is // never less than its minimum of about −20 s. if (@abs(jd) <= 7.8e11) { try testing.expect(math.isFinite(delta_t)); try testing.expect(delta_t > -25); try testing.expect(math.isFinite(tt.julian_day)); try testing.expect(math.isFinite(ut.julian_day));
for ([_]time.Calendar{ .gregorian, .julian }, [_]?time.CalendarDate{ gregorian, julian }) |calendar, maybe| { const date = maybe orelse return error.TestUnexpectedResult; try testing.expect(date.day >= 1 and date.day < 32); const back = time.Instant.fromCalendar(calendar, date, scale).julian_day; // Partly relative, because at 1e11 days an f64 resolves only // 1e-5 d. try testing.expectApproxEqAbs(jd, back, @max(1e-6, @abs(jd) * 1e-12)); } }}
test "fuzz: instants, ΔT and calendars" { for (instant_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzInstant({}, &smith); } try testing.fuzz({}, fuzzInstant, .{ .corpus = &instant_seeds });}
const calendar_seeds = [_][]const u8{ &(seed(1.0) ++ [_]u8{ 0, 0, 0, 0, 1, 0 }), &(seed(math.nan(f64)) ++ [_]u8{ 0xff, 0xff, 0xff, 0x7f, 12, 1 }), &(seed(-1e308) ++ [_]u8{ 0, 0, 0, 0x80, 2, 0 }),};
fn fuzzCalendar(_: void, smith: *Smith) anyerror!void { const day = smith.value(f64); const year = smith.value(i32); const month = smith.value(time.Month); const calendar = smith.value(time.Calendar); try calendarProperty(calendar, .{ .year = year, .month = month, .day = day });}
/// Any date, in either calendar, makes an instant without panicking, and a/// sensible one comes back as itself.fn calendarProperty(calendar: time.Calendar, date: time.CalendarDate) !void { const instant = time.Instant.fromCalendar(calendar, date, .ut); if (!math.isFinite(date.day)) { try testing.expect(!math.isFinite(instant.julian_day)); return; } if (date.day >= 1 and date.day < 29) { const back = instant.toCalendar(calendar) orelse return; try testing.expectEqual(date.year, back.year); try testing.expectEqual(date.month, back.month); try testing.expectApproxEqAbs(date.day, back.day, 1e-4); }}
test "fuzz: calendar dates" { for (calendar_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzCalendar({}, &smith); } try testing.fuzz({}, fuzzCalendar, .{ .corpus = &calendar_seeds });}
const angle_seeds = [_][]const u8{ &(seed(0.0) ++ seed(0.0) ++ seed(0.409)), &(seed(-1e-300) ++ seed(math.pi / 2.0) ++ seed(0.409)), &(seed(1e300) ++ seed(-math.pi / 2.0) ++ seed(-1e300)), &(seed(math.nan(f64)) ++ seed(1.0) ++ seed(0.409)), &(seed(2.0) ++ seed(math.inf(f64)) ++ seed(0.409)), &(seed(7.0) ++ seed(0.3) ++ seed(math.nan(f64))),};
fn fuzzAngles(_: void, smith: *Smith) anyerror!void { const a = smith.value(f64); const b = smith.value(f64); const e = smith.value(f64); try angleProperty(a, b, e);}
/// Normalization lands in range for any finite angle and gives NaN for any/// other; a conversion between frames never panics, and from a real/// direction gives a real direction in range.fn angleProperty(a: f64, b: f64, e: f64) !void { const angle = astro.Angle.fromRadians(a); const n = angle.normalized().radians(); const s = angle.signed().radians(); if (math.isFinite(a)) { try testing.expect(n >= 0 and n < math.tau); try testing.expect(s > -math.pi and s <= math.pi); } else { try testing.expect(math.isNan(n) and math.isNan(s)); }
const ecliptic: astro.coord.Ecliptic = .{ .longitude = angle, .latitude = astro.Angle.fromRadians(b) }; const equatorial: astro.coord.Equatorial = .{ .right_ascension = angle, .declination = astro.Angle.fromRadians(b) }; const obliquity = astro.Angle.fromRadians(e); const to_equatorial = ecliptic.toEquatorial(obliquity); const to_ecliptic = equatorial.toEcliptic(obliquity);
// A real direction: finite, latitude within the poles, and not so large // that its sine and cosine have lost every digit. const real = math.isFinite(a) and math.isFinite(b) and math.isFinite(e) and @abs(a) < 1e6 and @abs(b) <= math.pi / 2.0 and @abs(e) < 1e6; if (real) { for ([_]f64{ to_equatorial.right_ascension.radians(), to_ecliptic.longitude.radians() }) |x| { try testing.expect(x >= 0 and x < math.tau); } for ([_]f64{ to_equatorial.declination.radians(), to_ecliptic.latitude.radians() }) |x| { try testing.expect(x >= -math.pi / 2.0 and x <= math.pi / 2.0); } }}
test "fuzz: angles and coordinate conversions" { for (angle_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzAngles({}, &smith); } try testing.fuzz({}, fuzzAngles, .{ .corpus = &angle_seeds });}
fn fuzzEarth(_: void, smith: *Smith) anyerror!void { const jd = smith.value(f64); const scale = smith.value(time.TimeScale); try earthProperty(jd, scale);}
/// Nutation and obliquity from any instant: NaN from a non-finite one, and/// within a few thousand years of the present the values the theory allows.fn earthProperty(jd: f64, scale: time.TimeScale) !void { const instant = time.Instant.fromJulianDay(jd, scale); const n = astro.earth.nutation(instant); const mean = astro.earth.meanObliquity(instant); const true_ = astro.earth.trueObliquity(instant); if (!math.isFinite(jd)) { try testing.expect(math.isNan(n.longitude.radians())); try testing.expect(math.isNan(n.obliquity.radians())); try testing.expect(math.isNan(mean.radians())); try testing.expect(math.isNan(true_.radians())); return; } // Within Laskar's ±10,000 years. No sum of Table 22.A's terms can exceed // the sum of their magnitudes, which grows with T: 19.4316″ + 0.01823″|T| // for Δψ and 10.0733″ + 0.00141″|T| for Δε. The obliquity stays between // Laskar's extremes of 22°36′41″ and 24°14′07″ (Meeus, p. 147). if (@abs(jd - time.Instant.j2000_julian_day) < 3.6e6) { const t = @abs(instant.julianCenturies()); try testing.expect(@abs(n.longitude.arcseconds()) <= 19.4316 + 0.01823 * t); try testing.expect(@abs(n.obliquity.arcseconds()) <= 10.0733 + 0.00141 * t); try testing.expect(mean.degrees() > 22.6 and mean.degrees() < 24.25); }}
test "fuzz: nutation and obliquity" { for (instant_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzEarth({}, &smith); } try testing.fuzz({}, fuzzEarth, .{ .corpus = &instant_seeds });}
fn fuzzSidereal(_: void, smith: *Smith) anyerror!void { const jd = smith.value(f64); const scale = smith.value(time.TimeScale); try siderealProperty(jd, scale);}
/// Sidereal time from any instant: NaN from a non-finite one; otherwise in/// [0°, 360°) or, where the polynomial in T overflows, NaN; and within/// Laskar's ±10,000 years, apparent and mean differ by no more than the/// nutation in longitude.fn siderealProperty(jd: f64, scale: time.TimeScale) !void { const instant = time.Instant.fromJulianDay(jd, scale); const mean = astro.earth.greenwichMeanSiderealTime(instant).degrees(); const apparent = astro.earth.greenwichApparentSiderealTime(instant).degrees(); if (!math.isFinite(jd)) { try testing.expect(math.isNan(mean)); try testing.expect(math.isNan(apparent)); return; } for ([_]f64{ mean, apparent }) |theta| { try testing.expect(math.isNan(theta) or (theta >= 0 and theta < 360)); } if (@abs(jd - time.Instant.j2000_julian_day) < 3.6e6) { try testing.expect(!math.isNan(mean) and !math.isNan(apparent)); const dpsi = @abs(astro.earth.nutation(instant).longitude.degrees()); const difference = @abs(astro.Angle.fromDegrees(apparent - mean).signed().degrees()); try testing.expect(difference <= dpsi + 1e-9); }}
test "fuzz: sidereal time" { for (instant_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzSidereal({}, &smith); } try testing.fuzz({}, fuzzSidereal, .{ .corpus = &instant_seeds });}
fn fuzzSun(_: void, smith: *Smith) anyerror!void { const jd = smith.value(f64); const scale = smith.value(time.TimeScale); try sunProperty(jd, scale);}
/// The Sun from any instant: NaN from a non-finite one, and within the/// ±10,000 years the obliquity is good for, coordinates in their ranges and/// a distance the Earth's orbit allows.fn sunProperty(jd: f64, scale: time.TimeScale) !void { const p = astro.sun.position(time.Instant.fromJulianDay(jd, scale)); if (!math.isFinite(jd)) { try testing.expect(math.isNan(p.ecliptic.longitude.radians())); try testing.expect(math.isNan(p.equatorial.right_ascension.radians())); try testing.expect(math.isNan(p.equatorial.declination.radians())); try testing.expect(math.isNan(p.distance_au)); return; } if (@abs(jd - time.Instant.j2000_julian_day) < 3.6e6) { const lon = p.ecliptic.longitude.radians(); const ra = p.equatorial.right_ascension.radians(); try testing.expect(lon >= 0 and lon < math.tau); try testing.expect(ra >= 0 and ra < math.tau); // The Sun is never further from the equator than the obliquity. try testing.expect(@abs(p.equatorial.declination.degrees()) < 24.3); // Perihelion and aphelion, with room for the eccentricity's drift. try testing.expect(p.distance_au > 0.97 and p.distance_au < 1.03); }}
test "fuzz: the Sun" { for (instant_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzSun({}, &smith); } try testing.fuzz({}, fuzzSun, .{ .corpus = &instant_seeds });}
fn fuzzMoon(_: void, smith: *Smith) anyerror!void { const jd = smith.value(f64); const scale = smith.value(time.TimeScale); try moonProperty(jd, scale);}
/// The Moon from any instant: NaN from a non-finite one, and within a few/// thousand years of the present, coordinates in their ranges and a distance/// between the extremes of perigee and apogee.fn moonProperty(jd: f64, scale: time.TimeScale) !void { const p = astro.moon.position(time.Instant.fromJulianDay(jd, scale)); if (!math.isFinite(jd)) { try testing.expect(math.isNan(p.ecliptic.longitude.radians())); try testing.expect(math.isNan(p.ecliptic.latitude.radians())); try testing.expect(math.isNan(p.equatorial.right_ascension.radians())); try testing.expect(math.isNan(p.equatorial.declination.radians())); try testing.expect(math.isNan(p.distance_km)); return; } if (@abs(jd - time.Instant.j2000_julian_day) < 1.1e6) { const lon = p.ecliptic.longitude.radians(); try testing.expect(lon >= 0 and lon < math.tau); // The Moon's orbit is inclined about 5.1° to the ecliptic. try testing.expect(@abs(p.ecliptic.latitude.degrees()) < 5.5); try testing.expect(@abs(p.equatorial.declination.degrees()) < 29); // Perigee is never nearer than about 356,400 km nor apogee further // than about 406,700 km. try testing.expect(p.distance_km > 356000 and p.distance_km < 407000); }}
test "fuzz: the Moon" { for (instant_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzMoon({}, &smith); } try testing.fuzz({}, fuzzMoon, .{ .corpus = &instant_seeds });}
/// An instant, then an observer's latitude and longitude in degrees and/// height in meters, then the time scale.const topocentric_seeds = [_][]const u8{ &(seed(2451545.0) ++ seed(33.3561) ++ seed(-116.8625) ++ seed(1706) ++ [_]u8{0}), // The poles, where ρ cos φ′ is zero and the hour angle means nothing. &(seed(2460000.5) ++ seed(90) ++ seed(0) ++ seed(0) ++ [_]u8{1}), &(seed(2460000.5) ++ seed(-90) ++ seed(45) ++ seed(0) ++ [_]u8{0}), // The antimeridian, from both sides, and far past it. &(seed(2455000.25) ++ seed(-20) ++ seed(180) ++ seed(0) ++ [_]u8{0}), &(seed(2455000.25) ++ seed(-20) ++ seed(-180) ++ seed(0) ++ [_]u8{1}), &(seed(2455000.25) ++ seed(10) ++ seed(1e300) ++ seed(0) ++ [_]u8{0}), // Below sea level, and further from the Earth than the Moon is. &(seed(2440000.5) ++ seed(31.5) ++ seed(35.5) ++ seed(-430) ++ [_]u8{0}), &(seed(2440000.5) ++ seed(0) ++ seed(0) ++ seed(5e8) ++ [_]u8{1}), // A latitude just past the pole, and things that are not numbers. &(seed(2451545.0) ++ seed(90.000001) ++ seed(0) ++ seed(0) ++ [_]u8{0}), &(seed(2451545.0) ++ seed(math.nan(f64)) ++ seed(0) ++ seed(0) ++ [_]u8{0}), &(seed(2451545.0) ++ seed(0) ++ seed(math.inf(f64)) ++ seed(0) ++ [_]u8{0}), &(seed(2451545.0) ++ seed(0) ++ seed(0) ++ seed(-math.inf(f64)) ++ [_]u8{0}), &(seed(math.nan(f64)) ++ seed(0) ++ seed(0) ++ seed(0) ++ [_]u8{0}), &(seed(1e308) ++ seed(45) ++ seed(45) ++ seed(1e308) ++ [_]u8{1}),};
fn fuzzTopocentricSun(_: void, smith: *Smith) anyerror!void { try topocentricProperty(.sun, smith);}
fn fuzzTopocentricMoon(_: void, smith: *Smith) anyerror!void { try topocentricProperty(.moon, smith);}
/// The Sun or the Moon from anywhere at all: NaN from a non-finite instant/// or an observer that is not valid; otherwise every coordinate NaN or in/// its range; and within ±10,000 years, from no further than 10,000 km/// above or below the ellipsoid, the shift from geocentric is within the/// parallax and the distance changes by no more than the observer's own/// distance from the Earth's center.fn topocentricProperty(comptime body: enum { sun, moon }, smith: *Smith) !void { const jd = smith.value(f64); const observer: astro.observer.Observer = .{ .latitude = astro.Angle.fromDegrees(smith.value(f64)), .longitude = astro.Angle.fromDegrees(smith.value(f64)), .height_m = smith.value(f64), }; const scale = smith.value(time.TimeScale); const instant = time.Instant.fromJulianDay(jd, scale); const t, const geocentric, const distance_km = switch (body) { .sun => r: { const p = astro.sun.position(instant); break :r .{ astro.sun.topocentric(instant, observer), p.equatorial, p.distance_au * 149597870.7 }; }, .moon => r: { const p = astro.moon.position(instant); break :r .{ astro.moon.topocentric(instant, observer), p.equatorial, p.distance_km }; }, }; const values = [_]f64{ t.equatorial.right_ascension.degrees(), t.equatorial.declination.degrees(), t.ecliptic.longitude.degrees(), t.ecliptic.latitude.degrees(), t.hour_angle.degrees(), t.distance_km, t.horizontal.azimuth.degrees(), t.horizontal.altitude.degrees(), }; const valid = @abs(observer.latitude.degrees()) <= 90 and math.isFinite(observer.longitude.degrees()) and math.isFinite(observer.height_m); if (!math.isFinite(jd) or !valid) { for (values) |v| try testing.expect(math.isNan(v)); return; } const ra, const dec, const lon, const lat, const hour_angle, const distance, const azimuth, const altitude = values; try testing.expect(math.isNan(ra) or (ra >= 0 and ra < 360)); try testing.expect(math.isNan(dec) or @abs(dec) <= 90); try testing.expect(math.isNan(lon) or (lon >= 0 and lon < 360)); try testing.expect(math.isNan(lat) or @abs(lat) <= 90); try testing.expect(math.isNan(hour_angle) or (hour_angle > -180 and hour_angle <= 180)); try testing.expect(math.isNan(distance) or distance >= 0); try testing.expect(math.isNan(azimuth) or (azimuth >= 0 and azimuth < 360)); try testing.expect(math.isNan(altitude) or @abs(altitude) <= 90);
if (@abs(jd - time.Instant.j2000_julian_day) < 3.6e6 and @abs(observer.height_m) < 1e7) { const g = observer.geocentric(); const rho_km = @sqrt(g.rho_sin_phi * g.rho_sin_phi + g.rho_cos_phi * g.rho_cos_phi) * astro.observer.equatorial_radius_km; try testing.expect(@abs(distance - distance_km) <= rho_km * (1 + 1e-9) + 1e-6); const d1 = geocentric.declination.radians(); const d2 = t.equatorial.declination.radians(); const cos_shift = @sin(d1) * @sin(d2) + @cos(d1) * @cos(d2) * @cos(geocentric.right_ascension.radians() - t.equatorial.right_ascension.radians()); const shift = math.acos(@min(1.0, @max(-1.0, cos_shift))); try testing.expect(shift <= math.asin(rho_km / distance_km) + 1e-9); }}
test "fuzz: the topocentric Sun" { for (topocentric_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzTopocentricSun({}, &smith); } try testing.fuzz({}, fuzzTopocentricSun, .{ .corpus = &topocentric_seeds });}
test "fuzz: the topocentric Moon" { for (topocentric_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzTopocentricMoon({}, &smith); } try testing.fuzz({}, fuzzTopocentricMoon, .{ .corpus = &topocentric_seeds });}
/// A declination, an hour angle and a latitude, in degrees.const horizontal_seeds = [_][]const u8{ &(seed(-6.72) ++ seed(64.35) ++ seed(38.92)), // At the zenith and the nadir, where azimuth means nothing. &(seed(38.92) ++ seed(0) ++ seed(38.92)), &(seed(-38.92) ++ seed(180) ++ seed(38.92)), // At the poles of the sky, and from the Earth's poles. &(seed(90) ++ seed(123) ++ seed(45)), &(seed(-90) ++ seed(-45) ++ seed(-30)), &(seed(20) ++ seed(77) ++ seed(90)), &(seed(20) ++ seed(77) ++ seed(-90)), // On the horizon, and past every range. &(seed(0) ++ seed(-90) ++ seed(0)), &(seed(1e300) ++ seed(-1e300) ++ seed(1e300)), &(seed(math.nan(f64)) ++ seed(0) ++ seed(0)), &(seed(0) ++ seed(math.inf(f64)) ++ seed(0)),};
fn fuzzHorizontal(_: void, smith: *Smith) anyerror!void { const dec = smith.value(f64); const hour_angle = smith.value(f64); const latitude = smith.value(f64); try horizontalProperty(dec, hour_angle, latitude);}
/// Horizontal coordinates from any direction, hour angle and latitude:/// NaN from anything not finite; otherwise azimuth in [0°, 360°) and/// altitude in [−90°, +90°]; and, for a declination and latitude in their/// ranges, the altitude of eq. 13.6.fn horizontalProperty(dec: f64, hour_angle: f64, latitude: f64) !void { const body: astro.coord.Equatorial = .{ .right_ascension = astro.Angle.fromDegrees(0), .declination = astro.Angle.fromDegrees(dec), }; const h = body.toHorizontal(astro.Angle.fromDegrees(hour_angle), astro.Angle.fromDegrees(latitude)); const azimuth = h.azimuth.degrees(); const altitude = h.altitude.degrees(); if (!math.isFinite(dec) or !math.isFinite(hour_angle) or !math.isFinite(latitude)) { try testing.expect(math.isNan(azimuth)); try testing.expect(math.isNan(altitude)); return; } try testing.expect(azimuth >= 0 and azimuth < 360); try testing.expect(@abs(altitude) <= 90); if (@abs(dec) <= 90 and @abs(latitude) <= 90) { const d = dec * math.rad_per_deg; const phi = latitude * math.rad_per_deg; const sin_altitude = @sin(phi) * @sin(d) + @cos(phi) * @cos(d) * @cos(hour_angle * math.rad_per_deg); try testing.expectApproxEqAbs(sin_altitude, @sin(altitude * math.rad_per_deg), 1e-12); }}
test "fuzz: horizontal coordinates" { for (horizontal_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzHorizontal({}, &smith); } try testing.fuzz({}, fuzzHorizontal, .{ .corpus = &horizontal_seeds });}
/// An altitude in degrees, a pressure in millibars and a temperature in/// degrees Celsius.const refraction_seeds = [_][]const u8{ &(seed(0.5) ++ seed(1010) ++ seed(10)), // Either side of where refraction is held, and the poles of the // formulas, at −5.11° and −4.4°. &(seed(-1) ++ seed(1010) ++ seed(10)), &(seed(-0.3115) ++ seed(1010) ++ seed(10)), &(seed(-5.11) ++ seed(1010) ++ seed(10)), &(seed(-4.4) ++ seed(1010) ++ seed(10)), // The zenith and the nadir, a mountaintop in winter, and no air. &(seed(90) ++ seed(1010) ++ seed(10)), &(seed(-90) ++ seed(1010) ++ seed(10)), &(seed(2) ++ seed(600) ++ seed(-30)), &(seed(2) ++ seed(0) ++ seed(10)), // Absolute zero, and things that are not altitudes at all. &(seed(2) ++ seed(1010) ++ seed(-273)), &(seed(1e300) ++ seed(1010) ++ seed(10)), &(seed(math.nan(f64)) ++ seed(1010) ++ seed(10)), &(seed(math.inf(f64)) ++ seed(1010) ++ seed(10)),};
fn fuzzRefraction(_: void, smith: *Smith) anyerror!void { const altitude = smith.value(f64); const atmosphere: astro.observer.Atmosphere = .{ .pressure_mb = smith.value(f64), .temperature_c = smith.value(f64), }; try refractionProperty(altitude, atmosphere);}
/// Refraction of any altitude in any air: NaN from a non-finite altitude;/// and for an altitude in [−90°, +90°], air at a pressure between none and/// twice the standard and a temperature between −100 °C and +60 °C, a/// refraction between nothing and the 38.80′ held below the horizon,/// scaled for the air, either way, and there and back to a millionth of a/// second of arc.fn refractionProperty(altitude: f64, atmosphere: astro.observer.Atmosphere) !void { const h = astro.Angle.fromDegrees(altitude); const observed = astro.observer.observedAltitude(h, atmosphere).degrees(); const airless = astro.observer.airlessAltitude(h, atmosphere).degrees(); if (!math.isFinite(altitude)) { try testing.expect(math.isNan(observed)); try testing.expect(math.isNan(airless)); return; } // The altitude as the angle holds it: degrees to radians and back can // move the last bit, and with next to no air that is more than the // refraction. const a = h.degrees(); const p = atmosphere.pressure_mb; const t = atmosphere.temperature_c; if (@abs(a) <= 90 and p >= 0 and p <= 2020 and t >= -100 and t <= 60) { const most = 38.80 / 60.0 * p / 1010 * 283 / (273 + t) + 1e-9; try testing.expect(observed >= a and observed - a <= most); try testing.expect(airless <= a and a - airless <= most); const back = astro.observer.airlessAltitude(astro.Angle.fromDegrees(observed), atmosphere).degrees(); try testing.expect(@abs(back - a) * 3600 < 1e-6); }}
test "fuzz: refraction" { for (refraction_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzRefraction({}, &smith); } try testing.fuzz({}, fuzzRefraction, .{ .corpus = &refraction_seeds });}
fn fuzzPhase(_: void, smith: *Smith) anyerror!void { const jd = smith.value(f64); const scale = smith.value(time.TimeScale); try phaseProperty(jd, scale);}
/// The phase from any instant: NaN and no name from a non-finite one, and/// within a few thousand years of the present every quantity in its range,/// the illuminated fraction agreeing with the phase angle, and the name/// agreeing with the elongation in longitude.fn phaseProperty(jd: f64, scale: time.TimeScale) !void { const p = astro.moon.phase.at(time.Instant.fromJulianDay(jd, scale)); if (!math.isFinite(jd)) { try testing.expect(math.isNan(p.illuminated_fraction)); try testing.expect(p.named() == null and p.waxing() == null); return; } if (@abs(jd - time.Instant.j2000_julian_day) < 1.1e6) { const psi = p.elongation.degrees(); const i = p.phase_angle.degrees(); const l = p.elongation_in_longitude.degrees(); const chi = p.bright_limb.degrees(); try testing.expect(psi >= 0 and psi <= 180); try testing.expect(i >= 0 and i <= 180); try testing.expect(l >= 0 and l < 360); try testing.expect(chi >= 0 and chi < 360); try testing.expect(p.illuminated_fraction >= 0 and p.illuminated_fraction <= 1); try testing.expectApproxEqAbs((1 + @cos(p.phase_angle.radians())) / 2, p.illuminated_fraction, 1e-12); const named = p.named() orelse return error.TestUnexpectedResult; const center = @as(f64, @floatFromInt(@backingInt(named))) * 45; try testing.expect(@abs(astro.Angle.fromDegrees(l - center).signed().degrees()) <= 22.5); }}
test "fuzz: the Moon's phase" { for (instant_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzPhase({}, &smith); } try testing.fuzz({}, fuzzPhase, .{ .corpus = &instant_seeds });}
const phases_seeds = [_][]const u8{ &(seed(2451545.0) ++ seed(2451600.0) ++ [_]u8{ 0, 2 }), &(seed(1e300) ++ seed(-1e300) ++ [_]u8{ 1, 1 }), &(seed(-2.6e16) ++ seed(2.6e16) ++ [_]u8{ 0, 3 }), &(seed(math.nan(f64)) ++ seed(math.inf(f64)) ++ [_]u8{ 1, 0 }), &(seed(2451550.09766) ++ seed(math.inf(f64)) ++ [_]u8{ 0, 0 }), // Found by the fuzzer: 350 million years on, where the series gave // phases out of order. &(seed(130254491473.4171) ++ seed(1.1717362069003993e+209) ++ [_]u8{ 0, 0 }), // Either side of the last lunation computed. &(seed(2451545.0 + 999.99 * 36525.0) ++ seed(math.inf(f64)) ++ [_]u8{ 0, 1 }),};
fn fuzzPhases(_: void, smith: *Smith) anyerror!void { const a = smith.value(f64); const b = smith.value(f64); const scale = smith.value(time.TimeScale); const kind = smith.value(astro.moon.phase.PrincipalPhase); try phasesProperty(a, b, scale, kind);}
/// Principal phases and lunar age from any instant: never a panic — the/// float-to-integer steps are the danger — NaN from a non-finite instant,/// and near the present the next phase after, the previous one before, and/// an age in range. Iteration between any two instants yields phases in/// order and inside the interval.fn phasesProperty(a: f64, b: f64, scale: time.TimeScale, kind: astro.moon.phase.PrincipalPhase) !void { const phase = astro.moon.phase; const instant = time.Instant.fromJulianDay(a, scale); const next = phase.next(instant, kind).julian_day; const previous = phase.previous(instant, kind).julian_day; const age = phase.lunarAgeDays(instant); if (!math.isFinite(a)) { try testing.expect(math.isNan(next) and math.isNan(previous) and math.isNan(age)); } else if (@abs(a - time.Instant.j2000_julian_day) < 1.1e6) { const jd = instant.to(.tt, .model).julian_day; try testing.expect(next > jd and previous < jd); try testing.expect(age >= 0 and age < 30); }
// At most a few hundred steps, so that a vast interval cannot hang the // fuzzer; what is checked is order and bounds, not completeness. var it = phase.between(instant, time.Instant.fromJulianDay(b, scale)); var last = -math.inf(f64); var steps: usize = 0; while (it.next()) |event| : (steps += 1) { if (steps == 300) break; try testing.expect(event.instant.julian_day > last); last = event.instant.julian_day; }}
test "fuzz: principal phases and lunar age" { for (phases_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzPhases({}, &smith); } try testing.fuzz({}, fuzzPhases, .{ .corpus = &phases_seeds });}
fn fuzzLibration(_: void, smith: *Smith) anyerror!void { const jd = smith.value(f64); const scale = smith.value(time.TimeScale); try librationProperty(jd, scale);}
/// The libration from any instant: NaN from a non-finite one, and within a/// few thousand years of the present the ranges the geometry allows — the/// optical libration never more than about 8° in longitude and 7° in/// latitude, and the axis within 25.4° of north (25.36° is the most seen in/// ±3,000 years). The physical libration is held to 0.065°: Meeus says it/// never exceeds 0.04°, but his series reaches 0.061° within 1900–2100,/// which the sums of its coefficients allow.fn librationProperty(jd: f64, scale: time.TimeScale) !void { const l = astro.moon.libration.at(time.Instant.fromJulianDay(jd, scale)); if (!math.isFinite(jd)) { try testing.expect(math.isNan(l.total.longitude.radians())); try testing.expect(math.isNan(l.total.latitude.radians())); try testing.expect(math.isNan(l.position_angle_of_axis.radians())); return; } if (@abs(jd - time.Instant.j2000_julian_day) < 1.1e6) { try testing.expect(@abs(l.optical.longitude.degrees()) < 8.5); try testing.expect(@abs(l.optical.latitude.degrees()) < 7.5); try testing.expect(@abs(l.physical.longitude.degrees()) < 0.065); try testing.expect(@abs(l.physical.latitude.degrees()) < 0.065); const total = l.total.longitude.degrees(); try testing.expect(total > -180 and total <= 180); const p = l.position_angle_of_axis.degrees(); try testing.expect((p >= 0 and p <= 25.4) or (p >= 334.6 and p < 360)); }}
test "fuzz: the Moon's libration" { for (instant_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzLibration({}, &smith); } try testing.fuzz({}, fuzzLibration, .{ .corpus = &instant_seeds });}
fn fuzzTopocentricLibration(_: void, smith: *Smith) anyerror!void { const jd = smith.value(f64); const observer: astro.observer.Observer = .{ .latitude = astro.Angle.fromDegrees(smith.value(f64)), .longitude = astro.Angle.fromDegrees(smith.value(f64)), .height_m = smith.value(f64), }; const scale = smith.value(time.TimeScale); try topocentricLibrationProperty(jd, observer, scale);}
/// The topocentric libration from any instant and any observer: NaN from a/// non-finite instant or an observer that is not valid; and within a few/// thousand years of the present, from within 10 km of the ellipsoid, no/// further from the geocentric libration than the Moon's parallax from/// there, and so within the geocentric ranges widened by about a degree.fn topocentricLibrationProperty(jd: f64, observer: astro.observer.Observer, scale: time.TimeScale) !void { const instant = time.Instant.fromJulianDay(jd, scale); const l = astro.moon.libration.topocentric(instant, observer); const valid = @abs(observer.latitude.degrees()) <= 90 and math.isFinite(observer.longitude.degrees()) and math.isFinite(observer.height_m); if (!math.isFinite(jd) or !valid) { try testing.expect(math.isNan(l.total.longitude.radians())); try testing.expect(math.isNan(l.total.latitude.radians())); try testing.expect(math.isNan(l.position_angle_of_axis.radians())); return; } if (@abs(jd - time.Instant.j2000_julian_day) < 1.1e6 and @abs(observer.height_m) < 10_000) { const g = astro.moon.libration.at(instant); const distance_km = astro.moon.position(instant).distance_km; const geocentric = observer.geocentric(); const rho_km = @sqrt(geocentric.rho_sin_phi * geocentric.rho_sin_phi + geocentric.rho_cos_phi * geocentric.rho_cos_phi) * astro.observer.equatorial_radius_km; const b1 = g.total.latitude.radians(); const b2 = l.total.latitude.radians(); const cos_moved = @sin(b1) * @sin(b2) + @cos(b1) * @cos(b2) * @cos(l.total.longitude.radians() - g.total.longitude.radians()); const moved = math.acos(@min(1.0, @max(-1.0, cos_moved))); try testing.expect(moved <= math.asin(rho_km / distance_km) + 0.001 * math.rad_per_deg);
try testing.expect(@abs(l.optical.longitude.degrees()) < 9.6); try testing.expect(@abs(l.optical.latitude.degrees()) < 8.6); try testing.expect(@abs(l.physical.longitude.degrees()) < 0.07); try testing.expect(@abs(l.physical.latitude.degrees()) < 0.07); const total = l.total.longitude.degrees(); try testing.expect(total > -180 and total <= 180); const p = l.position_angle_of_axis.degrees(); try testing.expect((p >= 0 and p <= 26.5) or (p >= 333.5 and p < 360)); }}
test "fuzz: the Moon's topocentric libration" { for (topocentric_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzTopocentricLibration({}, &smith); } try testing.fuzz({}, fuzzTopocentricLibration, .{ .corpus = &topocentric_seeds });}
fn fuzzSubsolar(_: void, smith: *Smith) anyerror!void { const jd = smith.value(f64); const scale = smith.value(time.TimeScale); try subsolarProperty(jd, scale);}
/// The subsolar point and colongitude from any instant: NaN from a/// non-finite one, and within a few thousand years of the present the/// longitude in (−180°, +180°], the latitude within the Moon's 1.5°/// obliquity to the ecliptic and a little libration, and the colongitude in/// [0°, 360°) and equal to 90° minus the longitude.fn subsolarProperty(jd: f64, scale: time.TimeScale) !void { const instant = time.Instant.fromJulianDay(jd, scale); const p = astro.moon.libration.subsolarPoint(instant); const c = astro.moon.libration.colongitude(instant); if (!math.isFinite(jd)) { try testing.expect(math.isNan(p.longitude.radians())); try testing.expect(math.isNan(p.latitude.radians())); try testing.expect(math.isNan(c.radians())); return; } if (@abs(jd - time.Instant.j2000_julian_day) < 1.1e6) { const lon = p.longitude.degrees(); try testing.expect(lon > -180 and lon <= 180); try testing.expect(@abs(p.latitude.degrees()) < 1.7); try testing.expect(c.degrees() >= 0 and c.degrees() < 360); try testing.expectApproxEqAbs(0, c.sub(astro.Angle.fromDegrees(90 - lon)).signed().degrees(), 1e-9); }}
test "fuzz: the subsolar point and colongitude" { for (instant_seeds) |s| { var smith: Smith = .{ .in = s }; try fuzzSubsolar({}, &smith); } try testing.fuzz({}, fuzzSubsolar, .{ .corpus = &instant_seeds });}