Leap Seconds V: UTC Strikes Back

src/AstroDates.jl

@@ -7,7 +7,8 @@ import Dates
 
 export Date, Time, DateTime,
     year, month, day, j2000day, calendar,
-    hour, minute, second, date, time
+    hour, minute, second, date, time,
+    secondinday
 
 abstract type Calendar end
 
@@ -187,6 +188,7 @@ const GPS_EPOCH = Date(1980, 1, 6)
 const J2000_EPOCH = Date(2000, 1, 1)
 const MIN_EPOCH = Date(typemin(Int32))
 const MAX_EPOCH = Date(typemax(Int32))
+const UNIX_EPOCH = Date(1970, 1, 1)
 
 struct Time
     hour::Int
@@ -236,6 +238,8 @@ hour(s::Time) = s.hour
 minute(s::Time) = s.minute
 second(s::Time) = s.second
 
+secondinday(t::Time) = t.second + 60 * t.minute + 3600 * t.hour
+
 function show(io::IO, t::Time)
     h = lpad(hour(t), 2, '0')
     m = lpad(minute(t), 2, '0')

src/AstroTime.jl

@@ -58,6 +58,8 @@ macro timescale(scale::Symbol, ep::Symbol, args...)
     return quote
         struct $sc <: TimeScale end
         const $epoch = Epoch{$sc()}
+        Base.show(io::IO, $sc) = print(io, string(typeof($sc)))
+
         Dates.CONVERSION_TRANSLATIONS[$epoch] = (
             Dates.Year,
             Dates.Month,

src/Epochs.jl

@@ -12,7 +12,8 @@ using ..Periods
 
 export Epoch,
     JULIAN_EPOCH, J2000_EPOCH, MODIFIED_JULIAN_EPOCH,
-    FIFTIES_EPOCH, GALILEO_EPOCH, GPS_EPOCH, CCSDS_EPOCH
+    FIFTIES_EPOCH, GALILEO_EPOCH, GPS_EPOCH, CCSDS_EPOCH,
+    PAST_INFINITY, FUTURE_INFINITY, UNIX_EPOCH
 
 struct Epoch{S, T} <: Dates.AbstractDateTime
     epoch::Int64
@@ -77,8 +78,6 @@ Epoch{S}(d::Date) where {S} = Epoch{S}(d, AstroDates.H00)
 
 Epoch{S}(dt::DateTime) where {S} = Epoch{S}(date(dt), time(dt))
 
-include("leapseconds.jl")
-
 function Epoch{S}(year::Int, month::Int, day::Int, hour::Int=0,
                   minute::Int=0, second::Float64=0.0) where S
     Epoch{S}(Date(year, month, day), Time(hour, minute, second))
@@ -109,13 +108,19 @@ isless(ep1::Epoch, ep2::Epoch) = isless(get(ep1 - ep2), 0.0)
 -(ep::Epoch{S}, p::Period) where {S} = Epoch{S}(ep, -get(seconds(p)))
 -(a::Epoch, b::Epoch) = ((a.epoch - b.epoch) + (a.offset - b.offset)) * seconds
 
-const JULIAN_EPOCH = Epoch{TT}(AstroDates.JULIAN_EPOCH, AstroDates.H12)
-const J2000_EPOCH = Epoch{TT}(AstroDates.J2000_EPOCH, AstroDates.H12)
-const MODIFIED_JULIAN_EPOCH = Epoch{TT}(AstroDates.MODIFIED_JULIAN_EPOCH, AstroDates.H00)
-const FIFTIES_EPOCH = Epoch{TT}(AstroDates.FIFTIES_EPOCH, AstroDates.H00)
-const CCSDS_EPOCH = Epoch{TT}(AstroDates.CCSDS_EPOCH, AstroDates.H00)
-const GALILEO_EPOCH = Epoch{TT}(AstroDates.GALILEO_EPOCH, AstroDates.H00)
-const GPS_EPOCH = Epoch{TT}(AstroDates.GPS_EPOCH, AstroDates.H00)
+include("leapseconds.jl")
+
+const JULIAN_EPOCH = TTEpoch(AstroDates.JULIAN_EPOCH, AstroDates.H12)
+const J2000_EPOCH = TTEpoch(AstroDates.J2000_EPOCH, AstroDates.H12)
+const MODIFIED_JULIAN_EPOCH = TTEpoch(AstroDates.MODIFIED_JULIAN_EPOCH, AstroDates.H00)
+const FIFTIES_EPOCH = TTEpoch(AstroDates.FIFTIES_EPOCH, AstroDates.H00)
+const CCSDS_EPOCH = TTEpoch(AstroDates.CCSDS_EPOCH, AstroDates.H00)
+const GALILEO_EPOCH = TTEpoch(AstroDates.GALILEO_EPOCH, AstroDates.H00)
+const GPS_EPOCH = TTEpoch(AstroDates.GPS_EPOCH, AstroDates.H00)
+const UNIX_EPOCH = TAIEpoch(AstroDates.UNIX_EPOCH, Time(0, 0, 10.0))
+
+const PAST_INFINITY = TAIEpoch(UNIX_EPOCH, -Inf)
+const FUTURE_INFINITY = TAIEpoch(UNIX_EPOCH, Inf)
 
 const EPOCH_77 = Epoch{TAI}(1977, 1, 1)
 

src/TimeScales.jl

@@ -42,9 +42,9 @@ for (acronym, scale) in zip(ACRONYMS, SCALES)
         struct $scale <: TimeScale end
         const $acronym = $scale()
         export $scale, $acronym
+
+        Base.show(io::IO, ::$scale) = print(io, $name)
     end
 end
 
-Base.show(io::IO, x::TimeScale) = print(io, string(typeof(x)))
-
 end

src/accessors.jl

@@ -12,9 +12,9 @@ j2000(ep::Epoch) = get(days(ep - J2000_EPOCH))
 function DateTime(ep::Epoch)
     if !isfinite(ep.offset)
         if ep.offset < 0
-            DateTime(AstroDates.MIN_EPOCH, AstroDates.H00)
+            return DateTime(AstroDates.MIN_EPOCH, AstroDates.H00)
         else
-            DateTime(AstroDates.MAX_EPOCH, Time(23, 59, 59.999))
+            return DateTime(AstroDates.MAX_EPOCH, Time(23, 59, 59.999))
         end
     end
 
@@ -43,7 +43,7 @@ function DateTime(ep::Epoch)
     time_comp = Time(time, offset2000B)
 
     if insideleap(ep)
-        leap = 1.0
+        leap = getleap(ep)
         time_comp = Time(hour(time_comp), minute(time_comp), second(time_comp) + leap)
     end
 

src/leapseconds.jl

@@ -3,35 +3,56 @@ using LeapSeconds
 
 export insideleap
 
-struct TaiOffset
+insideleap(ep::Epoch{S}) where {S} = false
+getleap(ep::Epoch{S}) where {S} = 0.0
+
+struct TAIOffset
+    date::TAIEpoch{Float64}
+    jd::Float64
     start::TAIEpoch{Float64}
-    stop::TAIEpoch{Float64}
+    reference::TAIEpoch{Float64}
     leap::Float64
-    epoch::Float64
-    drift::Float64
     offset::Float64
+    slope_utc::Float64
+    slope_tai::Float64
 end
 
-isless(ep::Epoch, t::TaiOffset) = isless(ep, t.start)
+Base.isless(t::TAIOffset, ep::Epoch) = isless(ep, t.date)
+Base.isless(ep::Epoch, t::TAIOffset) = isless(ep, t.date)
 
-const OFFSETS = TaiOffset[]
+Base.isless(t::TAIOffset, jd::Float64) = isless(jd, t.jd)
+Base.isless(jd::Float64, t::TAIOffset) = isless(jd, t.jd)
 
-for (i, (epoch, offset, dep, rate)) in enumerate(zip(LeapSeconds.EPOCHS,
-                                                     LeapSeconds.OFFSETS,
-                                                     LeapSeconds.DRIFT_EPOCHS,
-                                                     LeapSeconds.DRIFT_RATES))
-    tai = TAIEpoch(Dates.julian2datetime(epoch))
-    previous = i == 1 ? 0.0 : OFFSETS[i-1].offset
-    start = TAIEpoch(tai, previous)
-    start_offset = offset + (epoch - dep) * rate
-    slope = rate / SECONDS_PER_DAY
+const TAI_OFFSETS = TAIOffset[]
+
+function findoffset(ep)
+    idx = searchsortedlast(TAI_OFFSETS, ep)
+    idx == 0 && return nothing
+
+    TAI_OFFSETS[idx]
 end
 
-for (epoch, offset) in zip(LeapSeconds.LS_EPOCHS, LeapSeconds.LEAP_SECONDS)
-    tai = TAIEpoch(Dates.julian2datetime(epoch))
-    push!(OFFSETS,
-          TaiOffset(TAIEpoch(tai, offset - 1), TAIEpoch(tai, offset),
-                    1.0, 0.0, 0.0, 0.0))
+const EPOCHS = [LeapSeconds.EPOCHS; LeapSeconds.LS_EPOCHS]
+const OFFSETS = [LeapSeconds.OFFSETS; LeapSeconds.LEAP_SECONDS]
+const DRIFT_EPOCHS = [LeapSeconds.DRIFT_EPOCHS;
+                      zeros(length(LeapSeconds.LS_EPOCHS))]
+const DRIFT_RATES = [LeapSeconds.DRIFT_RATES;
+                     zeros(length(LeapSeconds.LS_EPOCHS))]
+
+for (ep, offset, dep, rate) in zip(EPOCHS, OFFSETS, DRIFT_EPOCHS, DRIFT_RATES)
+    previous = isempty(TAI_OFFSETS) ? 0.0 : last(TAI_OFFSETS).offset
+    tai = TAIEpoch(Dates.julian2datetime(ep))
+    ref = TAIEpoch(TAIEpoch(Dates.julian2datetime(dep)), offset)
+    start = TAIEpoch(tai, previous)
+    start_offset = offset + (ep - dep) * rate
+    stop = TAIEpoch(tai, start_offset)
+    slope = rate / SECONDS_PER_DAY
+    leap = get(stop - start) / (1 + slope)
+    o = TAIOffset(start, ep,
+                  TAIEpoch(start, leap),
+                  ref,
+                  leap, offset, slope, slope / (1 + slope))
+    push!(TAI_OFFSETS, o)
 end
 
 const LEAP_STARTS = TAIEpoch{Float64}[]
@@ -44,12 +65,27 @@ for (epoch, offset) in zip(LeapSeconds.LS_EPOCHS, LeapSeconds.LEAP_SECONDS)
 end
 
 function insideleap(ep::UTCEpoch)
-    idx = searchsortedlast(LEAP_STARTS, ep)
-    if idx == 0
-        return 0.0
-    else
-        ep < LEAP_STOPS[idx]
-    end
+    offset = findoffset(ep)
+    offset === nothing && return false
+
+    ep < offset.start
 end
 
-insideleap(ep::Epoch{S}) where {S} = false
+function getleap(ep::UTCEpoch)
+    offset = findoffset(ep)
+    offset === nothing && return 0.0
+
+    offset.leap
+end
+
+function getoffset(t::TAIOffset, ep::Epoch)
+    t.slope_tai == 0.0 && return t.offset
+
+    t.offset + get(ep - t.reference) * t.slope_tai
+end
+
+function getoffset(to::TAIOffset, d::Date, t::Time)
+    days = AstroDates.julian(d) - to.jd
+    fraction = secondinday(t)
+    to.offset + days * (to.slope_utc * SECONDS_PER_DAY) + fraction * to.slope_utc
+end

src/offsets.jl

@@ -10,11 +10,17 @@ const LB_RATE = 1.550519768e-8
 
 tai_offset(::InternationalAtomicTime, ep) = 0.0
 tai_offset(::TerrestrialTime, ep) = OFFSET_TAI_TT
-tai_offset(::CoordinatedUniversalTime, ep) = -offset_tai_utc(julian(ep))
 tai_offset(::UniversalTime, ep) = tai_offset(UTC, ep) + getΔUT1(julian(ep))
 tai_offset(::GeocentricCoordinateTime, ep) = tai_offset(TT, ep) + LG_RATE * get(ep - EPOCH_77)
 tai_offset(::BarycentricCoordinateTime, ep) = tai_offset(TT, ep) + LB_RATE * get(ep - EPOCH_77)
 
+function tai_offset(::CoordinatedUniversalTime, ep)
+    offset = findoffset(ep)
+    offset === nothing && return 0.0
+
+    -getoffset(offset, ep)
+end
+
 """
     tai_offset(TDB, ep)
 
@@ -120,7 +126,10 @@ tai_offset(::InternationalAtomicTime, date, time) = 0.0
 function tai_offset(::CoordinatedUniversalTime, date, time)
     minute_in_day = hour(time) * 60 + minute(time)
     correction  = minute_in_day < 0 ? (minute_in_day - 1439) ÷ 1440 : minute_in_day ÷ 1440;
-    offset_tai_utc(AstroDates.julian(date) + correction)
+    offset = findoffset(AstroDates.julian(date) + correction)
+    offset === nothing && return 0.0
+
+    getoffset(offset, date, time)
 end
 
 function tai_offset(scale, date, time)

test/epochs.jl

@@ -17,5 +17,39 @@
     tai1 = Epoch{TAI}(tt)
     @test tai1.epoch == 100000
     @test tai1.offset ≈ 0.0
+
+    @testset "Leap Seconds" begin
+        @test string(UTCEpoch(2018, 8, 8, 0, 0, 0.0)) == "2018-08-08T00:00:00.000 UTC"
+
+        # Test transformation to calendar date during pre-leap second era
+        @test string(UTCEpoch(1961, 3, 5, 23, 4, 12.0)) == "1961-03-05T23:04:12.000 UTC"
+
+        let
+            before = UTCEpoch(2012, 6, 30, 23, 59, 59.0)
+            start = UTCEpoch(2012, 6, 30, 23, 59, 60.0)
+            during = UTCEpoch(2012, 6, 30, 23, 59, 60.5)
+            after = UTCEpoch(2012, 7, 1, 0, 0, 0.0)
+
+            @test before.epoch == 394372833
+            @test before.offset == 0.0
+            @test start.epoch == 394372834
+            @test start.offset == 0.0
+            @test during.epoch == 394372834
+            @test during.offset == 0.5
+            @test after.epoch == 394372835
+            @test after.offset == 0.0
+
+            @test !insideleap(before)
+            @test insideleap(start)
+            @test insideleap(during)
+            @test !insideleap(after)
+
+            # Test transformation to calendar date during leap second
+            @test string(before) == "2012-06-30T23:59:59.000 UTC"
+            @test string(start) == "2012-06-30T23:59:60.000 UTC"
+            @test string(during) == "2012-06-30T23:59:60.500 UTC"
+            @test string(after) == "2012-07-01T00:00:00.000 UTC"
+        end
+    end
 end