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.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>// SPDX-License-Identifier: MIT
//! `astro.observer` and the topocentric Sun and Moon, against Meeus's//! worked examples and JPL Horizons.
const std = @import("std");const math = std.math;const testing = std.testing;
const astro = @import("astro");const Angle = astro.Angle;const Instant = astro.time.Instant;const Observer = astro.observer.Observer;const horizons = @import("horizons.zig");
// Palomar Observatory as Meeus gives it in examples 11.a and 40.a: latitude// +33°21′22″, height 1706 m, and longitude 7h47m27s, which he counts west.const palomar: Observer = .{ .latitude = Angle.fromDms(.{ .degrees = 33, .minutes = 21, .seconds = 22 }), .longitude = Angle.fromHms(.{ .hours = 7, .minutes = 47, .seconds = 27 }).negate(), .height_m = 1706,};
test "example 11.a: ρ sin φ′ and ρ cos φ′ for Palomar" { // Printed to 0.000001, from the IAU 1976 ellipsoid; WGS 84's differs // from it by 3 m in the equatorial radius and moves neither digit. const g = palomar.geocentric(); try testing.expectApproxEqAbs(0.546861, g.rho_sin_phi, 0.0000005); try testing.expectApproxEqAbs(0.836339, g.rho_cos_phi, 0.0000005);}
// Example 40.a (pp. 280–281): Mars from Palomar on 2003 August 28 at// 3h17m00s UT, from the geocentric position the book gives, at 0.37276 au.test "example 40.a: Mars from Palomar" { const instant = Instant.fromCalendar(.gregorian, .init(2003, .august, 28, 3, 17, 0), .ut); // The book has θ₀ = 1h40m45s, to the second. const theta = Angle.fromHms(.{ .hours = 1, .minutes = 40, .seconds = 45 }); try testing.expectApproxEqAbs(theta.hours() * 3600, astro.earth.greenwichApparentSiderealTime(instant).hours() * 3600, 0.5);
const t = astro.observer.topocentric(instant, palomar, .{ .right_ascension = Angle.fromDegrees(339.530208), .declination = Angle.fromDegrees(-15.771083), }, 0.37276 * km_per_au); // α′ = 22h38m08.54s and δ′ = −15°46′30.0″, printed to 0.01 s and 0.1″. // This comes out at 08.545 s: the book rounds θ₀ to the second, and // takes π = 8.794″/Δ, from the IAU 1976 radius, where this takes the // ratio of WGS 84's radius to the distance. const alpha = Angle.fromHms(.{ .hours = 22, .minutes = 38, .seconds = 8.54 }); const delta = Angle.fromDms(.{ .negative = true, .degrees = 15, .minutes = 46, .seconds = 30.0 }); try testing.expectApproxEqAbs(alpha.hours() * 3600, t.equatorial.right_ascension.hours() * 3600, 0.01); try testing.expectApproxEqAbs(delta.arcseconds(), t.equatorial.declination.arcseconds(), 0.05);}
const km_per_au = 149597870.7;
const Site = struct { name: []const u8, observer: Observer, sun: horizons.Table, moon: horizons.Table,};
// The sites of tools/horizons.sh, as it gives them to Horizons.const sites = [_]Site{ .{ .name = "Palomar", .observer = .{ .latitude = Angle.fromDegrees(33.3561111), .longitude = Angle.fromDegrees(243.1375), .height_m = 1706 }, .sun = horizons.parse(@embedFile("reference/sun-palomar.txt")), .moon = horizons.parse(@embedFile("reference/moon-palomar.txt")), }, .{ .name = "Tromsø", .observer = .{ .latitude = Angle.fromDegrees(69.6492), .longitude = Angle.fromDegrees(18.9553), .height_m = 100 }, .sun = horizons.parse(@embedFile("reference/sun-tromso.txt")), .moon = horizons.parse(@embedFile("reference/moon-tromso.txt")), }, .{ .name = "Paranal", .observer = .{ .latitude = Angle.fromDegrees(-24.6272), .longitude = Angle.fromDegrees(289.5973), .height_m = 2635 }, .sun = horizons.parse(@embedFile("reference/sun-paranal.txt")), .moon = horizons.parse(@embedFile("reference/moon-paranal.txt")), },};
const geocentric_sun = horizons.parse(@embedFile("reference/sun.txt"));const geocentric_moon = horizons.parse(@embedFile("reference/moon.txt"));
/// The difference between two directions, in seconds of arc, however they/// wrap.fn separation(a: f64, b: f64) f64 { return @abs(Angle.fromDegrees(a - b).signed().arcseconds());}
/// The worst differences from Horizons, as right ascension times cos δ,/// declination, ecliptic longitude and ecliptic latitude in seconds of arc,/// distance in km, and then, for the position alone, azimuth times the/// cosine of the altitude and airless altitude in seconds of arc: of the/// shift from/// geocentric to topocentric, of the topocentric position from the generic/// conversion, and of the topocentric position from a body's own/// `topocentric` up to the last measured ΔT.const Worst = struct { shift: [5]f64 = @splat(0), position: [7]f64 = @splat(0), wrapper: [7]f64 = @splat(0),};
const Body = enum { sun, moon };
/// The last instant this library's ΔT is measured rather than predicted,/// 2026 October 1.const last_measured_delta_t = 2461314.5;
fn measure(topocentric: *const horizons.Table, geocentric: *const horizons.Table, observer: Observer, body: Body, worst: *Worst) void { const jd = topocentric.column("Date"); const ra = topocentric.column("R.A."); const dec = topocentric.column("DEC"); const delta = topocentric.column("delta"); const delta_t = topocentric.column("TDB-UT"); const lon = topocentric.column("ObsEcLon"); const lat = topocentric.column("ObsEcLat"); const azimuth = topocentric.column("Azimuth"); const elevation = topocentric.column("Elevation"); const g_jd = geocentric.column("Date"); const g_ra = geocentric.column("R.A."); const g_dec = geocentric.column("DEC"); const g_delta = geocentric.column("delta"); const g_lon = geocentric.column("ObsEcLon"); const g_lat = geocentric.column("ObsEcLat");
for (topocentric.rows[0..topocentric.row_count], geocentric.rows[0..geocentric.row_count]) |row, g_row| { std.debug.assert(row[jd] == g_row[g_jd]); const tt = Instant.fromJulianDay(row[jd], .tt); // The Earth turned by Horizons' own UT, so that what is compared is // the parallax and not two models of ΔT. const ut = Instant.fromJulianDay(row[jd] - row[delta_t] / 86400, .ut); const equatorial, const ecliptic, const distance_km = switch (body) { .sun => p: { const p = astro.sun.position(tt); break :p .{ p.equatorial, p.ecliptic, p.distance_au * km_per_au }; }, .moon => p: { const p = astro.moon.position(tt); break :p .{ p.equatorial, p.ecliptic, p.distance_km }; }, }; const t = astro.observer.topocentric(ut, observer, equatorial, distance_km); const cos_dec = @cos(row[dec] * math.rad_per_deg);
const ours_shift = [5]f64{ Angle.fromDegrees(t.equatorial.right_ascension.degrees() - equatorial.right_ascension.degrees()).signed().arcseconds() * cos_dec, (t.equatorial.declination.degrees() - equatorial.declination.degrees()) * 3600, Angle.fromDegrees(t.ecliptic.longitude.degrees() - ecliptic.longitude.degrees()).signed().arcseconds(), (t.ecliptic.latitude.degrees() - ecliptic.latitude.degrees()) * 3600, t.distance_km - distance_km, }; const their_shift = [5]f64{ Angle.fromDegrees(row[ra] - g_row[g_ra]).signed().arcseconds() * cos_dec, (row[dec] - g_row[g_dec]) * 3600, Angle.fromDegrees(row[lon] - g_row[g_lon]).signed().arcseconds(), (row[lat] - g_row[g_lat]) * 3600, (row[delta] - g_row[g_delta]) * km_per_au, }; for (&worst.shift, ours_shift, their_shift) |*w, a, b| w.* = @max(w.*, @abs(a - b)); const columns: Columns = .{ .ra = ra, .dec = dec, .lon = lon, .lat = lat, .delta = delta, .azimuth = azimuth, .elevation = elevation, }; for (&worst.position, errors(row, columns, t)) |*w, e| w.* = @max(w.*, e);
if (row[jd] < last_measured_delta_t) { const wrapped = switch (body) { .sun => astro.sun.topocentric(tt, observer), .moon => astro.moon.topocentric(tt, observer), }; for (&worst.wrapper, errors(row, columns, wrapped)) |*w, e| w.* = @max(w.*, e); } }}
const Columns = struct { ra: usize, dec: usize, lon: usize, lat: usize, delta: usize, azimuth: usize, elevation: usize,};
fn errors(row: [horizons.max_columns]f64, c: Columns, t: astro.observer.Topocentric) [7]f64 { return .{ separation(row[c.ra], t.equatorial.right_ascension.degrees()) * @cos(row[c.dec] * math.rad_per_deg), @abs(row[c.dec] - t.equatorial.declination.degrees()) * 3600, separation(row[c.lon], t.ecliptic.longitude.degrees()), @abs(row[c.lat] - t.ecliptic.latitude.degrees()) * 3600, @abs(row[c.delta] * km_per_au - t.distance_km), separation(row[c.azimuth], t.horizontal.azimuth.degrees()) * @cos(row[c.elevation] * math.rad_per_deg), @abs(row[c.elevation] - t.horizontal.altitude.degrees()) * 3600, };}
/// The worst of each kind of difference over the three sites.fn worstOver(body: Body) Worst { var worst: Worst = .{}; for (&sites) |*site| switch (body) { .sun => measure(&site.sun, &geocentric_sun, site.observer, .sun, &worst), .moon => measure(&site.moon, &geocentric_moon, site.observer, .moon, &worst), }; return worst;}
fn print(name: []const u8, w: Worst) void { std.debug.print("{s}: α δ λ β Δ: shift {any}, position {any}, through the body {any}\n", .{ name, w.shift, w.position, w.wrapper });}
/// Every difference below its tolerance, in the order of `Worst`.fn expectBelow(comptime n: usize, worst: [n]f64, tolerance: [n]f64) !void { for (worst, tolerance) |w, t| try testing.expect(w < t);}
// JPL Horizons (DE441) from Palomar, Tromsø and Cerro Paranal at 201// instants over 1900–2100: see tools/horizons.sh. Each list below is right// ascension (times cos δ), declination, ecliptic longitude and ecliptic// latitude in seconds of arc, then distance in kilometers, then azimuth// (times the cosine of the altitude) and airless altitude in seconds of// arc.//// The shift from geocentric to topocentric, compared with Horizons' own// shift between its geocentric and topocentric tables, tests the parallax// alone. The largest differences seen are 0.29″, 0.11″, 0.29″, 0.12″ and// 0.98 km. They are diurnal aberration, which Horizons applies and this// library leaves out: at most 0.32″ cos φ, eastwards, and the difference in// right ascension is 0.27″ at Palomar and 0.11″ at Tromsø, as cos φ has it.// They are held to 0.3″, 0.12″, 0.3″, 0.12″ and 1 km.//// The topocentric position itself is as good as the geocentric one: the// largest differences are 28.82″, 2.98″, 28.92″, 1.29″ and 10,893 km, the// 0.01° and 0.00007 au of Meeus's low-accuracy Sun, held to 29″, 3″, 29″,// 1.3″ and 10,900 km. The azimuth and altitude carry that, and before 1962// the larger of the two is 28.65″; from 1962 the table's UT is UTC, which// turns the Earth by up to 0.9 s, 13.5″, more or less than UT1 would, and// the largest are 36.94″ and 35.66″, held to 37″ and 36″.//// The Earth is turned by Horizons' UT in both, so that ΔT is not what is// being compared. `sun.topocentric` turns it by this library's ΔT, and is// compared up to the last measured value, 2026 October 1; after that the// prediction and Horizons' part ways, by more than a minute by 2100. Up to// then it does no worse, and, as this library's ΔT and UT1 agree better// than UTC does, its azimuth and altitude do a little better: 30.73″ and// 27.43″, held to 31″ and 28″.test "the topocentric Sun against JPL Horizons, 1900–2100" { const w = worstOver(.sun); errdefer print("sun", w); try expectBelow(5, w.shift, .{ 0.3, 0.12, 0.3, 0.12, 1 }); try expectBelow(7, w.position, .{ 29, 3, 29, 1.3, 10900, 37, 36 }); try expectBelow(7, w.wrapper, .{ 29, 3, 29, 1.3, 10900, 31, 28 });}
// As for the Sun. The shift differs from Horizons' by at most 0.46″,// 0.22″, 0.45″, 0.21″ and 0.90 km: diurnal aberration again, and the// parallax of a Moon whose geocentric distance is up to 46 km out. They are// held to 0.5″, 0.25″, 0.5″, 0.25″ and 1 km. The position: 10.94″, 3.59″,// 10.97″, 3.11″ and 46.5 km, the geocentric theory's own errors, held to// 11″, 3.6″, 11″, 3.2″ and 47 km. The azimuth and altitude: before 1962,// 5.64″ at most; from 1962, with UTC for UT1 as for the Sun, 11.66″ and// 11.88″, held to 12″ and 12″. Through `moon.topocentric` to 2026 October// 1: 11.10″, 3.59″, 11.13″, 3.11″, 44.7 km, 17.82″ and 16.37″, held to// 11.2″, 3.6″, 11.2″, 3.2″, 47 km, 18″ and 16.5″; the azimuth and altitude// carry, besides, the difference between this library's ΔT and the UT1// Horizons used. After 2026, ΔT's prediction would put the Moon tens of// seconds of arc out by 2100, against Horizons' different one.test "the topocentric Moon against JPL Horizons, 1900–2100" { const w = worstOver(.moon); errdefer print("moon", w); try expectBelow(5, w.shift, .{ 0.5, 0.25, 0.5, 0.25, 1 }); try expectBelow(7, w.position, .{ 11, 3.6, 11, 3.2, 47, 12, 12 }); try expectBelow(7, w.wrapper, .{ 11.2, 3.6, 11.2, 3.2, 47, 18, 16.5 });}
test "the shift is within the parallax, and the distance within the observer's" { var jd: f64 = 2451545.0; while (jd < 2451545.0 + 60) : (jd += 0.173) { const instant = Instant.fromJulianDay(jd, .tt); const geocentric = astro.moon.position(instant); for (&sites) |*site| { const t = astro.moon.topocentric(instant, site.observer); const g = site.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(t.distance_km - geocentric.distance_km) <= rho_km * (1 + 1e-12)); const shift = angularSeparation(geocentric.equatorial, t.equatorial); try testing.expect(shift <= math.asin(rho_km / geocentric.distance_km) * (1 + 1e-9)); // And the hour angle is the local apparent sidereal time less the // topocentric right ascension. const last = astro.earth.greenwichApparentSiderealTime(instant).add(site.observer.longitude); const hour_angle = last.sub(t.equatorial.right_ascension).signed(); try testing.expectApproxEqAbs(0, hour_angle.sub(t.hour_angle).signed().arcseconds(), 1e-6); } }}
fn angularSeparation(a: astro.coord.Equatorial, b: astro.coord.Equatorial) f64 { const d1 = a.declination.radians(); const d2 = b.declination.radians(); const c = @sin(d1) * @sin(d2) + @cos(d1) * @cos(d2) * @cos(a.right_ascension.radians() - b.right_ascension.radians()); return math.acos(@min(1.0, @max(-1.0, c)));}
test "longitude wraps, and any finite height is accepted" { const instant = Instant.fromJulianDay(2460000.5, .tt); const west: Observer = .{ .latitude = Angle.fromDegrees(31.5), .longitude = Angle.fromDegrees(-35.5), .height_m = -430 }; const east: Observer = .{ .latitude = Angle.fromDegrees(31.5), .longitude = Angle.fromDegrees(324.5), .height_m = -430 }; const a = astro.moon.topocentric(instant, west); const b = astro.moon.topocentric(instant, east); try testing.expectApproxEqAbs(a.equatorial.right_ascension.degrees(), b.equatorial.right_ascension.degrees(), 1e-9); try testing.expectApproxEqAbs(a.equatorial.declination.degrees(), b.equatorial.declination.degrees(), 1e-9); try testing.expectApproxEqAbs(a.distance_km, b.distance_km, 1e-6);}
test "at the poles and the equator, ρ sin φ′ and ρ cos φ′ are the ellipsoid's" { const north = (Observer{ .latitude = Angle.fromDegrees(90), .longitude = Angle.fromDegrees(0) }).geocentric(); try testing.expectApproxEqAbs(1 - 1.0 / 298.257223563, north.rho_sin_phi, 1e-15); try testing.expectApproxEqAbs(0, north.rho_cos_phi, 1e-15); const equator = (Observer{ .latitude = Angle.fromDegrees(0), .longitude = Angle.fromDegrees(0), .height_m = 6378.137 }).geocentric(); try testing.expectApproxEqAbs(0, equator.rho_sin_phi, 1e-15); try testing.expectApproxEqAbs(1.001, equator.rho_cos_phi, 1e-15);}
test "a latitude beyond ±90°, or anything not finite, gives NaN" { const instant = Instant.fromJulianDay(2460000.5, .tt); const nan = math.nan(f64); const bad = [_]Observer{ .{ .latitude = Angle.fromDegrees(90.000001), .longitude = Angle.fromDegrees(0) }, .{ .latitude = Angle.fromDegrees(-91), .longitude = Angle.fromDegrees(0) }, .{ .latitude = Angle.fromRadians(nan), .longitude = Angle.fromDegrees(0) }, .{ .latitude = Angle.fromDegrees(10), .longitude = Angle.fromRadians(math.inf(f64)) }, .{ .latitude = Angle.fromDegrees(10), .longitude = Angle.fromDegrees(0), .height_m = nan }, }; for (bad) |o| { try testing.expect(math.isNan(o.geocentric().rho_sin_phi)); try expectNan(astro.moon.topocentric(instant, o)); try expectNan(astro.sun.topocentric(instant, o)); } for ([_]f64{ nan, math.inf(f64), -math.inf(f64) }) |jd| { try expectNan(astro.moon.topocentric(Instant.fromJulianDay(jd, .tt), palomar)); try expectNan(astro.sun.topocentric(Instant.fromJulianDay(jd, .ut), palomar)); }}
fn expectNan(t: astro.observer.Topocentric) !void { try testing.expect(math.isNan(t.equatorial.right_ascension.radians())); try testing.expect(math.isNan(t.equatorial.declination.radians())); try testing.expect(math.isNan(t.ecliptic.longitude.radians())); try testing.expect(math.isNan(t.ecliptic.latitude.radians())); try testing.expect(math.isNan(t.hour_angle.radians())); try testing.expect(math.isNan(t.horizontal.azimuth.radians())); try testing.expect(math.isNan(t.horizontal.altitude.radians())); try testing.expect(math.isNan(t.distance_km));}
// The README's example. Horizons gives azimuth 152.961975° and airless// altitude 51.135899° for the Moon from this site at this instant; the// README rounds to 0.1°, and this is within the 18″ the comparison with// Horizons above holds `moon.topocentric` to.test "README: the Moon from Palomar" { const p: Observer = .{ .latitude = Angle.fromDegrees(33.3561), .longitude = Angle.fromDegrees(-116.8625), // east is positive .height_m = 1706, }; const now = Instant.fromCalendar( .gregorian, .init(2026, .october, 23, 4, 0, 0), .ut, ); const moon = astro.moon.topocentric(now, p); try testing.expectApproxEqAbs(152.961975, moon.horizontal.azimuth.degrees(), 18.0 / 3600.0); try testing.expectApproxEqAbs(51.135899, moon.horizontal.altitude.degrees(), 18.0 / 3600.0); try testing.expectApproxEqAbs(153.0, moon.horizontal.azimuth.degrees(), 0.05); try testing.expectApproxEqAbs(51.1, moon.horizontal.altitude.degrees(), 0.05);}