Reworked EOP parsing and loading

src/IERSConventions.jl

@@ -1,14 +1,16 @@
 module IERSConventions
 
-    using Tempo
-
+    using DelimitedFiles
     using PrecompileTools
     using ReferenceFrameRotations
     using StaticArrays
 
+    using Tempo
+    using JSMDInterfaces.Math: AbstractInterpolationMethod, interpolate
+    using JSMDUtils.Math: InterpAkima, arcsec2rad
+
     # Basic definitions
     include("models.jl")
-    include("angles.jl")
     include("poisson.jl")
 
     # Fundamental arguments
@@ -32,7 +34,10 @@ module IERSConventions
     include("polar.jl")
 
     # EOP functions 
-    include("eop.jl")
+    include("eop/data.jl")
+    include("eop/parsers.jl")
+    include("eop/loaders.jl")
+    include("eop/retrieval.jl")
 
     # Rotation functions 
     include("rotations.jl")

src/eop/data.jl

@@ -0,0 +1,204 @@
+
+export eop_filename
+
+# EOP Data 
+# ============================
+
+"""
+    NutationCorrections
+
+Container to hold the nutation corrections and celestial pole offsets associated to a  
+given IERS model (e.g., 2000A)
+
+### Fields    
+- `δX, δY`: Celestial pole offsets referred to the model IAU2000A (rad).
+- `δΔψ, δΔϵ`: Nutation corrections in longitude and obliquity (rad)
+
+### See also 
+See also [`EOPData`](@ref).
+"""
+struct NutationCorrections{T}
+    δX::Vector{T}
+    δY::Vector{T}
+    δΔψ::Vector{T}
+    δΔϵ::Vector{T}
+end
+
+function NutationCorrections(::Type{T}=Float64) where T 
+    return NutationCorrections(T[], T[], T[], T[])
+end
+
+
+"""
+    EOPData{T}
+
+EOP Data container. Data is parameterised by time expressed in Terrestrial Time (TT) days. 
+
+### Fields 
+- `filename` : File where the EOP data are stored. 
+- `days_UTC`: UTC days since J2000.
+- `days_TT` : TT days since J2000.
+- `UT1_TT`: UT1 minus TT offset, in seconds.
+- `LOD`: Length of day offset (s).
+- `xp, yp`: Polar motion with respect to the crust (rad).
+- `nut1996, nut2003, nut2006`: nutation and CIP corrections.  
+"""
+mutable struct EOPData{T}
+
+    filename::String 
+
+    days_UTC::Vector{T}
+    days_TT::Vector{T}
+
+    UT1_TT::Vector{T}
+    LOD::Vector{T}
+
+    xp::Vector{T}
+    yp::Vector{T}
+
+    nut1996::NutationCorrections{T}
+    nut2003::NutationCorrections{T}
+    nut2010::NutationCorrections{T}
+
+end
+
+function EOPData(::Type{T}=Float64) where T 
+    return EOPData(
+        "", T[], T[], T[], T[], T[], T[], 
+        NutationCorrections(T), NutationCorrections(T), NutationCorrections(T)
+    )
+end 
+
+function Base.show(io::IO, eop::EOPData)
+
+    if isempty(eop.days_UTC)
+        println(io, "EOPData()")
+    else 
+        println(
+            io, 
+            "EOPData(filename=\"$(eop.filename), \"beg=\"$(eop.days_UTC[1]) (UTC)\", "* 
+            "end=\"$(eop.days_UTC[end]) (UTC)\")"
+        )
+    end 
+
+end
+
+"""
+    set_nutation_corr!(eop::EOPData, m::IERSModel, nc::NutationCorrections)
+
+Set the nutation and cip corrections associated to the IERS model `m` to `nc`.
+"""
+set_nutation_corr!(eop::EOPData, ::IERS1996, nc::NutationCorrections) = eop.nut1996 = nc
+set_nutation_corr!(eop::EOPData, ::IERS2003, nc::NutationCorrections) = eop.nut2003 = nc
+set_nutation_corr!(eop::EOPData, ::IERS2010, nc::NutationCorrections) = eop.nut2010 = nc
+
+"""
+    eop_filename(eop::EOPData) 
+
+Get the loaded Earth Orientation Parameters (EOP) filename.
+"""
+function eop_filename(eop::EOPData)
+
+    if isempty(eop.filename) 
+        throw(ErrorException("Unable to retrieve filename, no EOP data has been loaded."))
+    end
+
+    return eop.filename
+end
+
+
+# Interpolators 
+# ============================
+
+"""
+    NutCorrectionsInterpolator
+
+Container to store the interpolators for the nutation corrections and celestial pole offsets 
+associated to a given IERS model. 
+
+### See also 
+See also [`EOPInterpolator`](@ref)
+"""
+struct NutCorrectionsInterpolator{T <: AbstractInterpolationMethod}
+    δX::T
+    δY::T
+    δΔψ::T
+    δΔϵ::T
+end
+
+function NutCorrectionsInterpolator(::Type{T}) where T
+    NutCorrectionsInterpolator(
+        _akima_init(T), _akima_init(T), _akima_init(T), _akima_init(T)
+    )
+end 
+
+function NutCorrectionsInterpolator(days_tt, nc::NutationCorrections) 
+    NutCorrectionsInterpolator(
+        InterpAkima(days_tt, nc.δX), 
+        InterpAkima(days_tt, nc.δY), 
+        InterpAkima(days_tt, nc.δΔψ), 
+        InterpAkima(days_tt, nc.δΔϵ)
+    )
+end
+
+"""
+    EOPInterpolator
+
+Container to store the interpolators for the loaded EOP data. 
+"""
+mutable struct EOPInterpolator{T <: AbstractInterpolationMethod}
+
+    init::Bool 
+
+    xp::T 
+    yp::T 
+
+    ut1_tt::T 
+    lod::T
+
+    nut1996::NutCorrectionsInterpolator{T}
+    nut2003::NutCorrectionsInterpolator{T}
+    nut2010::NutCorrectionsInterpolator{T}
+
+end
+
+function EOPInterpolator(::Type{T}=Float64) where T 
+    EOPInterpolator(
+        false, 
+        _akima_init(T), _akima_init(T), _akima_init(T), _akima_init(T),
+        NutCorrectionsInterpolator(T), 
+        NutCorrectionsInterpolator(T), 
+        NutCorrectionsInterpolator(T)
+    )
+end
+
+function Base.show(io::IO, i::EOPInterpolator)
+    println(io, "EOPInterpolator(init=$(i.init))")
+end
+
+# Function to initialise a dummy Akima interpolator 
+_akima_init(::Type{T}) where T = InterpAkima([1, 2, 3, 4, 5, 6], zeros(T, 6))
+
+
+# IERS Constants  
+# ============================
+
+"""
+    IERS_EOP_DATA 
+
+Earth Orientation Parameters Data. 
+
+### See also 
+See also [`EOPData`](@ref).
+"""
+const IERS_EOP_DATA = EOPData();
+
+"""
+    IERS_EOP 
+
+Earth Orientation Parameters (EOP) interpolators. 
+
+### See also 
+See also [`load_eop!`](@ref) and [`EOPInterpolator`](@ref).
+"""
+const IERS_EOP = EOPInterpolator();

src/eop/loaders.jl

@@ -0,0 +1,129 @@
+
+export eop_load_data!
+
+"""
+    eop_load_data!(filename, m::IERS2010)
+
+Initialise the Earth Orientation Parameters (EOP) from a dedicated JSMD `.eop.dat` file.
+
+!!! note 
+    Currently, only EOP data associated to the IAU2006/2000A precession-nutation model 
+    is supported.
+"""
+function eop_load_data!(filename, m::IERS2010)
+
+    # Set the EOP data 
+    eop_set_data!(filename, m)
+
+    # Initialise the interpolators 
+    IERS_EOP.init = true 
+
+    # Set polar motion
+    IERS_EOP.xp = InterpAkima(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.xp)
+    IERS_EOP.yp = InterpAkima(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.yp)
+
+    IERS_EOP.lod = InterpAkima(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.LOD)
+    IERS_EOP.ut1_tt = InterpAkima(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.UT1_TT)
+
+    # Set nutation correction interpolators
+    # TODO: ensure these are set! 
+    # nci1996 = NutCorrectionsInterpolator(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.nut1996)
+    # nci2003 = NutCorrectionsInterpolator(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.nut2003)
+    nci2010 = NutCorrectionsInterpolator(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.nut2010)
+
+    # IERS_EOP.nut1996 = nci1996
+    # IERS_EOP.nut2003 = nci2003 
+    IERS_EOP.nut2010 = nci2010
+
+    @info "EOP initialised from file `$(filename)`."
+    nothing 
+
+end
+
+
+"""
+    eop_set_data!(filename, m::IERSModel)  
+
+Set Earth Orientation Parameters (EOP) to be used for frames transformations from a JSMD 
+`.eop.dat` file, given a reference model `m`.
+"""
+function eop_set_data!(filename, m::IERSModel)
+
+    oldfile = IERS_EOP_DATA.filename
+    days_utc, days_tt, xp, yp, ut1_tt, lod, δX, δY, δΔψ, δΔϵ = eop_read_data(filename)
+
+    if (!isempty(IERS_EOP_DATA.days_UTC))
+        @warn "Existing EOP data from \'$(oldfile)\' will be overwritten by \'$(filename)\'."
+    end
+
+    # Set data time stamps
+    IERS_EOP_DATA.filename = filename
+    IERS_EOP_DATA.days_UTC = days_utc
+    IERS_EOP_DATA.days_TT  = days_tt 
+
+    # Set TT-to-UT1 offset and LOD
+    IERS_EOP_DATA.UT1_TT = ut1_tt 
+    IERS_EOP_DATA.LOD    = lod
+
+    # Set polar motion coordinates 
+    IERS_EOP_DATA.xp = xp 
+    IERS_EOP_DATA.yp = yp
+
+    # Set the nutation corrections 
+    nc = NutationCorrections(δX, δY, δΔψ, δΔϵ)
+    set_nutation_corr!(IERS_EOP_DATA, m, nc) 
+
+    # TODO: convert these corrections to those for the other models!
+    nothing
+
+end
+
+
+"""
+    eop_read_data(filename)  
+
+Read Earth Orientation Parameters (EOP) from a dedicated JSMD `.eop.dat` file. 
+
+!!! note 
+    JSMD's `.eop.dat` EOP files are meant to have a fixed structure to ease the 
+    retrieval of the relevant EOP data. 
+"""
+function eop_read_data(filename)
+
+    if !endswith(filename, "eop.dat")
+        throw(
+            ArgumentError(
+                "EOP reader support only '.eop.dat' files! Please prepare " * 
+                "the data with \'eop_parse_csv\' and retry."
+            )
+        )
+    end
+
+    # Load the file
+    data = readdlm(filename; header=false)
+
+    # Retrieve UT1-UTC 
+    Δut1 = @view(data[:, 4])
+
+    # Convert UTC days to TT seconds
+    days_utc = @view(data[:, 1])
+    days_tai = map(t->Tempo.utc2tai(Tempo.DJ2000, t)[2], days_utc)
+    days_tt  = days_tai .+ Tempo.OFFSET_TAI_TT/Tempo.DAY2SEC
+
+    # Compute UT1-TT, in seconds
+    sec_ut1 = days_utc/Tempo.DAY2SEC + Δut1
+    ut1_tt = sec_ut1 - days_tt*Tempo.DAY2SEC
+
+    # Retrieve other quantities
+    xp = @view(data[:, 2])
+    yp = @view(data[:, 3])
+    lod = @view(data[:, 5])
+
+    δX = data[:, 6]
+    δY = data[:, 7]
+    δΔψ = data[:, 6]
+    δΔϵ = data[:, 7]
+
+    return days_utc, days_tt, xp, yp, ut1_tt, lod, δX, δY, δΔψ, δΔϵ
+
+end