Implement SPICE (bsp) for planet ephemeris

Author: akoumjian

.github/workflows/c.yml

@@ -3,7 +3,7 @@ name: ASSIST Unit Tests (C)
 on: [push, pull_request]
 
 env:
-  JPL_PLANET_EPHEM: ../../data/linux_p1550p2650.440
+  JPL_PLANET_EPHEM: ../../data/de440.bsp
 
 jobs:
   build:
@@ -17,19 +17,11 @@ jobs:
         run: | 
           pwd
           git clone https://github.com/hannorein/rebound.git
-          #- name: Compile REBOUND
-          #  working-directory: ../rebound
-          #  run: | 
-          #    pwd
-          #    make
-          #- name: Compile ASSIST
-          #  run: | 
-          #    make
       - name: Download DE441 dataset
         working-directory: ./data/
         run: |
-          #wget --no-verbose https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de441/linux_m13000p17000.441  
-          wget --no-verbose https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de440/linux_p1550p2650.440
+          #wget --no-verbose https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de441.bsp
+          wget --no-verbose https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de440.bsp
           wget --no-verbose https://ssd.jpl.nasa.gov/ftp/eph/small_bodies/asteroids_de441/sb441-n16.bsp
       - name: Reproduce Holman trajectory
         working-directory: ./unit_tests/holman/
@@ -91,3 +83,18 @@ jobs:
         run: |
           make
           LD_LIBRARY_PATH=../../../rebound/src/ ./rebound
+      - name: Test Loading Spice Kernels
+        working-directory: ./unit_tests/spk_init/
+        run: |
+          make
+          LD_LIBRARY_PATH=../../../rebound/src/ ./rebound
+      - name: Test Loading Constants from Planet SPK
+        working-directory: ./unit_tests/spk_load_constants/
+        run: |
+          make
+          LD_LIBRARY_PATH=../../../rebound/src/ ./rebound
+      - name: Test Joining Mass Data from SPK
+        working-directory: ./unit_tests/spk_join_masses/
+        run: |
+          make
+          LD_LIBRARY_PATH=../../../rebound/src/ ./rebound

.github/workflows/python.yml

@@ -24,8 +24,8 @@ jobs:
       - name: Download DE441 dataset
         working-directory: ./data/
         run: |
-          #wget --no-verbose https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de441/linux_m13000p17000.441  
-          wget --no-verbose https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de440/linux_p1550p2650.440
+          #wget --no-verbose https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de441.bsp  
+          wget --no-verbose https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de440.bsp
           wget --no-verbose https://ssd.jpl.nasa.gov/ftp/eph/small_bodies/asteroids_de441/sb441-n16.bsp
       - name: Build package
         run: pip install -e .

.python-version

@@ -0,0 +1 @@
+assist

MANIFEST.in

@@ -1,10 +1,8 @@
 include src/assist.c
 include src/forces.c
-include src/planets.c
 include src/spk.c
 include src/assist.h
 include src/forces.h
-include src/planets.h
 include src/spk.h
 include README.md
 include LICENSE

Makefile

@@ -18,3 +18,25 @@ doc:
 	cd doc/doxygen && doxygen
 	$(MAKE) -C doc html
 
+
+# Iterate through each directory in the unit_tests directory and compile the test files
+# then run them
+test:
+	@echo "Running tests ..."
+	@set -e; \
+	for dir in $(wildcard unit_tests/*); do \
+		if [ -d $$dir ]; then \
+			echo "Entering directory $$dir"; \
+			if $(MAKE) -C $$dir; then \
+				if [ -f "$$dir/rebound" ]; then \
+					(cd $$dir && LD_LIBRARY_PATH=$(REB_DIR):. ./rebound); \
+				else \
+					echo "No rebound executable found in $$dir"; \
+					exit 1; \
+				fi \
+			else \
+				echo "Make failed in $$dir"; \
+				exit 1; \
+			fi \
+		fi \
+	done

README.md

@@ -25,18 +25,36 @@ Now we can install numpy, REBOUND, and ASSIST:
     pip install rebound 
     pip install assist
 
-To use use ASSIST, you also need to download ephemeris data files. One file for planet ephemeris and another suplementary file for asteroid ephemeris. The following commands download these files with curl. You can also manually download them using your browser. Note that these are large files, almost 1GB in size.
+To use use ASSIST, you also need to download ephemeris data files. One file for planet ephemeris and another suplementary file for asteroid ephemeris. You can do this with Python packages or by downloading files directly.  Note that these are large files, almost 1GB in size.
+
+    pip install naif-de440
+    pip install jpl-small-bodies-de441-n16
+
+
+You can initialize the ephemeris data from the packages like so:
+
+    python3
+
+    >>> import assist
+    >>> from naif_de440 import de440
+    >>> from jpl_small_bodies_de441_n16 import de441_n16
+    >>> ephem = assist.Ephem(de440, de441_n16)
+
+The variables to the ephemeris files are simply path strings to the files. Alternatively, you can download these files with curl or your browser.
 
     mkdir data
-    curl https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de440/linux_p1550p2650.440 -o data/linux_p1550p2650.440
+    curl https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de440.bsp -o data/de440.bsp
     curl https://ssd.jpl.nasa.gov/ftp/eph/small_bodies/asteroids_de441/sb441-n16.bsp -o data/sb441-n16.bsp
 
-Now you can try out if assist works.
+Now you can point assist to those files directly, like so:
 
     python3
 
     >>> import assist
-    >>> ephem = assist.Ephem("data/linux_p1550p2650.440", "data/sb441-n16.bsp")
+    >>> ephem = assist.Ephem("data/de440.bsp", "data/sb441-n16.bsp")
+
+Once you've initialized the ephemeris data, you can test that assist is working by running the following commands:
+
     >>> print(ephem.jd_ref)
     >>> ephem.get_particle("Earth", 0)
 
@@ -51,12 +69,12 @@ To install the C version of ASSIST, first clone the REBOUND and then the ASSIST
 
 To use use ASSIST, you also need to download ephemeris data files. One file for planet ephemeris and another suplementary file for asteroid ephemeris. The following commands download these files with curl. You can also manually download them using your browser. Note that these are large files, almost 1GB in size.
 
-    curl https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de440/linux_p1550p2650.440 -o assist/data/linux_p1550p2650.440
+    curl https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de440.bsp -o assist/data/de440.bsp
     curl https://ssd.jpl.nasa.gov/ftp/eph/small_bodies/asteroids_de441/sb441-n16.bsp -o assist/data/sb441-n16.bsp
 
 For some of the examples, you will also need the planet ephemeris file with an extended coverage. Note that this file is 2.6GB in size.
 
-    curl https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de441/linux_m13000p17000.441 -o assist/data/linux_m13000p17000.441
+    curl https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de441.bsp -o assist/data/de441.bsp
 
 Next, go to one of the example directories and compile the problem file. This will also trigger the installation of the REBOUND and ASSIST shared libraries.
 
@@ -81,5 +99,6 @@ ASSIST is open source, freely distributed under the [GNU General Public license,
 * Robert Weryk, University of Western Ontario
 * Dan Tamayo, Harvey Mudd College, <[email protected]>
 * David M. Hernandez, Center for Astrophysics | Harvard & Smithsonian 
+* Alec Koumjian, Asteroid Institute | B612 Foundation, <[email protected]>
 
 

assist/ephem.py

@@ -1,10 +1,14 @@
+import warnings
+from ctypes import (CFUNCTYPE, POINTER, Structure, byref, c_char, c_char_p,
+                    c_double, c_int, c_long, c_uint, c_uint32, c_ulong,
+                    c_void_p, cast)
 from typing import Optional, Union
 
-from . import clibassist, assist_error_messages
-from ctypes import Structure, c_double, POINTER, c_int, c_uint, c_long, c_ulong, c_void_p, c_char_p, CFUNCTYPE, byref, c_uint32, c_uint, cast, c_char
 import rebound
+
 import assist
-import warnings
+
+from . import assist_error_messages, clibassist
 
 ASSIST_BODY_IDS = {
         0: "Sun",
@@ -82,8 +86,10 @@ def get_particle(self, body: Union[int, str], t: float) -> rebound.Particle:
     def __del__(self) -> None:
         clibassist.assist_ephem_free_pointers(byref(self))
 
-    _fields_ =  [("jd_ref", c_double),
-                 ("_pl", c_void_p),
-                 ("_spl", c_void_p),
-            ]
+    _fields_ = [
+        ("jd_ref", c_double),
+        ("spk_global", c_void_p),
+        ("spk_planets", c_void_p),
+        ("spk_asteroids", c_void_p),
+    ]
 

assist/extras.py

@@ -1,13 +1,17 @@
-from typing import NoReturn, List
-
-from . import clibassist
-from ctypes import Structure, c_double, POINTER, c_int, c_uint, c_long, c_ulong, c_void_p, c_char_p, CFUNCTYPE, byref, c_uint32, c_uint, cast, c_char
-import rebound
 import warnings
-from .ephem import Ephem
+from ctypes import (CFUNCTYPE, POINTER, Structure, byref, c_char, c_char_p,
+                    c_double, c_int, c_long, c_uint, c_uint32, c_ulong,
+                    c_void_p, cast)
+from typing import List, NoReturn
+
 import numpy as np
 import numpy.typing as npt
 
+import rebound
+
+from . import clibassist
+from .ephem import Ephem
+
 ASSIST_FORCES = {
     "SUN"                : 0x01,
     "PLANETS"            : 0x02,

assist/test/test_apophis.py

@@ -1,13 +1,16 @@
-import rebound
-import assist
-import unittest
 import math
+import unittest
+
 import numpy as np
+import rebound
+
+import assist
+
 
 class TestAssist(unittest.TestCase):
     def test_apophis(self):
 
-        ephem = assist.Ephem("data/linux_p1550p2650.440", "data/sb441-n16.bsp")
+        ephem = assist.Ephem("data/de440.bsp", "data/sb441-n16.bsp")
 
         t_initial = 2.4621385359989386E+06 - ephem.jd_ref # Julian Days relative to jd_ref
 

assist/test/test_basic.py

@@ -1,7 +1,10 @@
+import math
+import unittest
+
 import rebound
+
 import assist
-import unittest
-import math
+
 
 class TestAssist(unittest.TestCase):
     def test_rebound(self):
@@ -12,18 +15,18 @@ def test_rebound(self):
         self.assertEqual(sim.t,1.0)
 
     def test_ephem(self):
-        ephem = assist.Ephem("data/linux_p1550p2650.440", "data/sb441-n16.bsp")
+        ephem = assist.Ephem("data/de440.bsp", "data/sb441-n16.bsp")
         self.assertEqual(ephem.jd_ref, 2451545.0)
         p = ephem.get_particle(0,0) # Sun
-        self.assertEqual(p.x, -0.007137179161607906)
+        self.assertEqual(p.x, -0.0071371791616079054)
         p = ephem.get_particle(1,100) # planet
         self.assertEqual(p.x, 0.12906301685045435)
         p = ephem.get_particle(20,200) #asteroid
         self.assertEqual(p.x, -2.62956381075119)
         del ephem
 
     def test_ephem_names(self):
-        ephem = assist.Ephem("data/linux_p1550p2650.440", "data/sb441-n16.bsp")
+        ephem = assist.Ephem("data/de440.bsp", "data/sb441-n16.bsp")
         p1 = ephem.get_particle(0,0) # Sun
         p2 = ephem.get_particle("Sun",0) # Also Sun
         self.assertEqual(p1.x, p2.x)
@@ -32,7 +35,7 @@ def test_ephem_names(self):
             ephem.get_particle("Planet 9",0) # Does not exist
 
     def test_holman(self):
-        ephem = assist.Ephem("data/linux_p1550p2650.440", "data/sb441-n16.bsp")
+        ephem = assist.Ephem("data/de440.bsp", "data/sb441-n16.bsp")
 
         sim = rebound.Simulation()
         extras = assist.Extras(sim, ephem)

assist/test/test_forces.py

@@ -1,12 +1,15 @@
+import math
+import unittest
+
 import rebound
+
 import assist
-import unittest
-import math
+
 
 class TestForces(unittest.TestCase):
     def test_errors(self):
         sim = rebound.Simulation()
-        ephem = assist.Ephem("data/linux_p1550p2650.440", "data/sb441-n16.bsp")
+        ephem = assist.Ephem("data/de440.bsp", "data/sb441-n16.bsp")
         extras = assist.Extras(sim, ephem)
 
         with self.assertRaises(AttributeError) as context:
@@ -33,7 +36,7 @@ def test_forces(self):
 
         sim2 = sim.copy()
 
-        ephem = assist.Ephem("data/linux_p1550p2650.440", "data/sb441-n16.bsp")
+        ephem = assist.Ephem("data/de440.bsp", "data/sb441-n16.bsp")
         extras = assist.Extras(sim, ephem)
         extras2 = assist.Extras(sim2, ephem)
 

assist/test/test_interpolate.py

@@ -1,7 +1,10 @@
+import math
+import unittest
+
 import rebound
+
 import assist
-import unittest
-import math
+
 
 class TestInterpolate(unittest.TestCase):
     def test_interpolate(self):
@@ -18,7 +21,7 @@ def test_interpolate(self):
 
         sim2 = sim.copy()
 
-        ephem = assist.Ephem("data/linux_p1550p2650.440", "data/sb441-n16.bsp")
+        ephem = assist.Ephem("data/de440.bsp", "data/sb441-n16.bsp")
         extras = assist.Extras(sim, ephem)
         extras2 = assist.Extras(sim2, ephem)
 
@@ -30,7 +33,7 @@ def test_interpolate(self):
             d = sim.particles[0] - sim2.particles[0]
             au2meter = 149597870700
 
-            self.assertLess(math.fabs(d.x*au2meter), 0.01) # 1cm accurary 
+            self.assertLess(math.fabs(d.x*au2meter), 0.02) # 1cm accurary 
             self.assertLess(math.fabs(d.y*au2meter), 0.01) # 1cm accurary 
             self.assertLess(math.fabs(d.z*au2meter), 0.01) # 1cm accurary 
 

data/readme.txt

@@ -1,7 +1,7 @@
 This is the standard location for the ephemeris files, i.e.:
-linux_m13000p17000.441 & sb441-n16.bsp
+de441.bsp & sb441-n16.bsp
 
 You can download them from: 
 
-https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de440/linux_p1550p2650.440
+https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de440.bsp
 https://ssd.jpl.nasa.gov/ftp/eph/small_bodies/asteroids_de441/sb441-n16.bsp

docs/installation.md

@@ -31,14 +31,14 @@ To use ASSIST, you also need to download ephemeris data files. One file for plan
 
 ```bash
 mkdir data
-curl https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de440/linux_p1550p2650.440 -o data/linux_p1550p2650.440
+curl https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de440.bsp -o data/de440.bsp
 curl https://ssd.jpl.nasa.gov/ftp/eph/small_bodies/asteroids_de441/sb441-n16.bsp -o data/sb441-n16.bsp
 ```
 
 For some of the examples, you will also need the planet ephemeris file with an extended coverage.
 
 ```bash
-curl https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de441/linux_m13000p17000.441 -o assist/data/linux_m13000p17000.441
+curl https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de441.bsp -o assist/data/de441.bsp
 ```
 
 !!! Info 
@@ -62,7 +62,7 @@ You can now verify that your installation of assist works.
     Then run the following code:
     ```python
     import assist
-    ephem = assist.Ephem("data/linux_p1550p2650.440", "data/sb441-n16.bsp")
+    ephem = assist.Ephem("data/de440.bsp", "data/sb441-n16.bsp")
     print(ephem.jd_ref)
     ephem.get_particle("Earth", 0)
     ```

examples/asteroid/problem.c

@@ -17,7 +17,7 @@ int main(int argc, char* argv[]){
 
     // Load the ephemeris data
     struct assist_ephem* ephem = assist_ephem_create(
-            "../../data/linux_p1550p2650.440",
+            "../../data/de440.bsp",
             "../../data/sb441-n16.bsp");
     if (!ephem){
         printf("Cannot create ephemeris structure.\n");

examples/ephemeris/problem.c

@@ -17,7 +17,7 @@ int main(int argc, char* argv[]){
     // download the files.
 
     struct assist_ephem* ephem = assist_ephem_create(
-            "../../data/linux_p1550p2650.440",
+            "../../data/de440.bsp",
             "../../data/sb441-n16.bsp");
     if (!ephem){
         printf("Cannot create ephemeris structure.\n");