LeapSeconds.hpp

Go to the documentation of this file.

#pragma once
 
#include <array>
#include <janus/math/Interpolate.hpp>
#include <vulcan/time/JulianDate.hpp>
#include <vulcan/time/TimeConstants.hpp>
 
namespace vulcan::time {
 
// =============================================================================
// Leap Second Table (IERS Bulletin C)
// =============================================================================

struct LeapSecondEntry {
    int year;
    int month;
    int day;
    int delta_at; 
};

inline constexpr std::array<LeapSecondEntry, 28> LEAP_SECOND_TABLE = {{
    // Year, Month, Day, TAI-UTC
    {1972, 1, 1, 10}, // Initial offset when leap seconds began
    {1972, 7, 1, 11}, {1973, 1, 1, 12}, {1974, 1, 1, 13}, {1975, 1, 1, 14},
    {1976, 1, 1, 15}, {1977, 1, 1, 16}, {1978, 1, 1, 17}, {1979, 1, 1, 18},
    {1980, 1, 1, 19}, {1981, 7, 1, 20}, {1982, 7, 1, 21}, {1983, 7, 1, 22},
    {1985, 7, 1, 23}, {1988, 1, 1, 24}, {1990, 1, 1, 25}, {1991, 1, 1, 26},
    {1992, 7, 1, 27}, {1993, 7, 1, 28}, {1994, 7, 1, 29}, {1996, 1, 1, 30},
    {1997, 7, 1, 31}, {1999, 1, 1, 32}, {2006, 1, 1, 33}, {2009, 1, 1, 34},
    {2012, 7, 1, 35}, {2015, 7, 1, 36}, {2017, 1, 1, 37}, // Most recent leap
                                                          // second
}};

[[nodiscard]] inline int leap_seconds_at_tai(double tai_jd) {
    // Before 1972: use 10 seconds as approximation
    double jd_1972 = calendar_to_jd(1972, 1, 1, 0, 0, 0.0) +
                     10.0 / constants::time::SECONDS_PER_DAY;
    if (tai_jd < jd_1972) {
        return 10;
    }
 
    // Search backwards through table (most queries will be recent dates)
    for (int i = static_cast<int>(LEAP_SECOND_TABLE.size()) - 1; i >= 0; --i) {
        const auto &entry = LEAP_SECOND_TABLE[i];
        // Convert entry date to TAI JD
        // Entry is when delta_at takes effect at 00:00:00 UTC
        // In TAI, this is 00:00:00 UTC + delta_at seconds
        double entry_tai_jd =
            calendar_to_jd(entry.year, entry.month, entry.day, 0, 0, 0.0) +
            entry.delta_at / constants::time::SECONDS_PER_DAY;
 
        if (tai_jd >= entry_tai_jd) {
            return entry.delta_at;
        }
    }
 
    return 10; // Fallback (shouldn't reach for valid inputs)
}

[[nodiscard]] inline int leap_seconds_at_utc(double utc_jd) {
    // Before 1972: use 10 seconds
    double jd_1972 = calendar_to_jd(1972, 1, 1, 0, 0, 0.0);
    if (utc_jd < jd_1972) {
        return 10;
    }
 
    // Search backwards
    for (int i = static_cast<int>(LEAP_SECOND_TABLE.size()) - 1; i >= 0; --i) {
        const auto &entry = LEAP_SECOND_TABLE[i];
        double entry_utc_jd =
            calendar_to_jd(entry.year, entry.month, entry.day, 0, 0, 0.0);
 
        if (utc_jd >= entry_utc_jd) {
            return entry.delta_at;
        }
    }
 
    return 10;
}

[[nodiscard]] inline bool is_leap_second(int year, int month, int day, int hour,
                                         int min, double sec) {
    // Leap seconds only occur at 23:59:60 on June 30 or December 31
    if (hour != 23 || min != 59 || sec < 60.0) {
        return false;
    }
 
    if (!((month == 6 && day == 30) || (month == 12 && day == 31))) {
        return false;
    }
 
    // Check if this date had a leap second
    int next_month = (month == 6) ? 7 : 1;
    int next_year = (month == 6) ? year : year + 1;
 
    for (const auto &entry : LEAP_SECOND_TABLE) {
        if (entry.year == next_year && entry.month == next_month) {
            return true;
        }
    }
 
    return false;
}

[[nodiscard]] inline double next_leap_second(double utc_jd) {
    for (const auto &entry : LEAP_SECOND_TABLE) {
        double entry_jd =
            calendar_to_jd(entry.year, entry.month, entry.day, 0, 0, 0.0);
        if (entry_jd > utc_jd) {
            return entry_jd;
        }
    }
    return 0.0; // No future leap second known
}
 
// =============================================================================
// Symbolic Leap Second Lookup (via janus::Interpolator)
// =============================================================================

[[nodiscard]] inline const janus::Interpolator &leap_second_interpolator() {
    static const janus::Interpolator interp = []() {
        janus::NumericVector jd_points(
            static_cast<int>(LEAP_SECOND_TABLE.size()));
        janus::NumericVector delta_at(
            static_cast<int>(LEAP_SECOND_TABLE.size()));
        for (size_t i = 0; i < LEAP_SECOND_TABLE.size(); ++i) {
            const auto &entry = LEAP_SECOND_TABLE[i];
            jd_points(static_cast<int>(i)) =
                calendar_to_jd(entry.year, entry.month, entry.day);
            delta_at(static_cast<int>(i)) = static_cast<double>(entry.delta_at);
        }
        return janus::Interpolator(jd_points, delta_at,
                                   janus::InterpolationMethod::Linear);
    }();
    return interp;
}

template <typename Scalar>
[[nodiscard]] Scalar leap_seconds_symbolic(const Scalar &utc_jd) {
    return leap_second_interpolator()(utc_jd);
}
 
} // namespace vulcan::time