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.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>// SPDX-License-Identifier: MIT
//! The orientation of the Earth: nutation and the obliquity of the ecliptic//! (Meeus, *Astronomical Algorithms*, 2nd ed., ch. 22), and its rotation,//! as sidereal time (ch. 12).//!//! Every calculation takes an `Instant` in either time scale; the ΔT model//! converts it to whichever one the calculation is written in: TT for//! nutation and obliquity, UT for sidereal time.
const std = @import("std");const math = std.math;
const Angle = @import("Angle.zig");const Instant = @import("time.zig").Instant;
/// The nutation at an instant: the short-period nodding of the Earth's axis.pub const Nutation = struct { /// Nutation in longitude, Δψ: the shift of the true equinox along the /// ecliptic. It adds to every apparent ecliptic longitude. longitude: Angle, /// Nutation in obliquity, Δε: the change in the angle between the /// equator and the ecliptic. It adds to the mean obliquity to give the /// true obliquity. obliquity: Angle,};
/// The nutation in longitude and in obliquity at an instant, from the 1980/// IAU theory as Meeus gives it in Table 22.A: every term with a coefficient/// of at least 0.0003″, which is all 63 of the table. NaN if the instant is/// not finite.pub fn nutation(instant: Instant) Nutation { const t = instant.julianCenturies(); if (!math.isFinite(t)) return .{ .longitude = nan, .obliquity = nan };
// The fundamental arguments, in degrees (Meeus, p. 144). const d = polynomial(t, .{ 297.85036, 445267.111480, -0.0019142, 1.0 / 189474.0 }); const m = polynomial(t, .{ 357.52772, 35999.050340, -0.0001603, -1.0 / 300000.0 }); const m_moon = polynomial(t, .{ 134.96298, 477198.867398, 0.0086972, 1.0 / 56250.0 }); const f = polynomial(t, .{ 93.27191, 483202.017538, -0.0036825, 1.0 / 327270.0 }); const omega = polynomial(t, .{ 125.04452, -1934.136261, 0.0020708, 1.0 / 450000.0 });
var dpsi: f64 = 0; var deps: f64 = 0; for (table_22a) |term| { const argument = (@as(f64, @floatFromInt(term.d)) * d + @as(f64, @floatFromInt(term.m)) * m + @as(f64, @floatFromInt(term.m_moon)) * m_moon + @as(f64, @floatFromInt(term.f)) * f + @as(f64, @floatFromInt(term.omega)) * omega) * math.rad_per_deg; dpsi += (term.sine + term.sine_t * t) * @sin(argument); deps += (term.cosine + term.cosine_t * t) * @cos(argument); } // The table's unit is 0.0001″. return .{ .longitude = Angle.fromArcseconds(dpsi * 0.0001), .obliquity = Angle.fromArcseconds(deps * 0.0001), };}
/// The mean obliquity of the ecliptic, ε₀: the angle between the ecliptic/// and the mean equator, without nutation. Laskar's expression (Meeus, eq./// 22.3), good to 0.01″ between 1000 and 3000 and to a few seconds of arc/// 10,000 years from 2000. Beyond 10,000 years it is a polynomial past its/// range and its value means nothing. NaN if the instant is not finite.pub fn meanObliquity(instant: Instant) Angle { const u = instant.julianCenturies() / 100; const seconds = polynomial(u, .{ 21.448, -4680.93, -1.55, 1999.25, -51.38, -249.67, -39.05, 7.12, 27.87, 5.79, 2.45, }); return Angle.fromDms(.{ .degrees = 23, .minutes = 26, .seconds = seconds });}
/// The true obliquity of the ecliptic, ε = ε₀ + Δε: the angle between the/// ecliptic and the true equator of the instant, the one to use with/// apparent coordinates. NaN if the instant is not finite.pub fn trueObliquity(instant: Instant) Angle { return meanObliquity(instant).add(nutation(instant).obliquity);}
/// The mean sidereal time at Greenwich at an instant: the hour angle of the/// mean vernal equinox from the Greenwich meridian, in [0°, 360°). Meeus's/// eq. 12.4, the IAU's 1982 expression, which is a function of UT; a TT/// instant is converted with the ΔT model, and a second of error in ΔT is/// 15″ of error here.////// UT is taken to be UT1. UTC, which civil time follows, is kept within/// 0.9 s of it, and that difference, up to 13.5″, is not applied. NaN if the/// instant is not finite.pub fn greenwichMeanSiderealTime(instant: Instant) Angle { const jd = instant.to(.ut, .model).julian_day; if (!math.isFinite(jd)) return nan; const days = jd - Instant.j2000_julian_day; const t = days / 36525; // 360.98564736629 d is split into 360 d and 0.98564736629 d. Whole days // of the first are whole turns, so only the fraction of a day is kept: // otherwise, a few centuries from 2000, the product runs to 10⁸ degrees // and loses milliseconds of time to rounding. const fraction = days - @floor(days); const degrees = 280.46061837 + 360 * fraction + 0.98564736629 * days + t * t * (0.000387933 - t / 38710000.0); return Angle.fromDegrees(degrees).normalized();}
/// The apparent sidereal time at Greenwich at an instant: the hour angle of/// the true vernal equinox, in [0°, 360°). It is the mean sidereal time plus/// the equation of the equinoxes, Δψ cos ε, with the true obliquity (Meeus,/// p. 88); the equation of the equinoxes never exceeds about 1.2 s of time./// This is the sidereal time to subtract apparent right ascensions from./// NaN if the instant is not finite.pub fn greenwichApparentSiderealTime(instant: Instant) Angle { const n = nutation(instant); const epsilon = meanObliquity(instant).add(n.obliquity); const equation_of_the_equinoxes = n.longitude.radians() * @cos(epsilon.radians()); return greenwichMeanSiderealTime(instant).add(Angle.fromRadians(equation_of_the_equinoxes)).normalized();}
const nan = Angle.fromRadians(math.nan(f64));
/// c[0] + c[1]x + c[2]x² + …, by Horner's method.fn polynomial(x: f64, comptime c: anytype) f64 { var sum: f64 = 0; comptime var i = c.len; inline while (i > 0) { i -= 1; sum = sum * x + c[i]; } return sum;}
/// A row of Meeus's Table 22.A: the multiples of D, M, M′, F and Ω that make/// the argument, and the coefficients of its sine in Δψ and of its cosine in/// Δε, in units of 0.0001″, each with a part proportional to T.const Term = struct { d: i8, m: i8, m_moon: i8, f: i8, omega: i8, sine: f64, sine_t: f64 = 0, cosine: f64 = 0, cosine_t: f64 = 0,};
fn row(d: i8, m: i8, m_moon: i8, f: i8, omega: i8, sine: f64, sine_t: f64, cosine: f64, cosine_t: f64) Term { return .{ .d = d, .m = m, .m_moon = m_moon, .f = f, .omega = omega, .sine = sine, .sine_t = sine_t, .cosine = cosine, .cosine_t = cosine_t };}
/// Meeus, Table 22.A (pp. 145–146), in the book's order.const table_22a = [_]Term{ row(0, 0, 0, 0, 1, -171996, -174.2, 92025, 8.9), row(-2, 0, 0, 2, 2, -13187, -1.6, 5736, -3.1), row(0, 0, 0, 2, 2, -2274, -0.2, 977, -0.5), row(0, 0, 0, 0, 2, 2062, 0.2, -895, 0.5), row(0, 1, 0, 0, 0, 1426, -3.4, 54, -0.1), row(0, 0, 1, 0, 0, 712, 0.1, -7, 0), row(-2, 1, 0, 2, 2, -517, 1.2, 224, -0.6), row(0, 0, 0, 2, 1, -386, -0.4, 200, 0), row(0, 0, 1, 2, 2, -301, 0, 129, -0.1), row(-2, -1, 0, 2, 2, 217, -0.5, -95, 0.3), row(-2, 0, 1, 0, 0, -158, 0, 0, 0), row(-2, 0, 0, 2, 1, 129, 0.1, -70, 0), row(0, 0, -1, 2, 2, 123, 0, -53, 0), row(2, 0, 0, 0, 0, 63, 0, 0, 0), row(0, 0, 1, 0, 1, 63, 0.1, -33, 0), row(2, 0, -1, 2, 2, -59, 0, 26, 0), row(0, 0, -1, 0, 1, -58, -0.1, 32, 0), row(0, 0, 1, 2, 1, -51, 0, 27, 0), row(-2, 0, 2, 0, 0, 48, 0, 0, 0), row(0, 0, -2, 2, 1, 46, 0, -24, 0), row(2, 0, 0, 2, 2, -38, 0, 16, 0), row(0, 0, 2, 2, 2, -31, 0, 13, 0), row(0, 0, 2, 0, 0, 29, 0, 0, 0), row(-2, 0, 1, 2, 2, 29, 0, -12, 0), row(0, 0, 0, 2, 0, 26, 0, 0, 0), row(-2, 0, 0, 2, 0, -22, 0, 0, 0), row(0, 0, -1, 2, 1, 21, 0, -10, 0), row(0, 2, 0, 0, 0, 17, -0.1, 0, 0), row(2, 0, -1, 0, 1, 16, 0, -8, 0), row(-2, 2, 0, 2, 2, -16, 0.1, 7, 0), row(0, 1, 0, 0, 1, -15, 0, 9, 0), row(-2, 0, 1, 0, 1, -13, 0, 7, 0), row(0, -1, 0, 0, 1, -12, 0, 6, 0), row(0, 0, 2, -2, 0, 11, 0, 0, 0), row(2, 0, -1, 2, 1, -10, 0, 5, 0), row(2, 0, 1, 2, 2, -8, 0, 3, 0), row(0, 1, 0, 2, 2, 7, 0, -3, 0), row(-2, 1, 1, 0, 0, -7, 0, 0, 0), row(0, -1, 0, 2, 2, -7, 0, 3, 0), row(2, 0, 0, 2, 1, -7, 0, 3, 0), row(2, 0, 1, 0, 0, 6, 0, 0, 0), row(-2, 0, 2, 2, 2, 6, 0, -3, 0), row(-2, 0, 1, 2, 1, 6, 0, -3, 0), row(2, 0, -2, 0, 1, -6, 0, 3, 0), row(2, 0, 0, 0, 1, -6, 0, 3, 0), row(0, -1, 1, 0, 0, 5, 0, 0, 0), row(-2, -1, 0, 2, 1, -5, 0, 3, 0), row(-2, 0, 0, 0, 1, -5, 0, 3, 0), row(0, 0, 2, 2, 1, -5, 0, 3, 0), row(-2, 0, 2, 0, 1, 4, 0, 0, 0), row(-2, 1, 0, 2, 1, 4, 0, 0, 0), row(0, 0, 1, -2, 0, 4, 0, 0, 0), row(-1, 0, 1, 0, 0, -4, 0, 0, 0), row(-2, 1, 0, 0, 0, -4, 0, 0, 0), row(1, 0, 0, 0, 0, -4, 0, 0, 0), row(0, 0, 1, 2, 0, 3, 0, 0, 0), row(0, 0, -2, 2, 2, -3, 0, 0, 0), row(-1, -1, 1, 0, 0, -3, 0, 0, 0), row(0, 1, 1, 0, 0, -3, 0, 0, 0), row(0, -1, 1, 2, 2, -3, 0, 0, 0), row(2, -1, -1, 2, 2, -3, 0, 0, 0), row(0, 0, 3, 2, 2, -3, 0, 0, 0), row(2, -1, 0, 2, 2, -3, 0, 0, 0),};
comptime { std.debug.assert(table_22a.len == 63);}