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.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>// SPDX-License-Identifier: MIT
//! `astro.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 horizons = @import("horizons.zig");
// Meeus, Astronomical Algorithms, 2nd ed., example 47.a (pp. 342–343): 1992// April 12 at 0h TD.const example_47a = Instant.fromJulianDay(2448724.5, .tt);
test "example 47.a: from the calendar date" { const instant = Instant.fromCalendar(.gregorian, .{ .year = 1992, .month = .april, .day = 12.0 }, .tt); try testing.expectEqual(example_47a.julian_day, instant.julian_day); try testing.expectApproxEqAbs(-0.077221081451, instant.julianCenturies(), 1e-12);}
// The book prints these to 0.000001° and the distance to 0.1 km. The last// digit can differ by one, because the book carries its nutation and// obliquity rounded to 0.001″ and 0.01″ into the apparent values.test "example 47.a: the apparent position" { const p = astro.moon.position(example_47a); try testing.expectApproxEqAbs(133.167265, p.ecliptic.longitude.degrees(), 2e-6); try testing.expectApproxEqAbs(-3.229126, p.ecliptic.latitude.degrees(), 1e-6); try testing.expectApproxEqAbs(368409.7, p.distance_km, 0.05); try testing.expectApproxEqAbs(134.688470, p.equatorial.right_ascension.degrees(), 2e-6); try testing.expectApproxEqAbs(13.768368, p.equatorial.declination.degrees(), 2e-6);}
// The same example's "exact values, obtained by using the complete// ELP-2000/82 theory", to the 10″ and 4″ Meeus gives for the method.test "example 47.a: within the method's accuracy of the complete theory" { const p = astro.moon.position(example_47a); const lon = Angle.fromDms(.{ .degrees = 133, .minutes = 10, .seconds = 0 }); const lat = Angle.fromDms(.{ .negative = true, .degrees = 3, .minutes = 13, .seconds = 45 }); const ra = Angle.fromHms(.{ .hours = 8, .minutes = 58, .seconds = 45.1 }); const dec = Angle.fromDms(.{ .degrees = 13, .minutes = 46, .seconds = 6 }); try testing.expectApproxEqAbs(lon.arcseconds(), p.ecliptic.longitude.arcseconds(), 10); try testing.expectApproxEqAbs(lat.arcseconds(), p.ecliptic.latitude.arcseconds(), 4); try testing.expectApproxEqAbs(ra.arcseconds(), p.equatorial.right_ascension.arcseconds(), 10); try testing.expectApproxEqAbs(dec.arcseconds(), p.equatorial.declination.arcseconds(), 4); try testing.expectApproxEqAbs(368405.6, p.distance_km, 5);}
const table = horizons.parse(@embedFile("reference/moon.txt"));const km_per_au = 149597870.7;
/// 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());}
// JPL Horizons (DE441) at 201 instants over 1900–2100: see// tools/horizons.sh. Meeus gives the method as "approximately 10″" in// longitude and 4″ in latitude. The largest errors seen are 10.66″ in// longitude, 3.03″ in latitude, 11.27″ in right ascension, 3.58″ in// declination and 46 km in distance, for which he states nothing; the// tolerances are set just above those.test "the Moon against JPL Horizons, 1900–2100" { const jd = table.column("Date"); const ra = table.column("R.A."); const dec = table.column("DEC"); const delta = table.column("delta"); const lon = table.column("ObsEcLon"); const lat = table.column("ObsEcLat");
var worst: [5]f64 = @splat(0); for (table.rows[0..table.row_count]) |row| { const p = astro.moon.position(Instant.fromJulianDay(row[jd], .tt)); const errors = [_]f64{ separation(row[lon], p.ecliptic.longitude.degrees()), @abs(row[lat] - p.ecliptic.latitude.degrees()) * 3600, separation(row[ra], p.equatorial.right_ascension.degrees()), @abs(row[dec] - p.equatorial.declination.degrees()) * 3600, @abs(row[delta] * km_per_au - p.distance_km), }; for (&worst, errors) |*w, e| w.* = @max(w.*, e); } errdefer std.debug.print("worst: lon {d:.2}″ lat {d:.2}″ ra {d:.2}″ dec {d:.2}″ distance {d:.1} km\n", .{ worst[0], worst[1], worst[2], worst[3], worst[4] }); try testing.expect(worst[0] < 11); try testing.expect(worst[1] < 4); try testing.expect(worst[2] < 12); try testing.expect(worst[3] < 4); try testing.expect(worst[4] < 50);}
test "non-finite instants give NaN" { for ([_]f64{ math.nan(f64), math.inf(f64), -math.inf(f64) }) |jd| { const p = astro.moon.position(Instant.fromJulianDay(jd, .tt)); 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)); }}