Redesign of AstroPeriod née Period and other small refactorings

.github/workflows/CI.yml

@@ -2,10 +2,10 @@ name: CI
 on:
   pull_request:
     branches:
-      - master
+      - main
   push:
     branches:
-      - master
+      - main
     tags: '*'
 jobs:
   test:
@@ -15,7 +15,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
+          - '1.3' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
           - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
           - 'nightly'
         os:

Project.toml

@@ -20,7 +20,7 @@ MacroTools = "0.5"
 Measurements = "2.2"
 MuladdMacro = "0.2"
 Reexport = "0.2, 1"
-julia = "1"
+julia = "1.3"
 
 [extras]
 ERFA = "17511681-8477-586a-8d98-4cfd5a1f2ec3"

docs/src/tutorial.md

@@ -83,7 +83,7 @@ julia> AstroTime.format(ep, "dd.mm.yyyy HH:MM ttt")
 
 ## Working with Epochs and Periods
 
-You can shift an `Epoch` in time by adding or subtracting a [`Period`](@ref) to it.
+You can shift an `Epoch` in time by adding or subtracting an [`AstroPeriod`](@ref) to it.
 
 AstroTime.jl provides a convenient way to construct periods by multiplying a value
 with a time unit.
@@ -105,7 +105,7 @@ The following time units are available:
 - `years`
 - `centuries`
 
-To shift an `Epoch` forward in time add a `Period` to it.
+To shift an `Epoch` forward in time add an `AstroPeriod` to it.
 
 ```julia
 julia> ep = UTCEpoch(2000, 1, 1)
@@ -125,7 +125,7 @@ julia> ep - 1days
 1999-12-31T00:00:00.000 UTC
 ```
 
-If you subtract two epochs you will receive the time between them as a `Period`.
+If you subtract two epochs you will receive the time between them as an `AstroPeriod`.
 
 ```julia
 julia> ep1 = UTCEpoch(2000, 1, 1)
@@ -138,8 +138,8 @@ julia> ep2 - ep1
 86400.0 seconds
 ```
 
-You can also construct a `Period` with a different time unit from
-another `Period`.
+You can also construct an `AstroPeriod` with a different time unit from
+another `AstroPeriod`.
 
 ```julia
 julia> dt = 86400.0seconds

src/AccurateArithmetic.jl

@@ -0,0 +1,57 @@
+module AccurateArithmetic
+
+# Adapted from AccurateArithmetic.jl
+
+export two_sum, apply_offset
+
+@inline function two_sum(a, b)
+    hi = a + b
+    a1 = hi - b
+    b1 = hi - a1
+    lo = (a - a1) + (b - b1)
+    return hi, lo
+end
+
+@inline function two_hilo_sum(a, b)
+    hi = a + b
+    lo = b - (hi - a)
+    return hi, lo
+end
+
+function four_sum(a, b, c, d)
+    t0, t1 = two_sum(a,  b)
+    t2, t3 = two_sum(c,  d)
+    hi, t4 = two_sum(t0, t2)
+    t5, lo = two_sum(t1, t3)
+    hm, ml = two_sum(t4, t5)
+    ml, lo = two_hilo_sum(ml, lo)
+    hm, ml = two_hilo_sum(hm, ml)
+    hi, hm = two_hilo_sum(hi,hm)
+    return hi, hm, ml, lo
+end
+
+function handle_infinity(fraction)
+    second = ifelse(fraction < 0, typemin(Int64), typemax(Int64))
+    return second, fraction, zero(fraction)
+end
+
+function apply_offset(s1::Int64, f1, e1, s2::Int64, f2, e2)
+    isfinite(f1 + f2) || return handle_infinity(f1 + f2)
+
+    sum, residual, _ = four_sum(f1, f2, e1, e2)
+    int_seconds = floor(Int64, sum)
+    second = s1 + s2 + int_seconds
+    fraction = sum - int_seconds
+    return second, fraction, residual
+end
+
+function apply_offset(s1::Int64, f1, e1, offset)
+    isfinite(f1 + offset) || return handle_infinity(f1 + offset)
+
+    s2 = floor(Int64, offset)
+    f2 = offset - s2
+    return apply_offset(s1, f1, e1, s2, f2, 0.0)
+end
+
+end
+

src/AstroTime.jl

@@ -40,6 +40,7 @@ const EPOCH_ISO_FORMAT = Ref{Dates.DateFormat{Symbol("yyyy-mm-ddTHH:MM:SS.fff tt
                                                     Dates.Delim{Char, 1},
                                                     Dates.DatePart{'t'}}}}()
 
+include("AccurateArithmetic.jl")
 include("TimeScales.jl")
 include("Periods.jl")
 include("AstroDates.jl")
@@ -235,7 +236,6 @@ macro timescale(scale::Symbol, parent=nothing, oneway=false)
             Dates.Second,
             Dates.Millisecond,
         )
-        Dates.default_format(::Type{$epoch_type}) = Dates.default_format(AstroDates.DateTime)
 
         $reg_expr
 

src/Epochs/Epochs.jl

@@ -18,6 +18,7 @@ import ..EPOCH_ISO_FORMAT
 import ..AstroDates: Date, DateTime, Time
 import ..AstroDates: calendar, fractionofday, fractionofsecond, secondinday, subsecond
 import ..AstroDates: j2000, julian, julian_twopart
+import ..AccurateArithmetic: apply_offset, two_sum
 
 using ..TimeScales: find_path
 

src/Epochs/aliases.jl

@@ -32,7 +32,7 @@ for (scale, acronym) in zip(TimeScales.NAMES, TimeScales.ACRONYMS)
         $epoch(::AbstractString)
 
         """
-            $($name)(jd1::T, jd2::T=zero(T); origin=:j2000) where T<:Period
+            $($name)(jd1::T, jd2::T=zero(T); origin=:j2000) where T<:AstroPeriod
 
         Construct a $($name) from a Julian date (optionally split into
         `jd1` and `jd2`). `origin` determines the variant of Julian

src/Epochs/dates.jl

@@ -7,7 +7,7 @@ function Epoch{S}(date::Date, time::Time, args...) where S
     return Epoch{S}(sec, time.fraction)
 end
 
-Dates.default_format(::Type{Epoch}) = EPOCH_ISO_FORMAT[]
+Dates.default_format(::Type{<:Epoch}) = EPOCH_ISO_FORMAT[]
 
 """
     Epoch(str[, format])
@@ -35,8 +35,6 @@ Epoch(str::AbstractString, format::Dates.DateFormat=Dates.default_format(Epoch))
 
 Epoch(str::AbstractString, format::AbstractString) = Epoch(str, Dates.DateFormat(format))
 
-Dates.default_format(::Type{Epoch{S}}) where {S} = Dates.default_format(AstroDates.DateTime)
-
 """
     Epoch{S}(str[, format]) where S
 

src/Epochs/julian.jl

@@ -2,7 +2,7 @@ const J2000_TO_JULIAN = 2.451545e6days
 const J2000_TO_MJD = 51544.5days
 
 """
-    Epoch{S}(jd1::T, jd2::T=zero(T); origin=:j2000) where {S, T<:Period}
+    Epoch{S}(jd1::T, jd2::T=zero(T); origin=:j2000) where {S, T<:AstroPeriod}
 
 Construct an `Epoch` with time scale `S` from a Julian date
 (optionally split into `jd1` and `jd2`). `origin` determines the
@@ -22,7 +22,7 @@ julia> Epoch{InternationalAtomicTime}(2.451545e6days, origin=:julian)
 2000-01-01T12:00:00.000 TAI
 ```
 """
-function Epoch{S}(jd1::T, jd2::T=zero(T), args...; origin=:j2000) where {S, T<:Period}
+function Epoch{S}(jd1::T, jd2::T=zero(T), args...; origin=:j2000) where {S, T<:AstroPeriod}
     if jd2 > jd1
         jd1, jd2 = jd2, jd1
     end
@@ -52,7 +52,7 @@ end
     julian_period([T,] ep::Epoch; origin=:j2000, scale=timescale(ep), unit=days)
 
 Return the period since Julian Epoch `origin` within the time scale `scale` expressed in
-`unit` for a given epoch `ep`. The result is a [`Period`](@ref) object by default.
+`unit` for a given epoch `ep`. The result is an [`AstroPeriod`](@ref) object by default.
 If the type argument `T` is present, the result is converted to `T` instead.
 
 ### Example ###

src/Epochs/offsets.jl

@@ -15,7 +15,7 @@ julia> TTEpoch(32.184, ep)
 ```
 """
 function Epoch{S2}(offset, ep::Epoch{S1}) where {S1<:TimeScale, S2<:TimeScale}
-    second, fraction = apply_offset(ep.second, ep.fraction, offset)
+    second, fraction = apply_offset(ep.second, ep.fraction, ep.error, offset)
     Epoch{S2}(second, fraction)
 end
 
@@ -36,8 +36,8 @@ julia> TAIEpoch(ep)
 ```
 """
 function Epoch{S2}(ep::Epoch{S1}) where {S1<:TimeScale, S2<:TimeScale}
-    second, fraction = apply_offset(ep.second, ep.fraction, S1(), S2())
-    Epoch{S2}(second, fraction)
+    second, fraction, error = apply_offset(ep.second, ep.fraction, ep.error, S1(), S2())
+    Epoch{S2}(second, fraction, error)
 end
 
 """
@@ -57,14 +57,14 @@ julia> Epoch(ep, TAI)
 ```
 """
 function Epoch(ep::Epoch{S1}, scale::S2) where {S1<:TimeScale, S2<:TimeScale}
-    second, fraction = apply_offset(ep.second, ep.fraction, S1(), S2())
-    Epoch{S2}(second, fraction)
+    second, fraction, error = apply_offset(ep.second, ep.fraction, ep.error, S1(), S2())
+    Epoch{S2}(second, fraction, error)
 end
 
 function Epoch{S2}(ep::Epoch{S1}, args...) where {S1<:TimeScale, S2<:TimeScale}
     offset = getoffset(S1(), S2(), ep.second, ep.fraction, args...)
-    second, fraction = apply_offset(ep.second, ep.fraction, offset)
-    Epoch{S2}(second, fraction)
+    second, fraction, error = apply_offset(ep.second, ep.fraction, ep.error, offset)
+    Epoch{S2}(second, fraction, error)
 end
 
 Epoch{S}(ep::Epoch{S}) where {S<:TimeScale} = ep
@@ -115,10 +115,11 @@ function getoffset(ep::Epoch{S}, scale::TimeScale) where S<:TimeScale
     total_offset = 0.0
     second = ep.second
     fraction = ep.fraction
+    error = ep.error
     for i in 1:length(path) - 1
         offset::Float64 = getoffset(path[i], path[i+1], second, fraction)
         total_offset += offset
-        second, fraction = apply_offset(second, fraction, offset)
+        second, fraction, error = apply_offset(second, fraction, error, offset)
     end
     return total_offset
 end
@@ -146,19 +147,21 @@ end
 
 @inline function apply_offset(second::Int64,
                               fraction::T,
+                              error::T,
                               from::S1,
-                              to::S2)::Tuple{Int64, T} where {T, S1<:TimeScale, S2<:TimeScale}
+                              to::S2)::Tuple{Int64, T, T} where {T, S1<:TimeScale, S2<:TimeScale}
     path = find_path(from, to)
     isempty(path) && throw(NoPathError(string(from), string(to)))
-    length(path) == 2 && return _apply_offset(second, fraction, from, to)
-    return _apply_offset((second, fraction), path...)
+    length(path) == 2 && return _apply_offset(second, fraction, error, from, to)
+    return _apply_offset((second, fraction, error), path...)
 end
 
 @inline function _apply_offset(second::Int64,
                                fraction::T,
+                               error::T,
                                from::S1,
-                               to::S2)::Tuple{Int64, T} where {T, S1<:TimeScale, S2<:TimeScale}
-    return apply_offset(second, fraction, getoffset(from, to, second, fraction))
+                               to::S2)::Tuple{Int64, T, T} where {T, S1<:TimeScale, S2<:TimeScale}
+    return apply_offset(second, fraction, error, getoffset(from, to, second, fraction))
 end
 
 @generated function _apply_offset(sf, path...)
@@ -554,7 +557,7 @@ julia> getoffset(TT, TDB, 0, 0.0, π, 6371.0, 0.0)
     offset = 0.0
     for _ in 1:3
         offset = -getoffset(TDB, TT, tdb1, tdb2, elong, u, v)
-        tdb1, tdb2 = apply_offset(tt1, tt2, offset)
+        tdb1, tdb2 = apply_offset(tt1, tt2, 0.0, offset)
     end
     return offset
 end

src/Epochs/operations.jl

@@ -8,8 +8,11 @@ end
 
 Base.isless(ep1::Epoch, ep2::Epoch) = isless(value(ep1 - ep2), 0.0)
 
-Base.:+(ep::Epoch{S}, p::Period) where {S} = Epoch{S}(ep, value(seconds(p)))
-Base.:-(ep::Epoch{S}, p::Period) where {S} = Epoch{S}(ep, -value(seconds(p)))
+function Base.:+(ep::Epoch{S}, p::AstroPeriod) where {S}
+    second, fraction, error = apply_offset(ep.second, ep.fraction, ep.error, p.second, p.fraction, p.error)
+    return Epoch{S}(second, fraction, error)
+end
+Base.:-(ep::Epoch, p::AstroPeriod) = ep + (-p)
 
 """
     -(a::Epoch, b::Epoch)
@@ -24,6 +27,8 @@ julia> TAIEpoch(2018, 2, 6, 20, 45, 20.0) - TAIEpoch(2018, 2, 6, 20, 45, 0.0)
 ```
 """
 function Base.:-(a::Epoch{S}, b::Epoch{S}) where S<:TimeScale
-    return ((a.second - b.second) + (a.fraction - b.fraction)) * seconds
+    second = a.second - b.second
+    fraction = (a.error - b.error) + (a.fraction - b.fraction)
+    return (fraction + second) * seconds
 end
 

src/Epochs/ranges.jl

@@ -1,6 +1,6 @@
 (::Base.Colon)(start::Epoch{S}, stop::Epoch{S}) where {S} = (:)(start, 1.0days, stop)
 
-function (::Base.Colon)(start::Epoch{S}, step::Period, stop::Epoch{S}) where S
+function (::Base.Colon)(start::Epoch{S}, step::AstroPeriod, stop::Epoch{S}) where S
     step = seconds(step)
     step = start < stop ? step : -step
     StepRangeLen(start, step, floor(Int, value(stop-start)/value(step))+1)

src/Epochs/types.jl

@@ -1,38 +1,18 @@
-@inline function two_sum(a, b)
-    hi = a + b
-    a1 = hi - b
-    b1 = hi - a1
-    lo = (a - a1) + (b - b1)
-    return hi, lo
-end
-
 struct Epoch{S<:TimeScale, T} <: Dates.AbstractDateTime
     scale::S
     second::Int64
     fraction::T
-    function Epoch{S}(second::Int64, fraction::T) where {S<:TimeScale, T}
-        return new{S, T}(S(), second, fraction)
+    error::T
+    function Epoch{S}(second::Int64, fraction::T, error::T=zero(T)) where {S<:TimeScale, T<:AbstractFloat}
+        return new{S, T}(S(), second, fraction, error)
     end
 end
 
 Epoch{S,T}(ep::Epoch{S,T}) where {S,T} = ep
 
-@inline function apply_offset(second::Int64, fraction, offset)
-    sum, residual = two_sum(fraction, offset)
-    if !isfinite(sum)
-        fraction′ = sum
-        second′ = ifelse(sum < 0, typemin(Int64), typemax(Int64))
-    else
-        int_secs = floor(Int64, sum)
-        second′ = second + int_secs
-        fraction′ = sum - int_secs + residual
-    end
-    return second′, fraction′
-end
-
 function Epoch{S}(ep::Epoch{S}, Δt) where {S<:TimeScale}
-    second, fraction = apply_offset(ep.second, ep.fraction, Δt)
-    return Epoch{S}(second, fraction)
+    second, fraction, error = apply_offset(ep.second, ep.fraction, ep.error, Δt)
+    return Epoch{S}(second, fraction, error)
 end
 
 Base.show(io::IO, ep::Epoch) = print(io, DateTime(ep), " ", timescale(ep))