Update selecteph

Project.toml

@@ -1,7 +1,7 @@
 name = "PlanetaryEphemeris"
 uuid = "d83715d0-7e5f-11e9-1a59-4137b20d8363"
 authors = ["Jorge A. Pérez Hernández", "Luis Benet", "Luis Eduardo Ramírez Montoya"]
-version = "0.8.0"
+version = "0.8.1"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"

src/interpolation.jl

@@ -186,29 +186,41 @@ function reverse(tinterp::TaylorInterpolant{T,U,N}) where {T<:Real, U<:Number, N
 end
 
 @doc raw"""
-    selecteph(eph::TaylorInterpolant{T, U, 2}, bodyind::AbstractVector{Int}; euler::Bool = false,
-              ttmtdb::Bool = false) where {T <: Real, U <: Number}
-    selecteph(eph::TaylorInterpolant{T, U, 2}, i::Int) where {T <: Real, U <: Number}
+    selecteph(eph::TaylorInterpolant, bodyind::Union{Int, AbstractVector{Int}}, t0::T = eph.t0,
+              tf::T = eph.t0 + eph.t[end]; euler::Bool = false, ttmtdb::Bool = false) where {T <: Real}
 
-Return a `TaylorInterpolant` with only the ephemeris of the bodies with indices `bodyind/i`. The keyword arguments allow to
-include lunar euler angles and/or TT-TDB.
+Return a subset of `eph` with only the ephemeris of bodies `bodyind` in timerange `[t0, tf]`.
+The keyword arguments allow to include lunar euler angles and/or TT-TDB.
 """
-function selecteph(eph::TaylorInterpolant{T, U, 2}, bodyind::AbstractVector{Int}; euler::Bool = false,
-                   ttmtdb::Bool = false) where {T, U}
+function selecteph(eph::TaylorInterpolant, bodyind::Union{Int, AbstractVector{Int}}, t0::T = eph.t0,
+                   tf::T = eph.t0 + eph.t[end]; euler::Bool = false, ttmtdb::Bool = false) where {T <: Real}
+    # Times
+    tmin, tmax = minmax(eph.t0, eph.t0 + eph.t[end])
+    @assert tmin ≤ t0 ≤ tf ≤ tmax "$tmin ≤ t0 ≤ tf ≤ $tmax"
+    if issorted(eph.t)
+        i_0 = searchsortedlast(eph.t, t0)
+        i_f = searchsortedfirst(eph.t, tf)
+    else
+        i_0 = searchsortedlast(eph.t, tf, rev = true)
+        i_f = searchsortedfirst(eph.t, t0, rev = true)
+    end
+    # Degrees of freedom
     N = numberofbodies(eph)
+    @assert all(bodyind .< N) "bodyind .< $N" 
     idxs = nbodyind(N, bodyind)
     if euler
-        idxs = union(idxs, 6N+1:6N+12)
+        idxs = vcat(idxs, 6N+1:6N+12)
     end
     if ttmtdb
-        idxs = union(idxs, 6N+13)
+        idxs = vcat(idxs, 6N+13)
     end
-    x = view(eph.x, :, idxs)
-    return TaylorInterpolant(eph.t0, eph.t, x)
+    # Views
+    t = view(eph.t, i_0:i_f)
+    x = view(eph.x, i_0:i_f-1, idxs)
+    # New TaylorInterpolant
+    return TaylorInterpolant(eph.t0, t, x)
 end
 
-selecteph(eph::TaylorInterpolant, i::Int) = selecteph(eph, i:i)
-
 function join(bwd::TaylorInterpolant{T, U, 2}, fwd::TaylorInterpolant{T, U, 2}) where {T, U}
     @assert bwd.t0 == fwd.t0 "Initial time must be the same for both TaylorInterpolant"
     order_bwd = get_order(bwd.x[1, 1])

test/propagation.jl

@@ -98,6 +98,24 @@ using LinearAlgebra: norm
 
         @test selecteph(sol64, bodyind, euler = true, ttmtdb = true) == recovered_sol64
 
+        # Test selecteph
+        t0 = sol64.t0 + sol64.t[end]/3
+        tf = sol64.t0 + 2*sol64.t[end]/3
+        idxs = vcat(nbodyind(N, [su, ea, mo]), 6N+1:6N+13)
+        i_0 = searchsortedlast(sol64.t, t0)
+        i_f = searchsortedfirst(sol64.t, tf)
+
+        subsol = selecteph(sol64, [su, ea, mo], t0, tf; euler = true, ttmtdb = true)
+
+        @test subsol.t0 == sol64.t0
+        @test subsol.t0 + subsol.t[1] ≤ t0
+        @test subsol.t0 + subsol.t[end] ≥ tf
+        @test size(subsol.x) == (i_f - i_0, length(idxs))
+        @test subsol.x == sol64.x[i_0:i_f-1, idxs]
+        @test subsol(t0) == sol64(t0)[idxs]
+        @test subsol(tf) == sol64(tf)[idxs]
+
+        # Kernels URLs
         LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls"
         TTmTDBK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/TTmTDB.de430.19feb2015.bsp"
         SPK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de430_1850-2150.bsp"