B612-Asteroid-Institute · akoumjian · Jul 26, 2024
diff --git a/.github/workflows/c.yml b/.github/workflows/c.yml
@@ -29,6 +29,7 @@ jobs:
         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/bsp/de440.bsp
           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/small_bodies/asteroids_de441/sb441-n16.bsp
       - name: Reproduce Holman trajectory
@@ -91,3 +92,38 @@ jobs:
         run: |
           make
           LD_LIBRARY_PATH=../../../rebound/src/ ./rebound
+      - name: Compare Linux Binary and BSP DE440 Ephem
+        working-directory: ./unit_tests/de440_compare_ephem/
+        run: |
+          make
+          LD_LIBRARY_PATH=../../../rebound/src/ ./rebound
+      - name: Compare Linux Binary and BSP DE440 Constants
+        working-directory: ./unit_tests/de440_compare_constants/
+        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
+      - name: Reproduce Apophis Encounter with bsp 
+        working-directory: ./unit_tests/apophis_bsp/
+        run: |
+          make
+          LD_LIBRARY_PATH=../../../rebound/src/ ./rebound
+      - name: Reproduce Ceres-5303 encounter with bsp
+        working-directory: ./unit_tests/5303_Ceres_bsp/
+        run: |
+          make
+          LD_LIBRARY_PATH=../../../rebound/src/ ./rebound
diff --git a/README.md b/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. You can also manually download them using 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/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")
+
+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)
 

diff --git a/assist/ephem.py b/assist/ephem.py
@@ -1,10 +1,13 @@
+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
+import rebound
+
+from . import assist_error_messages, clibassist
 
 ASSIST_BODY_IDS = {
         0: "Sun",
@@ -82,8 +85,11 @@ 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),
+        ("jpl_planets", c_void_p),
+        ("spk_global", c_void_p),
+        ("spk_planets", c_void_p),
+        ("spk_asteroids", c_void_p),
+    ]
 
diff --git a/assist/extras.py b/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,

diff --git a/assist/test/test_apophis.py b/assist/test/test_apophis.py
@@ -1,9 +1,12 @@
-import rebound
-import assist
-import unittest
 import math
+import unittest
+
 import numpy as np
 
+import assist
+import rebound
+
+
 class TestAssist(unittest.TestCase):
     def test_apophis(self):
 

diff --git a/assist/test/test_basic.py b/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):
@@ -15,11 +18,11 @@ def test_ephem(self):
         ephem = assist.Ephem("data/linux_p1550p2650.440", "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.0071371791616045835)
         p = ephem.get_particle(1,100) # planet
-        self.assertEqual(p.x, 0.12906301685045435)
+        self.assertEqual(p.x, 0.12906301685043747)
         p = ephem.get_particle(20,200) #asteroid
-        self.assertEqual(p.x, -2.62956381075119)
+        self.assertEqual(p.x, -2.629563810751188)
         del ephem
 
     def test_ephem_names(self):

diff --git a/src/assist.c b/src/assist.c
@@ -114,7 +114,18 @@ int assist_ephem_init(struct assist_ephem* ephem, char *user_planets_path, char
         strncpy(planets_path, user_planets_path, FNAMESIZE-1);	
     }
 
-    if ((ephem->jpl = assist_jpl_init(planets_path)) == NULL) {
+    // Initialize either the linux binary or
+    // the spice kernel file for the planets.
+    // Use the extension to detect which
+    if (strstr(planets_path, ".bsp")){
+        ephem->spk_planets = assist_spk_init(planets_path);
+        ephem->spk_global = assist_load_spk_constants(planets_path);
+        assist_spk_join_masses(ephem->spk_planets, ephem->spk_global);
+    }else{
+        ephem->jpl_planets = assist_jpl_init(planets_path);
+    }
+
+    if (ephem->jpl_planets == NULL && ephem->spk_planets == NULL){
         return ASSIST_ERROR_EPHEM_FILE;	  
     }
 
@@ -132,26 +143,15 @@ int assist_ephem_init(struct assist_ephem* ephem, char *user_planets_path, char
     }
 
     if (asteroids_path_not_found != 1){
-        if ((ephem->spl = assist_spk_init(asteroids_path)) == NULL) {
+        if ((ephem->spk_asteroids = assist_spk_init(asteroids_path)) == NULL) {
             asteroids_path_not_found = 1;
         }
     }
 
+    // Join the mass data from whichever planets file we are using.
     if (asteroids_path_not_found != 1){
-        // Try to find masses of bodies in spk file in ephemeris constants
-        for(int n=0; n<ephem->spl->num; n++){ // loop over all asteroids
-            int found = 0;
-            for(int c=0; c<ephem->jpl->num; c++){ // loop over all constants
-                if (strncmp(ephem->jpl->str[c], "MA", 2) == 0) {
-                    int cid = atoi(ephem->jpl->str[c]+2);
-                    int offset = 2000000;
-                    if (cid==ephem->spl->targets[n].code-offset){
-                        ephem->spl->targets[n].mass = ephem->jpl->con[c];
-                        found = 1;
-                        break;
-                    }
-                }
-            }
+
+        if (ephem->jpl_planets) {
             // Use lookup table for new KBO objects in DE440/441
             // Source: https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/README.txt
             int massmap[] = {
@@ -187,32 +187,56 @@ int assist_ephem_init(struct assist_ephem* ephem, char *user_planets_path, char
                 8029,  2528381,
                 8030,  3515022,
             };
+
+            // Try to find masses of asteroid bodies in spk file in ephemeris constants
+            for(int n=0; n<ephem->spk_asteroids->num; n++){ // loop over all asteroids
+                int found = 0;
+                // If we're using linux binary
+                for(int c=0; c<ephem->jpl_planets->num; c++){ // loop over all constants
+                    if (strncmp(ephem->jpl_planets->str[c], "MA", 2) == 0) {
+                        int cid = atoi(ephem->jpl_planets->str[c]+2);
+                        int offset = 2000000;
+                        if (cid==ephem->spk_asteroids->targets[n].code-offset){
+                            ephem->spk_asteroids->targets[n].mass = ephem->jpl_planets->con[c];
+                            found = 1;
+                            break;
+                        }
+                    }
+                }
+
             if (found==0){
                 int mapped = -1;
                 for (int m=0; m<sizeof(massmap); m+=2){
-                    if (massmap[m+1]==ephem->spl->targets[n].code){
+                        if (massmap[m+1]==ephem->spk_asteroids->targets[n].code){
                         mapped = massmap[m];
                         break;
                     }
                 }
                 if (mapped != -1){
-                    for(int c=0; c<ephem->jpl->num; c++){ // loop over all constants (again)
-                        if (strncmp(ephem->jpl->str[c], "MA", 2) == 0) {
-                            int cid = atoi(ephem->jpl->str[c]+2);
+                        for(int c=0; c<ephem->jpl_planets->num; c++){ // loop over all constants (again)
+                            if (strncmp(ephem->jpl_planets->str[c], "MA", 2) == 0) {
+                                int cid = atoi(ephem->jpl_planets->str[c]+2);
                             if (cid==mapped){
-                                ephem->spl->targets[n].mass = ephem->jpl->con[c];
+                                    ephem->spk_asteroids->targets[n].mass = ephem->jpl_planets->con[c];
                                 found = 1;
                                 break;
                             }
                         }
                     }
                 }
             }
+
             if (found==0){
-                fprintf(stderr,"WARNING: Cannot find mass for asteroid %d (NAIF ID Number %d).\n", n, ephem->spl->targets[n].code );
+                    fprintf(stderr,"WARNING: Cannot find mass for asteroid %d (NAIF ID Number %d).\n", n, ephem->spk_asteroids->targets[n].code );
             }
 
         }
+
+        } else if (ephem->spk_global) {
+            assist_spk_join_masses(ephem->spk_asteroids, ephem->spk_global);
+        }
+
+
     }else{
         fprintf(stderr, "(ASSIST) %s\n", assist_error_messages[ASSIST_ERROR_AST_FILE]);
     }
@@ -233,13 +257,20 @@ struct assist_ephem* assist_ephem_create(char *user_planets_path, char *user_ast
 }
 
 void assist_ephem_free_pointers(struct assist_ephem* ephem){
-    if (ephem->jpl){
-        assist_jpl_free(ephem->jpl);
+    if (ephem->jpl_planets != NULL){
+        assist_jpl_free(ephem->jpl_planets);
+    }
+    if (ephem->spk_planets != NULL){
+        assist_spk_free(ephem->spk_planets);
+    }
+    if (ephem->spk_asteroids != NULL){
+        assist_spk_free(ephem->spk_asteroids);
     }
-    if (ephem->spl){
-        assist_spk_free(ephem->spl);
+    if (ephem->spk_global != NULL){
+        assist_free_spk_constants(ephem->spk_global);
     }
 }
+
 void assist_ephem_free(struct assist_ephem* ephem){
     assist_ephem_free_pointers(ephem);
     free(ephem);
@@ -278,8 +309,8 @@ void assist_init(struct assist_extras* assist, struct reb_simulation* sim, struc
     assist->sim = sim;
     assist->ephem_cache = calloc(1, sizeof(struct assist_ephem_cache));
     int N_total = ASSIST_BODY_NPLANETS;
-    if (ephem->spl){
-        N_total += ephem->spl->num;
+    if (ephem->spk_asteroids){
+        N_total += ephem->spk_asteroids->num;
     }
     assist->gr_eih_sources = 1; // Only include Sun by default
     assist->ephem_cache->items = calloc(N_total*7, sizeof(struct assist_cache_item));
@@ -578,3 +609,50 @@ static void assist_pre_timestep_modifications(struct reb_simulation* sim){
     memcpy(assist->last_state, sim->particles, sizeof(struct reb_particle)*sim->N);
 }
 
+double assist_get_constant(struct assist_ephem* ephem, const char* constant_name) {
+    // First, check if the constant is in jpl_s struct
+    if (ephem->jpl_planets != NULL) {
+        if (strcmp(constant_name, "AU") == 0) return ephem->jpl_planets->AU;
+        if (strcmp(constant_name, "EMRAT") == 0) return ephem->jpl_planets->cem;
+        if (strcmp(constant_name, "J2E") == 0) return ephem->jpl_planets->J2E;
+        if (strcmp(constant_name, "J3E") == 0) return ephem->jpl_planets->J3E;
+        if (strcmp(constant_name, "J4E") == 0) return ephem->jpl_planets->J4E;
+        if (strcmp(constant_name, "J2SUN") == 0) return ephem->jpl_planets->J2SUN;
+        if (strcmp(constant_name, "RE") == 0) return ephem->jpl_planets->RE;
+        if (strcmp(constant_name, "CLIGHT") == 0) return ephem->jpl_planets->CLIGHT;
+        if (strcmp(constant_name, "ASUN") == 0) return ephem->jpl_planets->ASUN;
+    }
+
+    // Then, check if the constant is in spk_global struct
+    if (ephem->spk_global != NULL) {
+        struct spk_constants* con = &ephem->spk_global->con;
+        if (strcmp(constant_name, "AU") == 0) return con->AU;
+        if (strcmp(constant_name, "EMRAT") == 0) return con->EMRAT;
+        if (strcmp(constant_name, "J2E") == 0) return con->J2E;
+        if (strcmp(constant_name, "J3E") == 0) return con->J3E;
+        if (strcmp(constant_name, "J4E") == 0) return con->J4E;
+        if (strcmp(constant_name, "J2SUN") == 0) return con->J2SUN;
+        if (strcmp(constant_name, "RE") == 0) return con->RE;
+        if (strcmp(constant_name, "CLIGHT") == 0) return con->CLIGHT;
+        if (strcmp(constant_name, "ASUN") == 0) return con->ASUN;
+    }
+
+    // If the constant is not found, return some error value (NaN)
+    return NAN;
+}
+
+
+// Function to truncate a double value to a specific number of bits
+double truncate_double(double value, int precision) {
+    union {
+        double d;
+        uint64_t u;
+    } u;
+    u.d = value;
+
+    // Mask to keep only the desired precision bits
+    uint64_t mask = ~((1ULL << (52 - precision)) - 1);
+    u.u &= mask;
+
+    return u.d;
+}