// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie // SPDX-License-Identifier: MIT //! `astro.moon.libration` against Meeus's worked examples and JPL Horizons, //! from the Earth's center and from three observatories. 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 libration = astro.moon.libration; const horizons = @import("horizons.zig"); // Meeus, Astronomical Algorithms, 2nd ed., example 53.a (p. 374): 1992 April // 12 at 0h TD, the instant of example 47.a. The book prints the librations // to 0.001° (optical, physical) and 0.01° (total), and P to 0.01°. const example_53a = Instant.fromJulianDay(2448724.5, .tt); test "example 53.a: the optical libration" { const l = libration.at(example_53a); try testing.expectApproxEqAbs(-1.206, l.optical.longitude.degrees(), 0.0005); try testing.expectApproxEqAbs(4.194, l.optical.latitude.degrees(), 0.0005); } test "example 53.a: the physical libration" { const l = libration.at(example_53a); try testing.expectApproxEqAbs(-0.025, l.physical.longitude.degrees(), 0.0005); try testing.expectApproxEqAbs(0.006, l.physical.latitude.degrees(), 0.0005); } test "example 53.a: the total libration and the position angle of the axis" { const l = libration.at(example_53a); try testing.expectApproxEqAbs(-1.23, l.total.longitude.degrees(), 0.005); try testing.expectApproxEqAbs(4.20, l.total.latitude.degrees(), 0.005); try testing.expectApproxEqAbs(15.08, l.position_angle_of_axis.degrees(), 0.005); } test "the total is the sum of the optical and the physical" { var jd: f64 = 2451545.0; while (jd < 2451545.0 + 400) : (jd += 3.3) { const l = libration.at(Instant.fromJulianDay(jd, .tt)); const lon = l.optical.longitude.add(l.physical.longitude).sub(l.total.longitude).signed(); try testing.expectApproxEqAbs(0, lon.degrees(), 1e-12); try testing.expectApproxEqAbs(l.optical.latitude.degrees() + l.physical.latitude.degrees(), l.total.latitude.degrees(), 1e-12); } } const table = horizons.parse(@embedFile("reference/moon-libration.txt")); /// The difference between two directions, in degrees, however they wrap. fn separation(a: f64, b: f64) f64 { return @abs(Angle.fromDegrees(a - b).signed().degrees()); } // JPL Horizons (DE441) at 201 instants over 1900–2100: see // tools/horizons.sh. Horizons gives the sub-observer point — the libration — // in the Moon's mean-Earth frame from DE441's integrated rotation, with // longitude positive to the east, the same convention as Meeus's, which // these signs check. The largest differences seen are 0.0081° in longitude // and 0.0081° in the position angle of the axis, held to 0.01°. In latitude // Horizons is 0.0224° below this throughout, scattered by only 0.0025° about // that, with a largest difference of 0.0273°: a steady offset, as a // difference between DE441's lunar frame and the IAU mean lunar equator // Meeus uses would give, and not the scatter of an error in the series. It // is held to 0.03°. test "the libration against JPL Horizons, 1900–2100" { const jd = table.column("Date"); const lon = table.column("ObsSub-LON"); const lat = table.column("ObsSub-LAT"); const pole = table.column("NP.ang"); var worst: [3]f64 = @splat(0); for (table.rows[0..table.row_count]) |row| { const l = libration.at(Instant.fromJulianDay(row[jd], .tt)); const errors = [_]f64{ separation(row[lon], l.total.longitude.degrees()), @abs(row[lat] - l.total.latitude.degrees()), separation(row[pole], l.position_angle_of_axis.degrees()), }; for (&worst, errors) |*w, e| w.* = @max(w.*, e); } errdefer std.debug.print("worst: l {d:.4}° b {d:.4}° P {d:.4}°\n", .{ worst[0], worst[1], worst[2] }); try testing.expect(worst[0] < 0.01); try testing.expect(worst[1] < 0.03); try testing.expect(worst[2] < 0.01); } test "non-finite instants give NaN" { for ([_]f64{ math.nan(f64), math.inf(f64), -math.inf(f64) }) |jd| { const l = libration.at(Instant.fromJulianDay(jd, .tt)); for ([_]libration.Selenographic{ l.optical, l.physical, l.total }) |s| { try testing.expect(math.isNan(s.longitude.radians())); try testing.expect(math.isNan(s.latitude.radians())); } try testing.expect(math.isNan(l.position_angle_of_axis.radians())); } } // Example 53.b (p. 377): the selenographic position of the Sun for the // instant of example 53.a, printed to 0.01°. The example takes the Sun and // the Moon from VSOP87 and the complete ELP-2000/82, so these carry the // 0.01° of this library's Sun. test "example 53.b: the subsolar point and the colongitude" { const p = libration.subsolarPoint(example_53a); try testing.expectApproxEqAbs(67.89, p.longitude.degrees(), 0.01); try testing.expectApproxEqAbs(1.46, p.latitude.degrees(), 0.01); try testing.expectApproxEqAbs(22.11, libration.colongitude(example_53a).degrees(), 0.01); } // Example 53.c (p. 377): sunrise on Copernicus, at selenographic longitude // −20.0°, on 1992 April 11.8065 TD, when the book finds c₀ = 19.75° and // b₀ = +1.46°: the Sun rises on a crater at longitude η when the colongitude // is −η, near enough. test "example 53.c: sunrise on Copernicus" { const instant = Instant.fromCalendar(.gregorian, .{ .year = 1992, .month = .april, .day = 11.8065 }, .tt); try testing.expectApproxEqAbs(19.75, libration.colongitude(instant).degrees(), 0.01); try testing.expectApproxEqAbs(1.46, libration.subsolarPoint(instant).latitude.degrees(), 0.01); } test "the colongitude is 90° minus the subsolar longitude, and grows about 12.19° a day" { var jd: f64 = 2451545.0; while (jd < 2451545.0 + 400) : (jd += 2.9) { const instant = Instant.fromJulianDay(jd, .tt); const c = libration.colongitude(instant); try testing.expect(c.degrees() >= 0 and c.degrees() < 360); const expected = Angle.fromDegrees(90).sub(libration.subsolarPoint(instant).longitude); try testing.expectApproxEqAbs(0, c.sub(expected).signed().degrees(), 1e-9); const tomorrow = libration.colongitude(Instant.fromJulianDay(jd + 1, .tt)); // The mean daily motion is 12.190749°; the libration and the // Moon's uneven motion move it by up to a degree either way. try testing.expectApproxEqAbs(12.19, tomorrow.sub(c).signed().degrees(), 1.2); } } // Horizons' subsolar point, from the same table as the libration. In // longitude the largest difference seen is 0.0079°, held to 0.01°. In // latitude the largest is 0.0260°, held to 0.03°; but it is not scatter. It // is −0.0224° × cos l₀, to a residual of 0.0013° rms: the same 0.0224° by // which Horizons' sub-Earth latitude, always near longitude 0°, sits below // this library's. The two lunar frames differ by a tilt of 0.0224° about // the axis through selenographic longitudes ±90°. With that taken out, the // largest difference in the subsolar latitude is 0.0059°, held to 0.006°. test "the subsolar point against JPL Horizons, 1900–2100" { const jd = table.column("Date"); const lon = table.column("SunSub-LON"); const lat = table.column("SunSub-LAT"); var worst: [3]f64 = @splat(0); for (table.rows[0..table.row_count]) |row| { const p = libration.subsolarPoint(Instant.fromJulianDay(row[jd], .tt)); const lat_difference = row[lat] - p.latitude.degrees(); const tilt = -0.0224 * @cos(p.longitude.radians()); const errors = [_]f64{ separation(row[lon], p.longitude.degrees()), @abs(lat_difference), @abs(lat_difference - tilt), }; for (&worst, errors) |*w, e| w.* = @max(w.*, e); } errdefer std.debug.print("worst: l₀ {d:.4}° b₀ {d:.4}°, less the tilt {d:.4}°\n", .{ worst[0], worst[1], worst[2] }); try testing.expect(worst[0] < 0.01); try testing.expect(worst[1] < 0.03); try testing.expect(worst[2] < 0.006); } const Site = struct { name: []const u8, observer: astro.observer.Observer, table: 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 }, .table = horizons.parse(@embedFile("reference/moon-libration-palomar.txt")), }, .{ .name = "Tromsø", .observer = .{ .latitude = Angle.fromDegrees(69.6492), .longitude = Angle.fromDegrees(18.9553), .height_m = 100 }, .table = horizons.parse(@embedFile("reference/moon-libration-tromso.txt")), }, .{ .name = "Paranal", .observer = .{ .latitude = Angle.fromDegrees(-24.6272), .longitude = Angle.fromDegrees(289.5973), .height_m = 2635 }, .table = horizons.parse(@embedFile("reference/moon-libration-paranal.txt")), }, }; /// The last instant this library's ΔT is measured rather than predicted, /// 2026 October 1. const last_measured_delta_t = 2461314.5; // The Moon's topocentric libration from Palomar, Tromsø and Cerro Paranal, // against Horizons' sub-observer point and north pole angle from each, at // the same 201 instants (see tools/horizons.sh). Up to the last measured // ΔT, the largest differences seen are 0.0082° in longitude and 0.0081° in // the position angle of the axis, held to 0.01°, as the geocentric ones // are. In latitude Horizons is again 0.0224° below this on average, the // tilt between DE441's lunar frame and the IAU mean lunar equator found in // the geocentric libration, and the largest difference is 0.0272°, held to // 0.03°; less −0.0224° × cos l for the tilt, it is 0.0049°, held to 0.005°. // // After 2026 the Earth is turned by this library's prediction of ΔT, and // Horizons by its own; the two have parted by more than a minute by 2100, // which moves the topocentric Moon and so its libration, by up to 0.0133° // at worst among the three. Those instants are held to 0.014°, in // longitude, in latitude less the tilt, and in the position angle. test "the topocentric libration against JPL Horizons, 1900–2100" { var worst: [5]f64 = @splat(0); for (&sites) |*site| { const jd = site.table.column("Date"); const lon = site.table.column("ObsSub-LON"); const lat = site.table.column("ObsSub-LAT"); const pole = site.table.column("NP.ang"); for (site.table.rows[0..site.table.row_count]) |row| { const l = libration.topocentric(Instant.fromJulianDay(row[jd], .tt), site.observer); const lat_difference = row[lat] - l.total.latitude.degrees(); const tilt = -0.0224 * @cos(l.total.longitude.radians()); const lon_error = separation(row[lon], l.total.longitude.degrees()); const pole_error = separation(row[pole], l.position_angle_of_axis.degrees()); if (row[jd] < last_measured_delta_t) { const errors = [_]f64{ lon_error, @abs(lat_difference), @abs(lat_difference - tilt), pole_error }; for (worst[0..4], errors) |*w, e| w.* = @max(w.*, e); } else { worst[4] = @max(worst[4], @max(lon_error, @max(@abs(lat_difference - tilt), pole_error))); } } } errdefer std.debug.print("worst: l {d:.4}° b {d:.4}°, less the tilt {d:.4}°, P {d:.4}°; after 2026 {d:.4}°\n", .{ worst[0], worst[1], worst[2], worst[3], worst[4] }); try testing.expect(worst[0] < 0.01); try testing.expect(worst[1] < 0.03); try testing.expect(worst[2] < 0.005); try testing.expect(worst[3] < 0.01); try testing.expect(worst[4] < 0.014); } // From the Earth's surface the Moon is seen from up to about 1.03° away // from its center's view, the parallax at perigee, and the libration // moves by as much: over 1900–2100 the three sites see it up to 0.989° // from the geocentric. Its two parts are still found the same way, and // still add up. test "the topocentric libration is within the parallax of the geocentric" { var jd: f64 = 2451545.0; while (jd < 2451545.0 + 400) : (jd += 0.37) { const instant = Instant.fromJulianDay(jd, .tt); const g = libration.at(instant); for (&sites) |*site| { const t = libration.topocentric(instant, site.observer); const dl = t.total.longitude.sub(g.total.longitude).signed().radians(); const moved = math.acos(@min(1.0, @sin(t.total.latitude.radians()) * @sin(g.total.latitude.radians()) + @cos(t.total.latitude.radians()) * @cos(g.total.latitude.radians()) * @cos(dl))); try testing.expect(moved * math.deg_per_rad < 1.03); const sum = t.optical.longitude.add(t.physical.longitude).sub(t.total.longitude).signed(); try testing.expectApproxEqAbs(0, sum.degrees(), 1e-12); try testing.expectApproxEqAbs(t.optical.latitude.degrees() + t.physical.latitude.degrees(), t.total.latitude.degrees(), 1e-12); } } } test "the topocentric libration of a non-finite instant, or from a latitude past the pole, is NaN" { const past_the_pole: astro.observer.Observer = .{ .latitude = Angle.fromDegrees(91), .longitude = Angle.fromDegrees(0) }; const cases = [_]struct { Instant, astro.observer.Observer }{ .{ Instant.fromJulianDay(math.nan(f64), .tt), sites[0].observer }, .{ Instant.fromJulianDay(math.inf(f64), .ut), sites[1].observer }, .{ Instant.fromJulianDay(2451545.0, .tt), past_the_pole }, }; for (cases) |c| { const l = libration.topocentric(c[0], c[1]); for ([_]libration.Selenographic{ l.optical, l.physical, l.total }) |p| { try testing.expect(math.isNan(p.longitude.radians())); try testing.expect(math.isNan(p.latitude.radians())); } try testing.expect(math.isNan(l.position_angle_of_axis.radians())); } }