Fix spiceInterface GPS time and Julian date

docs/source/Support/bskKnownIssues.rst

@@ -13,6 +13,14 @@ Version |release|
   ``PlanetRadiationBase``, used by :ref:`albedo` and :ref:`earthRadiationModel`, was also corrected, even though
   that case is not reached for patches on the planet surface.
 
+- GitHub issue 1612: The :ref:`spiceInterface` module computed the GPS time of ``spiceTimeOutMsg`` as the
+  elapsed TDB time since the GPS epoch instead of atomic time. ``GPSSeconds`` was in error by a periodic term of
+  a few milliseconds, and ``GPSWeek`` and ``GPSRollovers`` could be one week early at a week boundary. ``JulianDateCurrent``
+  was read from a SPICE Julian date string, which has at most seven decimals, and was therefore quantized to
+  1e-7 day. ``GPSSeconds`` overflowed a 32-bit integer after GPS week 3550, and a simulation start
+  before the GPS epoch produced invalid GPS outputs.  The corrected UTC Julian date preserves the SPICE
+  convention of holding at the following midnight during an inserted leap second, avoiding backward jumps.
+
 - GitHub issue 1609: :ref:`keplerianOrbit` returned incorrect eccentric and mean anomalies,
   lost the sign of the flight-path angle, and returned zero rates after construction or stale
   rates after ``set_a()``. These outputs are now consistent with the current elliptic orbit.

docs/source/Support/bskReleaseNotesSnippets/1612-spice-gps-time.rst

@@ -0,0 +1,2 @@
+- Fixed :ref:`spiceInterface` GPS week and seconds of ``spiceTimeOutMsg``, which were computed from TDB instead of atomic time, and its UTC Julian date, which was quantized to 1e-7 day, while preserving the SPICE convention of holding the Julian date at midnight during an inserted leap second.
+- Fixed :ref:`spiceInterface` ``GPSSeconds`` integer overflow after GPS week 3550, and held the GPS time outputs at zero with a warning when the simulation starts before the GPS epoch.

src/architecture/msgPayloadDefC/SpiceTimeMsgPayload.h

@@ -25,9 +25,9 @@
 
 //! The SPICE time output structure outputs time information to the rest of the system
 typedef struct {
-    double J2000Current;        //!< s Current J2000 elapsed time
-    double JulianDateCurrent;   //!< d Current JulianDate
-    double GPSSeconds;          //!< s Current GPS seconds
+    double J2000Current;        //!< s Current ephemeris time (TDB) elapsed since J2000
+    double JulianDateCurrent;   //!< d Current UTC Julian date
+    double GPSSeconds;          //!< s Current GPS seconds of week
     uint16_t GPSWeek;           //!< -- Current GPS week value
     uint64_t GPSRollovers;      //!< -- Count on the number of GPS rollovers
 }SpiceTimeMsgPayload;

src/simulation/environment/spiceInterface/_UnitTest/test_spiceBufferOwnership.py

@@ -24,7 +24,7 @@ def test_time_conversion_after_repeated_reset(picture, expected):
     elapsed_seconds = 1.0  # [s]
     initial_julian_date = 2451545.0  # [day]
     seconds_per_day = 86400.0  # [s/day]
-    date_tolerance = 5.0e-8  # [day] Half a unit in the final SPICE Julian-date decimal place.
+    date_tolerance = 1.0e-9  # [day] Double precision of the Julian date, with margin.
     for _ in range(3):
         model.Reset(0)
         model.UpdateState(macros.sec2nano(elapsed_seconds))
@@ -44,7 +44,7 @@ def test_time_output_rejects_insufficient_capacity(picture_length):
         model.getCurrentTimeString()
 
 
-@pytest.mark.parametrize("field", ["spiceBuffer", "charBufferSize"])
+@pytest.mark.parametrize("field", ["charBufferSize"])
 def test_scratch_storage_is_private(field):
     """Python callers cannot replace the buffer or change its capacity."""
     model = spiceInterface.SpiceInterface(auto_configure_kernels=False)

src/simulation/environment/spiceInterface/_UnitTest/test_spiceTimeOutput.py

@@ -0,0 +1,157 @@
+# ISC License
+#
+# Copyright (c) 2026, PIC4SeR & AVS Lab, Politecnico di Torino & Argotec S.R.L., University of Colorado Boulder
+#
+# Permission to use, copy, modify, and/or distribute this software for any
+# purpose with or without fee is hereby granted, provided that the above
+# copyright notice and this permission notice appear in all copies.
+#
+# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+# WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+# MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+# ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+# ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+# OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+
+"""Check the GPS time and UTC Julian date of the SPICE time output message against their definitions."""
+
+import datetime
+
+import pytest
+
+from Basilisk import __path__
+from Basilisk.simulation import spiceInterface
+from Basilisk.utilities import SimulationBaseClass, macros
+
+bsk_path = __path__[0]
+
+GPS_EPOCH = datetime.datetime(1980, 1, 6)  # GPS epoch, 1980 January 6, 00:00:00 UTC
+UNIX_EPOCH = datetime.datetime(1970, 1, 1)  # Unix epoch, 1970 January 1, 00:00:00 UTC
+UNIX_EPOCH_JULIAN_DATE = 2440587.5  # [d] Julian date of the Unix epoch
+SECONDS_PER_DAY = 86400.0  # [s]
+SECONDS_PER_WEEK = 604800.0  # [s]
+WEEKS_PER_ROLLOVER = 1024  # GPS weeks between rollovers
+
+
+def run_spice_time_output(epoch):
+    """Run ``spiceInterface`` for a single update at ``epoch`` and return the ``spiceTimeOutMsg`` recorder.
+
+    :param epoch: UTC time at which the simulation starts
+    :return: recorder of ``spiceTimeOutMsg``
+    """
+    sim = SimulationBaseClass.SimBaseClass()
+    process = sim.CreateNewProcess("process")
+    process.addTask(sim.CreateNewTask("task", macros.sec2nano(1.0)))
+
+    spice = spiceInterface.SpiceInterface()
+    spice.SPICEDataPath = bsk_path + "/supportData/EphemerisData/"
+    spice.UTCCalInit = epoch.strftime("%Y %B %d, %H:%M:%S.%f") + " (UTC)"
+    sim.AddModelToTask("task", spice)
+
+    time_log = spice.spiceTimeOutMsg.recorder()
+    sim.AddModelToTask("task", time_log)
+
+    sim.InitializeSimulation()
+    sim.ConfigureStopTime(0)
+    sim.ExecuteSimulation()
+    return time_log
+
+
+# gps_utc_offset is GPS - UTC, the leap seconds introduced since the GPS epoch, i.e. TAI - UTC - 19 s [s]
+@pytest.mark.parametrize("epoch, gps_utc_offset", [
+    (datetime.datetime(1999, 8, 21, 23, 59, 47), 13),  # GPS week 1024 starts at this UTC time
+    (datetime.datetime(2015, 7, 1, 0, 0, 0), 17),  # just after the leap second of 2015 June 30
+    (datetime.datetime(2024, 4, 4, 0, 0, 0), 18),  # TDB - TT near its maximum of about +1.7 ms
+    (datetime.datetime(2026, 3, 20, 12, 0, 0, 123000), 18),  # fractional UTC seconds
+    (datetime.datetime(2048, 6, 1, 0, 0, 0), 18),  # past GPS week 3551, where week*604800 s overflows a 32-bit int
+])
+def test_gps_time_and_julian_date(epoch, gps_utc_offset):
+    """The GPS week, GPS seconds of week and UTC Julian date of ``spiceTimeOutMsg`` equal their definitions.
+
+    GPS time is the UTC time elapsed since the GPS epoch plus the leap seconds introduced since then, and the
+    Julian date is the UTC Julian date.  The epochs include a GPS week rollover, a date just after a leap second,
+    a date near the maximum of the periodic difference between TDB and atomic time, and a date after GPS week 3551,
+    where the GPS week expressed in seconds no longer fits in a signed 32-bit integer.
+    """
+    time_log = run_spice_time_output(epoch)
+
+    gps_time = (epoch - GPS_EPOCH).total_seconds() + gps_utc_offset  # [s]
+    julian_date = UNIX_EPOCH_JULIAN_DATE + (epoch - UNIX_EPOCH).total_seconds()/SECONDS_PER_DAY  # [d]
+
+    gps_week = int(time_log.GPSWeek[-1]) + WEEKS_PER_ROLLOVER*int(time_log.GPSRollovers[-1])
+    assert gps_week == int(gps_time//SECONDS_PER_WEEK)
+    assert time_log.GPSSeconds[-1] == pytest.approx(gps_time % SECONDS_PER_WEEK, rel=0.0, abs=1.0e-6)  # [s]
+    assert time_log.JulianDateCurrent[-1] == pytest.approx(julian_date, rel=0.0, abs=1.0e-9)  # [d]
+
+
+def test_gps_time_before_gps_epoch():
+    """Before the GPS epoch, the GPS week, seconds and rollovers are held at zero."""
+    time_log = run_spice_time_output(datetime.datetime(1979, 6, 1, 0, 0, 0))
+
+    assert int(time_log.GPSWeek[-1]) == 0
+    assert int(time_log.GPSRollovers[-1]) == 0
+    assert time_log.GPSSeconds[-1] == 0.0  # [s]
+
+
+@pytest.mark.parametrize("epoch, gps_utc_offset, has_leap_second", [
+    (datetime.datetime(1998, 12, 31, 23, 59, 59), 12, True),  # before J2000
+    (datetime.datetime(2015, 6, 30, 23, 59, 59), 16, True),  # June leap second
+    (datetime.datetime(2016, 12, 31, 23, 59, 59), 17, True),  # December leap second
+    (datetime.datetime(2026, 3, 20, 23, 59, 59), 18, False),  # ordinary midnight
+])
+def test_time_output_across_midnight(epoch, gps_utc_offset, has_leap_second):
+    """Keep SPICE JDUTC fixed during an inserted second while GPS time advances.
+
+    Sample before, within and after the leap second, and also cross an ordinary
+    midnight to ensure fractional UTC seconds are preserved on normal days.
+    """
+    spice = spiceInterface.SpiceInterface()
+    spice.UTCCalInit = epoch.strftime("%Y-%m-%d %H:%M:%S") + " UTC"
+    spice.Reset(0)
+
+    midnight = epoch + datetime.timedelta(seconds=1)
+    midnight_julian_date = UNIX_EPOCH_JULIAN_DATE + (midnight - UNIX_EPOCH).total_seconds()/SECONDS_PER_DAY  # [d]
+    initial_gps_time = (epoch - GPS_EPOCH).total_seconds() + gps_utc_offset  # [s]
+    sample_times = (0.0, 0.5, 0.999, 1.0, 1.001, 1.5, 1.999, 2.0, 2.001, 2.5, 3.0)  # [s]
+    previous_julian_date = None
+    for elapsed in sample_times:
+        spice.UpdateState(macros.sec2nano(elapsed))
+        time_output = spice.spiceTimeOutMsg.read()
+        if has_leap_second:
+            seconds_from_midnight = min(elapsed - 1.0, 0.0) + max(elapsed - 2.0, 0.0)  # [s]
+        else:
+            seconds_from_midnight = elapsed - 1.0  # [s]
+        expected_julian_date = midnight_julian_date + seconds_from_midnight/SECONDS_PER_DAY  # [d]
+        expected_gps_time = initial_gps_time + elapsed  # [s]
+        gps_week = int(time_output.GPSWeek) + WEEKS_PER_ROLLOVER*int(time_output.GPSRollovers)
+
+        assert time_output.JulianDateCurrent == pytest.approx(expected_julian_date, rel=0.0, abs=1.0e-9)  # [d]
+        if has_leap_second and 1.0 <= elapsed <= 2.0:  # [s] Inserted second through the following midnight.
+            assert time_output.JulianDateCurrent == midnight_julian_date
+        if previous_julian_date is not None:
+            assert time_output.JulianDateCurrent >= previous_julian_date
+        previous_julian_date = time_output.JulianDateCurrent
+        assert gps_week == int(expected_gps_time//SECONDS_PER_WEEK)
+        assert time_output.GPSSeconds == pytest.approx(expected_gps_time % SECONDS_PER_WEEK,
+                                                     rel=0.0, abs=1.0e-6)  # [s]
+
+
+def test_gps_time_crosses_gps_epoch():
+    """Release the zero GPS outputs when a simulation crosses the GPS epoch."""
+    spice = spiceInterface.SpiceInterface()
+    spice.UTCCalInit = "1980 JAN 05 23:59:59 UTC"
+    spice.Reset(0)
+
+    sample_times = (0.0, 0.5, 1.0, 1.5, 2.0)  # [s]
+    for elapsed in sample_times:
+        spice.UpdateState(macros.sec2nano(elapsed))
+        time_output = spice.spiceTimeOutMsg.read()
+        expected_seconds = max(elapsed - 1.0, 0.0)  # [s]
+        assert time_output.GPSWeek == 0
+        assert time_output.GPSRollovers == 0
+        assert time_output.GPSSeconds == pytest.approx(expected_seconds, rel=0.0, abs=1.0e-6)  # [s]
+
+
+if __name__ == "__main__":
+    test_gps_time_and_julian_date(datetime.datetime(1999, 8, 21, 23, 59, 47), 13)

src/simulation/environment/spiceInterface/_UnitTest/test_unitSpice.py

@@ -250,6 +250,11 @@ def unitSpice(testPlottingFixture, show_plots, DateSpice, DatePlot, MarsTruthPos
         GPSEndTime = GPSTimeSeconds + 17 + 60.0 - 68.184  # 17 GPS skip seconds passed that date
     else:
         GPSEndTime = GPSTimeSeconds + 16 + 60.0 - 67.184  # 16 GPS skip seconds before
+    # GPS time is atomic time, so remove the periodic TDB - TT term of the SPICE leap seconds kernel
+    ETEnd = (date - datetime.datetime(2000, 1, 1, 12)).total_seconds() + 60.0  # [s] TDB past J2000
+    meanAnomaly = 6.239996 + 1.99096871E-7*ETEnd  # [rad] Earth-Moon barycenter mean anomaly
+    eccentricAnomaly = meanAnomaly + 1.671E-2*numpy.sin(meanAnomaly)  # [rad]
+    GPSEndTime -= 1.657E-3*numpy.sin(eccentricAnomaly)  # [s]
 
     GPSWeek = int(GPSEndTime / (86400 * 7))
     GPSSecondAssumed = GPSEndTime - GPSWeek * 86400 * 7

src/simulation/environment/spiceInterface/spiceInterface.cpp

@@ -34,6 +34,9 @@
 #include "spiceInterface.h"
 
 namespace {
+    constexpr double julianDateJ2000 = 2451545.0; //!< [d] Julian date of the J2000 epoch
+    constexpr double secondsPerDay = 86400.0;     //!< [s] Seconds in a Julian day
+
     /**
      * RAII guard for SPICE error mode.
      *
@@ -584,7 +587,7 @@ SpiceInterface::SpiceInterface()
     GPSWeek = 0;
     GPSRollovers = 0;
     timeDataInit = false;
-    JDGPSEpoch = 0.0;
+    TAIGPSEpoch = 0.0;
     GPSEpochTime = "1980 January 6, 00:00:00.0";
 
     referenceBase = "j2000";
@@ -716,8 +719,6 @@ void SpiceInterface::Reset(uint64_t CurrenSimNanos)
  */
 void SpiceInterface::initTimeData()
 {
-    double EpochDelteET;
-
     /* set epoch information.  If provided, then the epoch message information should be used.  */
     if (this->epochInMsg.isLinked()) {
         // Read in the epoch message and set the internal time structure
@@ -733,34 +734,38 @@ void SpiceInterface::initTimeData()
         }
     }
 
-    //! -Get the time value associated with the GPS epoch
-    str2et_c(this->GPSEpochTime.c_str(), &this->JDGPSEpoch);
+    //! - Get the TAI of the GPS epoch, which is fixed, so the conversion is done once here
+    double etGPSEpoch;
+    str2et_c(this->GPSEpochTime.c_str(), &etGPSEpoch);
+    this->TAIGPSEpoch = unitim_c(etGPSEpoch, "TDB", "TAI"); // [s]
     //! - Get the time value associate with the requested UTC date
     str2et_c(this->UTCCalInit.c_str(), &this->J2000ETInit);
-    //! - Take the JD epoch and get the elapsed time for it
-    deltet_c(this->JDGPSEpoch, "ET", &EpochDelteET);
-
+    if (this->J2000ETInit < etGPSEpoch) {
+        bskLogger.bskLog(BSK_WARNING, "The initial time is before the GPS epoch.  "
+                                      "The GPS time outputs are held at zero until the GPS epoch.");
+    }
 }
 
-/*! This method computes the GPS time data for the current elapsed time.  It uses
- the total elapsed times at both the GPS epoch time and the current time to
- compute the GPS time (week, seconds, rollovers)
+/*! This method computes the GPS time data for the current elapsed time.  GPS time
+ is atomic time (TAI - 19 s), so the elapsed time since the GPS epoch is the
+ difference of the TAI times at the current time and at the GPS epoch.  This
+ elapsed time gives the GPS time (week, seconds, rollovers).
 
  */
 void SpiceInterface::computeGPSData()
 {
-    double JDDifference;
-
-    //! - The difference between the epochs in julian date terms is the total
-    JDDifference = this->J2000Current - this->JDGPSEpoch;
-    //! - Scale the elapsed by a week's worth of seconds to get week
-    this->GPSWeek = (uint16_t) (JDDifference/(7*86400));
-    //! - Subtract out the GPS week scaled up to seconds to get time in week
-    this->GPSSeconds = JDDifference - this->GPSWeek*7*86400;
-
-    //! - Maximum GPS week is 1024 so get rollovers and subtract out those weeks
-    this->GPSRollovers = this->GPSWeek/1024;
-    this->GPSWeek = (uint16_t)(this->GPSWeek-this->GPSRollovers*1024);
+    //! - Take the elapsed atomic time since the GPS epoch, held at zero before the epoch
+    double gpsElapsed = std::max(unitim_c(this->J2000Current, "TDB", "TAI") - this->TAIGPSEpoch, 0.0); // [s]
+    const double secondsPerWeek = 604800.0; // [s]
+    const uint64_t weeksPerRollover = 1024;
+    //! - Scale the elapsed by a week's worth of seconds to get the total number of weeks
+    uint64_t totalWeeks = static_cast<uint64_t>(gpsElapsed/secondsPerWeek);
+    //! - Subtract out the weeks scaled up to seconds to get time in week, in floating point to avoid integer overflow
+    this->GPSSeconds = gpsElapsed - static_cast<double>(totalWeeks)*secondsPerWeek; // [s]
+
+    //! - Maximum GPS week is 1024 so get the rollovers and the week modulo 1024
+    this->GPSRollovers = totalWeeks/weeksPerRollover;
+    this->GPSWeek = static_cast<uint16_t>(totalWeeks % weeksPerRollover);
 }
 
 /*! This method takes the values computed in the model and outputs them.
@@ -821,10 +826,19 @@ void SpiceInterface::UpdateState(uint64_t CurrentSimNanos)
     //! - Increment the J2000 elapsed time based on init value and Current sim
     this->J2000Current = this->J2000ETInit + static_cast<double>(CurrentSimNanos) * NANO2SEC;
 
-    //! - Compute the current Julian Date string and cast it over to the double
-    et2utc_c(this->J2000Current, "J", 14, this->charBufferSize - 1, this->spiceBuffer.data());
-    std::string localString = this->spiceBuffer.data() + 3;
-    this->julianDateCurrent = std::stod(localString);
+    //! - Compute the current UTC Julian date from the UTC seconds past J2000
+    double deltaET; // [s] ET - UTC
+    deltet_c(this->J2000Current, "ET", &deltaET);
+    double utcSeconds = this->J2000Current - deltaET; // [s] UTC seconds past J2000
+    double utcDeltaET; // [s] ET - UTC when interpreting utcSeconds as a UTC epoch
+    deltet_c(utcSeconds, "UTC", &utcDeltaET);
+    const double leapSecondThreshold = 0.5; // [s] Distinguish a leap from roundoff in the periodic ET - TAI term.
+    if (utcDeltaET - deltaET > leapSecondThreshold) {
+        // During an inserted second, utcSeconds falls in the first second of the next day,
+        // whose UTC offset is larger. Preserve SPICE JDUTC by holding at that midnight.
+        utcSeconds = std::floor(utcSeconds);
+    }
+    this->julianDateCurrent = julianDateJ2000 + utcSeconds / secondsPerDay; // [d]
     //! Get GPS and Planet data and then write the message outputs. Reconstruction is allowed only
     //! on the planet pass (the raw sim time drives the knot grid; the spacecraft pass stays exact).
     this->computeGPSData();

src/simulation/environment/spiceInterface/spiceInterface.h

@@ -230,19 +230,18 @@ class SpiceInterface: public SysModel {
     bool timeDataInit;          //!< -- Flag indicating whether time has been init
     double J2000ETInit;         //!< s Seconds elapsed since J2000 at init
     double J2000Current;        //!< s Current J2000 elapsed time
-    double julianDateCurrent;   //!< s Current JulianDate
-    double GPSSeconds;          //!< s Current GPS seconds
+    double julianDateCurrent;   //!< d Current UTC Julian date
+    double GPSSeconds;          //!< s Current GPS seconds of week
     uint16_t GPSWeek;           //!< -- Current GPS week value
     uint64_t GPSRollovers;      //!< -- Count on the number of GPS rollovers
 
     BSKLogger bskLogger;                      //!< -- BSK Logging
 
 private:
     static constexpr int charBufferSize = 512; //!< Capacity of SPICE scratch buffers, including the null terminator.
-    std::array<char, charBufferSize> spiceBuffer{}; //!< Scratch storage for Julian date conversion.
 
     std::string GPSEpochTime;   //!< -- String for the GPS epoch
-    double JDGPSEpoch;          //!< s Epoch for GPS time.  Saved for efficiency
+    double TAIGPSEpoch;         //!< s TAI of the GPS epoch, seconds past J2000.  Saved for efficiency
 
     std::vector<SpicePlanetStateMsgPayload> planetData;
     std::vector<SpicePlanetStateMsgPayload> scData;