// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie // SPDX-License-Identifier: MIT //! Instants, the time scales they are counted in, and ΔT. //! //! Every calculation in this library takes an `Instant`: a Julian Day together //! with the time scale it is counted in. The theories are written in TT //! (Terrestrial Time), a uniform scale; civil time follows UT (Universal //! Time), which follows the Earth's rotation. The two differ by ΔT, about 69 s //! in 2025 and hours in antiquity, so an instant always says which one it is. //! //! UTC, the scale of computer clocks, is counted as UT. The two never differ //! by more than 0.9 s, which moves the Moon about 0.5″ — far inside the //! accuracy of anything here. //! //! Julian Days and calendar dates follow Meeus, *Astronomical Algorithms*, //! 2nd ed., chapter 7. ΔT is described at `Instant.deltaTSeconds`: the //! Espenak–Meeus polynomials for the past, as given in Espenak, "Polynomial //! Expressions for Delta T", on NASA's eclipse web site //! (), and the U.S. //! Naval Observatory's measured values since 2005 //! (). const std = @import("std"); const math = std.math; /// The clock an `Instant` is counted in. pub const TimeScale = enum { /// Terrestrial Time: the uniform scale the theories are written in. tt, /// Universal Time, which follows the Earth's rotation. UTC is counted as /// UT. ut, }; /// The calendar a `CalendarDate` is written in. /// /// There is no automatic switch between the two: Meeus switches at 1582 /// October 15, but countries adopted the Gregorian calendar at different /// times (Britain in 1752), so the caller says which calendar a date is in. /// Both are proleptic, extending indefinitely backwards and forwards. pub const Calendar = enum { gregorian, julian, }; /// A month of the year, numbered from 1 as calendars number them. pub const Month = enum(u4) { january = 1, february, march, april, may, june, july, august, september, october, november, december, }; /// A date in a calendar, with the time of day carried as the fraction of the /// day, as Meeus writes them: 1957 October 4.81 is `day = 4.81`. /// /// Years are astronomical: the year before 1 is 0, and the year before that /// is −1 (1 BCE is 0, 2 BCE is −1). pub const CalendarDate = struct { year: i32, month: Month, /// The day of the month, starting at 1.0, with the time of day as its /// fractional part. A value past the end of the month runs on into the /// next. day: f64, /// A date from its day of the month and its time of day. pub fn init(year: i32, month: Month, day: u5, hour: u5, minute: u6, second: f64) CalendarDate { const seconds = @as(f64, @floatFromInt(hour)) * 3600.0 + @as(f64, @floatFromInt(minute)) * 60.0 + second; return .{ .year = year, .month = month, .day = @as(f64, @floatFromInt(day)) + seconds / seconds_per_day, }; } }; /// Where a conversion between TT and UT takes ΔT from. pub const DeltaT = union(enum) { /// The Espenak–Meeus polynomials; see `Instant.deltaTSeconds`. model, /// A value supplied by the caller, in seconds: a measured value, or a /// more recent prediction than the model's. seconds: f64, }; /// A single moment: a Julian Day and the time scale it is counted in. pub const Instant = struct { /// The Julian Day: days since noon on −4712 January 1 in the proleptic /// Julian calendar, counted in `scale`. julian_day: f64, scale: TimeScale, /// The Julian Day of the Unix epoch, 1970 January 1 at 0h. pub const unix_epoch_julian_day = 2440587.5; /// An instant from a Julian Day counted in `scale`. Any `f64` is /// accepted; calculations on a non-finite one give NaN. pub fn fromJulianDay(julian_day: f64, scale: TimeScale) Instant { return .{ .julian_day = julian_day, .scale = scale }; } /// An instant from seconds since the Unix epoch, counted as UT. pub fn fromUnixSeconds(seconds: f64) Instant { return .{ .julian_day = seconds / seconds_per_day + unix_epoch_julian_day, .scale = .ut, }; } /// An instant from a date in the given calendar (Meeus, eq. 7.1). pub fn fromCalendar(calendar: Calendar, date: CalendarDate, scale: TimeScale) Instant { var y: f64 = @floatFromInt(date.year); var m: f64 = @floatFromInt(@backingInt(date.month)); if (m <= 2) { y -= 1; m += 12; } const b: f64 = switch (calendar) { .julian => 0, .gregorian => b: { const a = @floor(y / 100); break :b 2 - a + @floor(a / 4); }, }; // `@floor` rather than Meeus's INT, which truncates: the two agree for // the positive Julian Days he restricts himself to, and only `@floor` // keeps counting correctly before −4712. const jd = @floor(365.25 * (y + 4716)) + @floor(30.6001 * (m + 1)) + date.day + b - 1524.5; return .{ .julian_day = jd, .scale = scale }; } /// This instant as a date in the given calendar (Meeus, ch. 7), in the /// instant's own time scale. Null if the Julian Day is not finite or the /// year does not fit in an `i32`. pub fn toCalendar(self: Instant, calendar: Calendar) ?CalendarDate { const jd = self.julian_day; if (!math.isFinite(jd)) return null; // About ±2 billion years, past which the year would not fit. if (@abs(jd) > 7.8e11) return null; const z = @floor(jd + 0.5); const f = jd + 0.5 - z; const a = switch (calendar) { .julian => z, .gregorian => a: { const alpha = @floor((z - 1867216.25) / 36524.25); break :a z + 1 + alpha - @floor(alpha / 4); }, }; const b = a + 1524; const c = @floor((b - 122.1) / 365.25); const d = @floor(365.25 * c); const e = @floor((b - d) / 30.6001); const day = b - d - @floor(30.6001 * e) + f; const month = if (e < 14) e - 1 else e - 13; const year = if (month > 2) c - 4716 else c - 4715; if (!(month >= 1 and month <= 12)) return null; if (!(year >= math.minInt(i32) and year <= math.maxInt(i32))) return null; return .{ .year = @intFromFloat(year), .month = @fromBackingInt(@as(u4, @intFromFloat(month))), .day = day, }; } /// Seconds since the Unix epoch, in UT. A TT instant is converted first, /// taking ΔT from `delta_t`. pub fn toUnixSeconds(self: Instant, delta_t: DeltaT) f64 { const ut = self.to(.ut, delta_t); return (ut.julian_day - unix_epoch_julian_day) * seconds_per_day; } /// This instant counted in `scale`, taking ΔT from `delta_t` if the scale /// changes. NaN if the Julian Day is not finite. pub fn to(self: Instant, scale: TimeScale, delta_t: DeltaT) Instant { if (!math.isFinite(self.julian_day)) return .{ .julian_day = math.nan(f64), .scale = scale }; if (self.scale == scale) return self; const jd = switch (scale) { // ΔT is a function of UT here, which is what we have. .tt => self.julian_day + deltaTAt(self.julian_day, delta_t) / seconds_per_day, // ΔT is a function of UT here, which is what we want, so it is // found by iteration. ΔT changes by at most a few seconds a year, // so each pass shrinks the error by about seven orders of // magnitude and three are ample. .ut => ut: { var ut = self.julian_day; for (0..3) |_| ut = self.julian_day - deltaTAt(ut, delta_t) / seconds_per_day; break :ut ut; }, }; return .{ .julian_day = jd, .scale = scale }; } /// Julian centuries of TT since J2000.0 (JDE 2451545.0), the T in which /// the theories are written (Meeus, eq. 22.1). A UT instant is converted /// with the ΔT model. NaN if the Julian Day is not finite. pub fn julianCenturies(self: Instant) f64 { return (self.to(.tt, .model).julian_day - j2000_julian_day) / 36525; } /// The Julian Day of J2000.0, 2000 January 1.5 TT. pub const j2000_julian_day = 2451545.0; /// ΔT = TT − UT at this instant, in seconds, from the model this library /// uses unless a caller supplies a value as `DeltaT.seconds`. It is built /// from four pieces: /// /// - **Before 2005**, the Espenak–Meeus polynomials. From −500 to 1950 /// they follow the historical record of Morrison and Stephenson (2004), /// and from 1955 to 2005 they reproduce measured values to within about /// 0.1 s. Before −500 they are the long-term parabola −20 + 32u², /// u = (y − 1820)/100. /// - **From 2005 to the last measured value** (2026 October 1), the U.S. /// Naval Observatory's measurements, interpolated linearly between the /// start of each year. /// - **From then to 2150**, a prediction: a cubic that leaves the last /// measured value at its measured rate of change and arrives at the /// long-term parabola in 2150 with the parabola's value and slope. It /// gives about 88 s for 2050 and 199 s for 2100. No prediction of ΔT /// is good to better than tens of seconds a few decades out. /// - **After 2150**, the long-term parabola. /// /// Outside the measured span, 1955 to the last measured value, the /// values are those for use with the ELP-2000/82 lunar theory this /// library computes the Moon from: Morrison and Stephenson's values /// assume a tidal acceleration of the Moon of −26″/cy², ELP-2000/82 uses /// −25.858″/cy², and the difference is the correction /// c = −0.000012932 (y − 1955)² s that Espenak and Meeus give for it. /// Meeus applies the same kind of correction to his own table of ΔT /// before 1955 (*Astronomical Algorithms*, 2nd ed., p. 78). It is under /// 0.05 s across 1900–2100 and about 200 s at −2000. Measured values come /// from atomic clocks and need no correction. /// /// The instant's own Julian Day is used whichever scale it is in: the /// difference is under a day, and ΔT does not change measurably in a day. pub fn deltaTSeconds(self: Instant) f64 { return modelDeltaT(self.julian_day); } }; const seconds_per_day = 86400.0; fn deltaTAt(julian_day: f64, delta_t: DeltaT) f64 { return switch (delta_t) { .model => modelDeltaT(julian_day), .seconds => |s| s, }; } /// The decimal year Espenak and Meeus evaluate their polynomials at. They /// define it as y = year + (month − 0.5)/12, the middle of a month; counting /// continuously in mean Gregorian years from 2000 January 1.0 agrees with /// that at the middle of every month and is smooth in between. fn decimalYear(julian_day: f64) f64 { return 2000.0 + (julian_day - 2451544.5) / 365.2425; } /// ΔT measured by the U.S. Naval Observatory: the Julian Day (UT) and ΔT in /// seconds, at the start of each year from 2005, and the last value /// published. Extending it is a matter of adding rows from /// ; the prediction picks up /// from wherever the table ends. const measured = [_][2]f64{ .{ 2453371.5, 64.6876 }, // 2005 January 1 .{ 2453736.5, 64.8452 }, // 2006 January 1 .{ 2454101.5, 65.1464 }, // 2007 January 1 .{ 2454466.5, 65.4573 }, // 2008 January 1 .{ 2454832.5, 65.7768 }, // 2009 January 1 .{ 2455197.5, 66.0699 }, // 2010 January 1 .{ 2455562.5, 66.3246 }, // 2011 January 1 .{ 2455927.5, 66.6030 }, // 2012 January 1 .{ 2456293.5, 66.9069 }, // 2013 January 1 .{ 2456658.5, 67.2810 }, // 2014 January 1 .{ 2457023.5, 67.6439 }, // 2015 January 1 .{ 2457388.5, 68.1024 }, // 2016 January 1 .{ 2457754.5, 68.5927 }, // 2017 January 1 .{ 2458119.5, 68.9676 }, // 2018 January 1 .{ 2458484.5, 69.2202 }, // 2019 January 1 .{ 2458849.5, 69.3612 }, // 2020 January 1 .{ 2459215.5, 69.3594 }, // 2021 January 1 .{ 2459580.5, 69.2945 }, // 2022 January 1 .{ 2459945.5, 69.2039 }, // 2023 January 1 .{ 2460310.5, 69.1752 }, // 2024 January 1 .{ 2460676.5, 69.1377 }, // 2025 January 1 .{ 2461041.5, 69.1099 }, // 2026 January 1 .{ 2461314.5, 69.2065 }, // 2026 October 1 }; /// Where the prediction ends and the long-term parabola takes over. const prediction_end_year = 2150.0; fn modelDeltaT(julian_day: f64) f64 { if (!math.isFinite(julian_day)) return math.nan(f64); if (julian_day >= measured[0][0]) { const last = measured[measured.len - 1]; if (julian_day <= last[0]) return interpolateMeasured(julian_day); return predicted(decimalYear(julian_day)); } const y = decimalYear(julian_day); // From 1955 the polynomials follow atomic-clock measurements, which need // no correction. return pastPolynomials(y) + if (y < 1955) lunarCorrection(y) else 0; } /// ΔT between two rows of `measured`, linearly. fn interpolateMeasured(julian_day: f64) f64 { var i: usize = 1; while (i < measured.len - 1 and measured[i][0] < julian_day) i += 1; const a = measured[i - 1]; const b = measured[i]; return a[1] + (b[1] - a[1]) * (julian_day - a[0]) / (b[0] - a[0]); } /// Beyond the measurements: a cubic Hermite curve from the last measured /// value, leaving it at the rate between the last two rows, to the corrected /// long-term parabola at `prediction_end_year`, meeting it with its slope. fn predicted(y: f64) f64 { if (y >= prediction_end_year) return parabola(y) + lunarCorrection(y); const a = measured[measured.len - 2]; const b = measured[measured.len - 1]; const y0 = decimalYear(b[0]); const v0 = b[1]; const m0 = (b[1] - a[1]) / (y0 - decimalYear(a[0])); const y1 = prediction_end_year; const v1 = parabola(y1) + lunarCorrection(y1); const m1 = parabolaSlope(y1) + lunarCorrectionSlope(y1); const h = y1 - y0; const t = (y - y0) / h; const t2 = t * t; const t3 = t2 * t; return (2 * t3 - 3 * t2 + 1) * v0 + (t3 - 2 * t2 + t) * h * m0 + (-2 * t3 + 3 * t2) * v1 + (t3 - t2) * h * m1; } /// Espenak and Meeus's correction for using ΔT with ELP-2000/82, which takes /// a slightly different tidal acceleration of the Moon (−25.858″/cy²) from /// the −26″/cy² Morrison and Stephenson's values assume. Zero in 1955; the /// measured values after that need none. fn lunarCorrection(y: f64) f64 { const d = y - 1955; return -0.000012932 * d * d; } fn lunarCorrectionSlope(y: f64) f64 { return -2 * 0.000012932 * (y - 1955); } /// The Espenak–Meeus polynomials up to 2005, without the lunar correction. fn pastPolynomials(y: f64) f64 { if (y < -500) return parabola(y); if (y < 500) { const u = y / 100; return poly(u, &.{ 10583.6, -1014.41, 33.78311, -5.952053, -0.1798452, 0.022174192, 0.0090316521 }); } if (y < 1600) { const u = (y - 1000) / 100; return poly(u, &.{ 1574.2, -556.01, 71.23472, 0.319781, -0.8503463, -0.005050998, 0.0083572073 }); } if (y < 1700) { const t = y - 1600; return poly(t, &.{ 120, -0.9808, -0.01532, 1.0 / 7129.0 }); } if (y < 1800) { const t = y - 1700; return poly(t, &.{ 8.83, 0.1603, -0.0059285, 0.00013336, -1.0 / 1174000.0 }); } if (y < 1860) { const t = y - 1800; return poly(t, &.{ 13.72, -0.332447, 0.0068612, 0.0041116, -0.00037436, 0.0000121272, -0.0000001699, 0.000000000875 }); } if (y < 1900) { const t = y - 1860; return poly(t, &.{ 7.62, 0.5737, -0.251754, 0.01680668, -0.0004473624, 1.0 / 233174.0 }); } if (y < 1920) { const t = y - 1900; return poly(t, &.{ -2.79, 1.494119, -0.0598939, 0.0061966, -0.000197 }); } if (y < 1941) { const t = y - 1920; return poly(t, &.{ 21.20, 0.84493, -0.076100, 0.0020936 }); } if (y < 1961) { const t = y - 1950; return poly(t, &.{ 29.07, 0.407, -1.0 / 233.0, 1.0 / 2547.0 }); } if (y < 1986) { const t = y - 1975; return poly(t, &.{ 45.45, 1.067, -1.0 / 260.0, -1.0 / 718.0 }); } const t = y - 2000; return poly(t, &.{ 63.86, 0.3345, -0.060374, 0.0017275, 0.000651814, 0.00002373599 }); } /// The long-term trend, −20 + 32u² with u = (y − 1820)/100. fn parabola(y: f64) f64 { const u = (y - 1820) / 100; return -20 + 32 * u * u; } /// The parabola's rate of change, in seconds a year. fn parabolaSlope(y: f64) f64 { return 64 * (y - 1820) / 10000; } /// c[0] + c[1]x + c[2]x² + …, by Horner's method. fn poly(x: f64, c: []const f64) f64 { var sum: f64 = 0; var i = c.len; while (i > 0) { i -= 1; sum = sum * x + c[i]; } return sum; }