Merge pull request #46 from JuliaSpaceMissionDesign/dev

Update to v1.2.1

Author: MicheleCeresoli

.github/workflows/ci.yml

@@ -3,7 +3,7 @@ on:
   - push
   - pull_request
 env:
-  JULIA_NUM_THREADS: 1
+  JULIA_NUM_THREADS: 4
 jobs:
   test:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}

Artifacts.toml

@@ -1,7 +1,7 @@
 [testdata]
-git-tree-sha1 = "b58880275ed655c31bcb6be63df686bc1539af6e"
-lazy = true
+git-tree-sha1 = "d97cce4a3c744330b363c292097a96a1aab69e6d"
+lazy = true 
 
     [[testdata.download]]
-    sha256 = "4fc7a9ef011b2f90379f5129e3015d30ff313dd0493724cec5c55fc560180e13"
-    url = "https://gist.github.com/MicheleCeresoli/ff9f8eeea149022b32979e94ef969d96/raw/b58880275ed655c31bcb6be63df686bc1539af6e.tar.gz"
+    sha256 = "4d8ca82289ffe9161cdb49b68b0a373b472fc393255544d5d43ca6ead1774baa"
+    url = "https://gist.github.com/MicheleCeresoli/141622fbf94e6da4c9fec7ff815192da/raw/d97cce4a3c744330b363c292097a96a1aab69e6d.tar.gz"

Project.toml

@@ -1,22 +1,22 @@
 name = "Ephemerides"
 uuid = "6a9c3322-c8fe-4c26-8ad6-14a6f8acd2a0"
 authors = ["JSMD Team"]
-version = "1.2.0"
+version = "1.2.1"
 
 [deps]
 JSMDInterfaces = "6b30ee2f-618e-4a15-bf4e-7df7b496e609"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
 [compat]
 JSMDInterfaces = "1.4"
 LazyArtifacts = "1"
 PreallocationTools = "0.4"
 PrecompileTools = "1"
-StaticArrays = "1"
+StaticArraysCore = "1.4"
 julia = "1.7"
 
 [extras]

README.md

@@ -67,9 +67,9 @@ angles = ephem_rotation3(eph, 1, 31006, time)
 ```
 
 ## Current Limitations
-- The supported JPL binary SPK/PCK segments types are: 1, 2, 3, 8, 9, 12, 13, 14, 18, 19, 20, 21. 
+- The supported JPL binary SPK/PCK segments types are: 1, 2, 3, 5, 8, 9, 12, 13, 14, 15, 17, 18, 19, 20, 21. 
 - Binary INPOP kernels are not supported. 
-- Acceleration and jerk computations are natively unavailable for SPK segments of type 1 and 21, but can be computed with Automatic Differentiation (AD) of the position and velocity components.
+- Acceleration and jerk computations are natively unavailable for SPK segments of type 1, 5, 15, 17 and 21, but can be computed with Automatic Differentiation (AD) of the position and velocity components.
 
 ## Documentation 
 For further information on this package please refer to the [stable documentation](https://juliaspacemissiondesign.github.io/Ephemerides.jl/stable/)

docs/src/api/lapi.md

@@ -125,6 +125,13 @@ Ephemerides.SPKSegmentCache2
 Ephemerides.SPKSegmentType2
 ```
 
+### SPK Type 5
+```@docs 
+Ephemerides.SPKSegmentHeader5
+Ephemerides.SPKSegmentCache5
+Ephemerides.SPKSegmentType5
+```
+
 ### SPK Type 8 and 12
 ```@docs 
 Ephemerides.SPKSegmentHeader8
@@ -148,6 +155,27 @@ Ephemerides.SPKSegmentType14
 !!! note 
     SPK segments of type 14 have the same cache structure of SPK type 2 and 3.
 
+
+### SPK Type 15
+```@docs
+Ephemerides.SPKSegmentHeader15
+Ephemerides.SPKSegmentType15
+```
+
+!!! note 
+    The cache of SPK Type 15 segments is made of [`Ephemerides.TwoBodyUniversalCache`](@ref) objects.
+
+
+### SPK Type 17
+```@docs
+Ephemerides.SPKSegmentHeader17
+Ephemerides.SPKSegmentType17
+```
+
+!!! note 
+    SPK segments of type 17 do not require a cache structure.
+
+
 ### SPK Type 18 and 19
 ```@docs
 Ephemerides.SPKSegmentHeader18
@@ -181,8 +209,6 @@ Ephemerides.find_logical_record
 Ephemerides.get_coefficients!
 ```
 
-### Parsing 
-
 ## Interpolating Functions
 
 ### Caches 
@@ -227,5 +253,12 @@ Ephemerides.final_times
 Ephemerides.analyse_timespan 
 ```
 
+## TwoBody Routines 
+```@docs 
+Ephemerides.TwoBodyUniversalCache 
+Ephemerides.update_cache!
+Ephemerides.propagate_twobody
+Ephemerides.stumpff
+```
 
 

docs/src/index.md

@@ -62,6 +62,6 @@ angles = ephem_rotation3(eph, 1, 31006, time)
 ```
 
 ## Current Limitations
-- The supported JPL binary SPK/PCK segments types are: 1, 2, 3, 8, 9, 12, 13, 14, 18, 19, 20, 21. 
+- The supported JPL binary SPK/PCK segments types are: 1, 2, 3, 5, 8, 9, 12, 13, 14, 15, 17, 18, 19, 20, 21. 
 - Binary INPOP kernels are not supported. 
-- Acceleration and jerk computations are unavailable for SPK segments of type 1 and 21.
+- Acceleration and jerk computations are unavailable for SPK segments of type 1, 5, 15, 17 and 21.

docs/src/tutorials/load.md

@@ -2,7 +2,7 @@
 
 This tutorial will walk you through the basic features and interfaces that allow you to load binary ephemeris kernels.
 
-The supported sources of ephemerides are currently limited to binary PCK and SPK segments of type: 1, 2, 3, 8, 9, 12, 13, 18, 19, 20 and 21. 
+The supported sources of ephemerides are currently limited to binary PCK and SPK segments of type: 1, 2, 3, 5, 8, 9, 12, 13, 14, 15, 17, 18, 19, 20, and 21. 
 
 !!! note 
     Support for IMCCE INPOP ephemerides is yet to be implemented.

docs/src/tutorials/position.md

@@ -50,7 +50,7 @@ pvaj = ephem_vector12(eph, 3, 399, 0)
 In all these examples, the returned data is always in the form of a `StaticArray` in order to minimise memory allocations.
 
 !!! warning 
-    SPK segments of types 1 and 21 do not natively support acceleration and jerk computations. However, these values can be computed by Automatic Differentiation (AD) of the position and/or velocity components.
+    SPK segments of types 1 5, 15, 17 and 21 do not natively support acceleration and jerk computations. However, these values can be computed by Automatic Differentiation (AD) of the position and/or velocity components.
 
 
 ## Computing orientation angles

src/Ephemerides.jl

@@ -5,7 +5,7 @@ using LazyArtifacts
 using Mmap
 using PreallocationTools
 using PrecompileTools: PrecompileTools
-using StaticArrays
+using StaticArraysCore
 
 import JSMDInterfaces.Ephemeris as jEph
 
@@ -18,6 +18,9 @@ include("interp/hermite.jl")
 include("interp/lagrange.jl")
 include("interp/chebyshev.jl")
 
+# Twobody utilities 
+include("twobody.jl")
+
 # SPK segment types definitions
 include("spk/spktypes.jl")
 
@@ -27,9 +30,12 @@ include("daf.jl")
 # SPK segments implementations
 include("spk/spk1.jl")
 include("spk/spk2.jl")
+include("spk/spk5.jl")
 include("spk/spk8.jl")
 include("spk/spk9.jl")
 include("spk/spk14.jl")
+include("spk/spk15.jl")
+include("spk/spk17.jl")
 include("spk/spk18.jl")
 include("spk/spk19.jl")
 include("spk/spk20.jl")

src/daf.jl

@@ -588,25 +588,34 @@ function create_spk_segment(daf::DAF, desc::DAFSegmentDescriptor)
 
     mapped_spktype = SPK_SEGMENTLIST_MAPPING[segment_type(desc)]
     if mapped_spktype == 1
-        SPKSegmentType1(daf, desc)
+        SPKSegmentType2(daf, desc)
 
     elseif mapped_spktype == 2
-        SPKSegmentType2(daf, desc)
+        SPKSegmentType9(daf, desc)
 
     elseif mapped_spktype == 3 
-        SPKSegmentType8(daf, desc)
+        SPKSegmentType1(daf, desc)
 
     elseif mapped_spktype == 4 
-        SPKSegmentType9(daf, desc)
+        SPKSegmentType14(daf, desc)
 
     elseif mapped_spktype == 5
-        SPKSegmentType14(daf, desc)
+        SPKSegmentType15(daf, desc)
 
     elseif mapped_spktype == 6 
+        SPKSegmentType8(daf, desc)
+
+    elseif mapped_spktype == 7
         SPKSegmentType19(daf, desc)
 
-    else
+    elseif mapped_spktype == 8
         SPKSegmentType20(daf, desc)
+        
+    elseif mapped_spktype == 9
+        SPKSegmentType5(daf, desc)
+        
+    else 
+        SPKSegmentType17(daf, desc)
     end
 
 end

src/precompile.jl

@@ -4,18 +4,21 @@ PrecompileTools.@setup_workload begin
     test_dir = artifact"testdata"
 
     spks = [
-        "example1spk_seg1.bsp",
-        "example1spk_seg2.bsp",
-        "example1spk_seg3.bsp",
-        "example1spk_seg8.bsp",
-        "example1spk_seg9.bsp",
-        "example1spk_seg12.bsp",
-        "example1spk_seg13.bsp",
-        "example1spk_seg14.bsp",
-        "example1spk_seg18.bsp",
-        "example1spk_seg19.bsp",
-        "example1spk_seg20.bsp",
-        "example1spk_seg21.bsp",
+        "spk1_ex1.bsp",
+        "spk2_ex1.bsp",
+        "spk3_ex1.bsp",
+        "spk5_ex1.bsp",
+        "spk8_ex1.bsp",
+        "spk9_ex1.bsp",
+        "spk12_ex1.bsp",
+        "spk13_ex1.bsp",
+        "spk14_ex1.bsp",
+        "spk15_ex1.bsp",
+        "spk17_ex1.bsp",
+        "spk18_ex1.bsp",
+        "spk19_ex1.bsp",
+        "spk20_ex1.bsp",
+        "spk21_ex1.bsp",
     ]
 
     pcks = [
@@ -37,7 +40,7 @@ PrecompileTools.@setup_workload begin
             _ = ephem_vector3(ephj, cid, tid, tj);
             _ = ephem_vector6(ephj, cid, tid, tj);
 
-            if !(segment_type(desc) in (1, 21))
+            if !(segment_type(desc) in (1, 5, 15, 17, 21))
                 _ = ephem_vector9(ephj, cid, tid, tj);
                 _ = ephem_vector12(ephj, cid, tid, tj);
             end

src/spk/spk1.jl

@@ -58,8 +58,8 @@ Initialise the cache for an SPK segment of type 1.
 """
 function SPKSegmentCache1(head::SPKSegmentHeader1)
     SPKSegmentCache1(
-        Int[0.0], zeros(head.maxdim), zeros(3), zeros(3), zeros(head.maxdim, 3), 
-        Int[0.0], zeros(Int, 3), MVector(-1), 
+        0.0, zeros(head.maxdim), zeros(3), zeros(3), zeros(head.maxdim, 3), 
+        0, zeros(Int, 3), -1, 
         DiffCache(zeros(head.maxdim-1)), DiffCache(zeros(head.maxdim-2)), 
         DiffCache(zeros(head.maxdim+2)), DiffCache(zeros(3))
     )
@@ -89,7 +89,7 @@ function spk_vector3(daf::DAF, seg::SPKSegmentType1, time::Number)
     # Retrieve the MDA coefficients 
     get_coefficients!(daf, header(seg), cache(seg), index)
 
-    Δ = time - cache(seg).tl[1]
+    Δ = time - cache(seg).tl
 
     # Compute MDA position coefficients
     compute_mda_pos_coefficients!(cache(seg), Δ)
@@ -105,7 +105,7 @@ function spk_vector6(daf::DAF, seg::SPKSegmentType1, time::Number)
     # Retrieve the MDA coefficients 
     get_coefficients!(daf, header(seg), cache(seg), index)
 
-    Δ = time - cache(seg).tl[1]
+    Δ = time - cache(seg).tl
 
     # Compute MDA position coefficients
     compute_mda_pos_coefficients!(cache(seg), Δ)
@@ -115,10 +115,10 @@ function spk_vector6(daf::DAF, seg::SPKSegmentType1, time::Number)
     compute_mda_vel_coefficients!(cache(seg), Δ)
     vel = compute_mda_velocity(cache(seg), Δ)
 
-    return @inbounds SA[
+    return @inbounds SVector{6}(
         pos[1], pos[2], pos[3], 
         vel[1], vel[2], vel[3]
-    ]
+    )
 end
 
 function spk_vector9(::DAF, ::SPKSegmentType1, ::Number)
@@ -179,12 +179,12 @@ function get_coefficients!(daf::DAF, head::SPKSegmentHeader1, cache::SPKSegmentC
             index::Int)
 
     # Check whether the coefficients for this record are already loaded
-    index == cache.id[1] && return nothing
-    cache.id[1] = index 
+    index == cache.id && return nothing
+    cache.id = index 
 
     # Initial record address 
     i0 = 8*(head.iaa - 1) + 8*head.recsize*index
-    cache.tl[1] = get_float(array(daf), i0, endian(daf))
+    cache.tl = get_float(array(daf), i0, endian(daf))
 
     @inbounds for k = 1:head.maxdim 
         cache.g[k] = get_float(array(daf), i0 + 8k, endian(daf))    
@@ -207,7 +207,7 @@ function get_coefficients!(daf::DAF, head::SPKSegmentHeader1, cache::SPKSegmentC
     end
 
     i4 = i3 + 3*head.maxdim*8
-    cache.kqmax[1] = Int(get_float(array(daf), i4, endian(daf)))
+    cache.kqmax = Int(get_float(array(daf), i4, endian(daf)))
 
     @inbounds for j = 1:3 
         cache.kq[j] = get_float(array(daf), i4 + 8j, endian(daf))
@@ -228,11 +228,11 @@ function compute_mda_position(cache::SPKSegmentCache1, Δ::Number)
     end
 
     # At this point by definition of the above while loop, ks will always be = 1
-    return SA[
+    return SVector{3}(
         cache.refpos[1] + Δ*(cache.refvel[1] + Δ*get_tmp(cache.vct, Δ)[1]),
         cache.refpos[2] + Δ*(cache.refvel[2] + Δ*get_tmp(cache.vct, Δ)[2]),
         cache.refpos[3] + Δ*(cache.refvel[3] + Δ*get_tmp(cache.vct, Δ)[3])
-    ]
+    )
 
 end
 
@@ -249,11 +249,11 @@ function compute_mda_velocity(cache::SPKSegmentCache1, Δ::Number)
     end
 
     # At this point by definition of the above while loop, ks will always be = 0
-    return SA[
+    return SVector{3}(
         cache.refvel[1] + Δ*get_tmp(cache.vct, Δ)[1],
         cache.refvel[2] + Δ*get_tmp(cache.vct, Δ)[2],
         cache.refvel[3] + Δ*get_tmp(cache.vct, Δ)[3]
-    ]
+    )
 
 end
 
@@ -262,8 +262,8 @@ end
 @inbounds function compute_mda_pos_coefficients!(cache::SPKSegmentCache1, Δ::Number)
 
     tp = Δ
-    mq2 = cache.kqmax[1] - 2
-    ks = cache.kqmax[1] - 1
+    mq2 = cache.kqmax - 2
+    ks = cache.kqmax - 1
 
     get_tmp(cache.fc, Δ)[1] = 1.0
     for j = 1:mq2
@@ -273,7 +273,7 @@ end
     end
 
     # compute inverse coefficients 
-    for j = 1:cache.kqmax[1]
+    for j = 1:cache.kqmax
         get_tmp(cache.w, Δ)[j] = 1.0 / j
     end
 
@@ -297,7 +297,7 @@ end
 # This function is such that the output ks = 0
 @inbounds function compute_mda_vel_coefficients!(cache::SPKSegmentCache1, Δ::Number)
 
-    for j = 1:cache.kqmax[1]
+    for j = 1:cache.kqmax
         get_tmp(cache.w, Δ)[j+1] = get_tmp(cache.fc, Δ)[j+1]*get_tmp(cache.w, Δ)[j] - get_tmp(cache.wc, Δ)[j]*get_tmp(cache.w, Δ)[j+1]
     end