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.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>// SPDX-License-Identifier: MIT
//! `astro.Angle` and `astro.coord` against Meeus's worked examples.
const std = @import("std");const math = std.math;const testing = std.testing;
const astro = @import("astro");const Angle = astro.Angle;const coord = astro.coord;
test "angles read back in every unit" { const a = Angle.fromDegrees(180); try testing.expectApproxEqAbs(math.pi, a.radians(), 1e-15); try testing.expectApproxEqAbs(12, a.hours(), 1e-12); try testing.expectApproxEqAbs(648000, a.arcseconds(), 1e-9); try testing.expectApproxEqAbs(1.0 / 3600.0, Angle.fromArcseconds(1).degrees(), 1e-18); try testing.expectApproxEqAbs(15, Angle.fromHours(1).degrees(), 1e-12);}
test "degrees–minutes–seconds and hours–minutes–seconds" { try testing.expectApproxEqAbs(-0.5, Angle.fromDms(.{ .negative = true, .degrees = 0, .minutes = 30 }).degrees(), 1e-15); try testing.expectApproxEqAbs(-7.5, Angle.fromHms(.{ .negative = true, .hours = 0, .minutes = 30 }).degrees(), 1e-15);}
test "normalized is [0°, 360°) and signed is (−180°, +180°]" { const cases = [_][3]f64{ // degrees, normalized, signed .{ 0, 0, 0 }, .{ 360, 0, 0 }, .{ -1, 359, -1 }, .{ 180, 180, 180 }, .{ -180, 180, 180 }, .{ 181, 181, -179 }, .{ 725, 5, 5 }, .{ -725, 355, -5 }, }; for (cases) |c| { try testing.expectApproxEqAbs(c[1], Angle.fromDegrees(c[0]).normalized().degrees(), 1e-9); try testing.expectApproxEqAbs(c[2], Angle.fromDegrees(c[0]).signed().degrees(), 1e-9); } // A negative angle too small to survive `@mod` must not come back as 360°. const tiny = Angle.fromRadians(-1e-300).normalized().radians(); try testing.expect(tiny >= 0 and tiny < math.tau); try testing.expect(math.isNan(Angle.fromRadians(math.inf(f64)).normalized().radians())); try testing.expect(math.isNan(Angle.fromRadians(math.nan(f64)).signed().radians()));}
// Meeus, Astronomical Algorithms, 2nd ed., example 13.a (p. 95): Pollux,// whose coordinates the book gives both in sexagesimal and in degrees.const pollux: coord.Equatorial = .{ .right_ascension = Angle.fromHms(.{ .hours = 7, .minutes = 45, .seconds = 18.946 }), .declination = Angle.fromDms(.{ .degrees = 28, .minutes = 1, .seconds = 34.26 }),};const obliquity_j2000 = Angle.fromDegrees(23.4392911);
test "example 13.a: the sexagesimal coordinates of Pollux" { try testing.expectApproxEqAbs(116.328942, pollux.right_ascension.degrees(), 5e-7); try testing.expectApproxEqAbs(28.026183, pollux.declination.degrees(), 5e-7);}
test "example 13.a: Pollux from equatorial to ecliptic coordinates" { const ecliptic = pollux.toEcliptic(obliquity_j2000); // Printed to 1e-6°. try testing.expectApproxEqAbs(113.215630, ecliptic.longitude.degrees(), 1e-6); try testing.expectApproxEqAbs(6.684170, ecliptic.latitude.degrees(), 1e-6);}
// The exercise after example 13.a: back again.test "example 13.a's exercise: Pollux from ecliptic back to equatorial" { const ecliptic: coord.Ecliptic = .{ .longitude = Angle.fromDegrees(113.215630), .latitude = Angle.fromDegrees(6.684170), }; const equatorial = ecliptic.toEquatorial(obliquity_j2000); try testing.expectApproxEqAbs(116.328942, equatorial.right_ascension.degrees(), 1e-6); try testing.expectApproxEqAbs(28.026183, equatorial.declination.degrees(), 1e-6);}
test "conversions keep every quadrant and both poles" { var lon: f64 = 0; while (lon < 360) : (lon += 7.5) { var lat: f64 = -90; while (lat <= 90) : (lat += 7.5) { const ecliptic: coord.Ecliptic = .{ .longitude = Angle.fromDegrees(lon), .latitude = Angle.fromDegrees(lat) }; const equatorial = ecliptic.toEquatorial(obliquity_j2000); try testing.expect(equatorial.right_ascension.degrees() >= 0 and equatorial.right_ascension.degrees() < 360); try testing.expect(@abs(equatorial.declination.degrees()) <= 90); const back = equatorial.toEcliptic(obliquity_j2000); try testing.expectApproxEqAbs(lat, back.latitude.degrees(), 1e-9); // At a pole the longitude is undefined. if (@abs(lat) < 90) { try testing.expectApproxEqAbs(0, back.longitude.sub(Angle.fromDegrees(lon)).signed().degrees(), 1e-9); } } }}
// Example 13.b (pp. 95–96): Venus from the U.S. Naval Observatory in// Washington on 1987 April 10 at 19h21m00s UT, from its apparent right// ascension and declination. Meeus gives the longitude as +77°03′56″,// counted west; here it is east, and negative. He works with it as// 5h08m15.7s, though 77°03′56″ is 5h08m15.733s, and his hour angle carries// that rounding, 0.00014°, so the longitude is his rounded one here. His// azimuth counts from the south, so his 68.0337° is 248.0337° here.test "example 13.b: Venus from Washington" { const instant = astro.time.Instant.fromCalendar(.gregorian, .init(1987, .april, 10, 19, 21, 0), .ut); const longitude = Angle.fromHms(.{ .hours = 5, .minutes = 8, .seconds = 15.7 }).negate(); const latitude = Angle.fromDms(.{ .degrees = 38, .minutes = 55, .seconds = 17 }); const venus: coord.Equatorial = .{ .right_ascension = Angle.fromHms(.{ .hours = 23, .minutes = 9, .seconds = 16.641 }), .declination = Angle.fromDms(.{ .negative = true, .degrees = 6, .minutes = 43, .seconds = 11.61 }), }; // The apparent sidereal time, 8h34m56.853s, printed to 0.001 s. const theta = astro.earth.greenwichApparentSiderealTime(instant); const expected_theta = Angle.fromHms(.{ .hours = 8, .minutes = 34, .seconds = 56.853 }); try testing.expectApproxEqAbs(expected_theta.hours() * 3600, theta.hours() * 3600, 0.0005);
const hour_angle = theta.add(longitude).sub(venus.right_ascension).normalized(); // Printed to 0.000001°; it carries the sidereal time's rounding to // 0.001 s, 0.000004°. try testing.expectApproxEqAbs(64.352133, hour_angle.degrees(), 0.000005);
const horizontal = venus.toHorizontal(hour_angle, latitude); try testing.expectApproxEqAbs(68.0337 + 180.0, horizontal.azimuth.degrees(), 0.00005); try testing.expectApproxEqAbs(15.1249, horizontal.altitude.degrees(), 0.00005);}
test "on the meridian a body is due south or due north, and on the equator it rises due east" { const latitude = Angle.fromDegrees(40); var dec: f64 = -85; while (dec <= 85) : (dec += 5) { const body: coord.Equatorial = .{ .right_ascension = Angle.fromDegrees(0), .declination = Angle.fromDegrees(dec) }; const upper = body.toHorizontal(Angle.fromDegrees(0), latitude); // Upper culmination: south of the zenith when δ < φ, north of it // when δ > φ, at an altitude of 90° − |φ − δ|. // At δ = φ it is at the zenith, where azimuth means nothing. if (dec != 40) try testing.expectApproxEqAbs(@as(f64, if (dec < 40) 180 else 0), azimuthDegrees(upper), 1e-9); try testing.expectApproxEqAbs(90 - @abs(40 - dec), upper.altitude.degrees(), 1e-9); // Lower culmination is due north, at δ + φ − 90°, for a body that // passes between the pole and the nadir; one with δ < −φ passes // beyond the nadir, and is due south. if (dec > -40) { const lower = body.toHorizontal(Angle.fromDegrees(180), latitude); try testing.expectApproxEqAbs(0, azimuthDegrees(lower), 1e-9); try testing.expectApproxEqAbs(dec + 40 - 90, lower.altitude.degrees(), 1e-9); } } // A body on the equator six hours east of the meridian is on the // horizon, due east, from any latitude short of a pole. const equator: coord.Equatorial = .{ .right_ascension = Angle.fromDegrees(0), .declination = Angle.fromDegrees(0) }; var phi: f64 = -85; while (phi <= 85) : (phi += 5) { const rising = equator.toHorizontal(Angle.fromDegrees(-90), Angle.fromDegrees(phi)); try testing.expectApproxEqAbs(90, rising.azimuth.degrees(), 1e-9); try testing.expectApproxEqAbs(0, rising.altitude.degrees(), 1e-9); }}
/// The azimuth in degrees, in [−90°, 270°), so that one rounding has put/// just below 360° reads as just below 0°.fn azimuthDegrees(h: coord.Horizontal) f64 { const a = h.azimuth.signed().degrees(); return if (a < -90) a + 360 else a;}
test "horizontal coordinates of a non-finite direction are NaN" { const nan = Angle.fromRadians(math.nan(f64)); const body: coord.Equatorial = .{ .right_ascension = Angle.fromDegrees(0), .declination = nan }; const h = body.toHorizontal(Angle.fromDegrees(0), Angle.fromDegrees(40)); try testing.expect(math.isNan(h.azimuth.radians()) and math.isNan(h.altitude.radians())); const h2 = (coord.Equatorial{ .right_ascension = Angle.fromDegrees(0), .declination = Angle.fromDegrees(0) }).toHorizontal(nan, Angle.fromDegrees(40)); try testing.expect(math.isNan(h2.azimuth.radians()) and math.isNan(h2.altitude.radians()));}