From aaaad63d077c403565810ad4de356e420247c3d6 Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Tue, 27 May 2025 14:08:04 -0400 Subject: [PATCH 1/6] Use .bsp for DE440 and refactor constants loading --- .github/workflows/c.yml | 29 +- .github/workflows/python.yml | 4 +- MANIFEST.in | 2 - Makefile | 32 +- README.md | 31 +- assist/ephem.py | 34 +- assist/extras.py | 16 +- assist/test/test_apophis.py | 11 +- assist/test/test_basic.py | 15 +- assist/test/test_forces.py | 11 +- assist/test/test_interpolate.py | 11 +- data/readme.txt | 4 +- docs/installation.md | 6 +- examples/asteroid/problem.c | 2 +- examples/ephemeris/problem.c | 2 +- examples/interpolation/problem.c | 2 +- jupyter_examples/5303_Ceres.ipynb | 2 +- jupyter_examples/Apophis.ipynb | 4 +- jupyter_examples/GettingStarted.ipynb | 6 +- jupyter_examples/VariationalEquations.ipynb | 2 +- jupyter_examples/publication_figures.ipynb | 4 +- setup.py | 14 +- src/Makefile | 2 +- src/assist.c | 126 ++--- src/assist.h | 30 +- src/forces.c | 113 ++-- src/forces.h | 2 +- src/get_mass | Bin 53464 -> 0 bytes src/get_mass.c | 33 -- src/planets.c | 363 ------------ src/planets.h | 74 --- src/spk.c | 523 ++++++++++++++++-- src/spk.h | 74 ++- src/tools.c | 2 +- unit_tests/5303_Ceres/problem.c | 2 +- unit_tests/apophis/problem.c | 6 +- unit_tests/benchmark/problem.c | 2 +- unit_tests/convert_simulation/problem.c | 2 +- unit_tests/ephem_cache/problem.c | 2 +- unit_tests/holman/problem.c | 2 +- unit_tests/holman_reverse/problem.c | 2 +- unit_tests/interpolation/problem.c | 2 +- .../problem.c | 2 +- unit_tests/onthefly_interpolation/problem.c | 4 +- unit_tests/roundtrip/problem.c | 4 +- unit_tests/roundtrip_adaptive/problem.c | 2 +- unit_tests/spk_init/Makefile | 46 ++ unit_tests/spk_init/problem.c | 65 +++ unit_tests/spk_join_masses/Makefile | 46 ++ unit_tests/spk_join_masses/problem.c | 64 +++ unit_tests/spk_load_constants/Makefile | 46 ++ unit_tests/spk_load_constants/problem.c | 33 ++ unit_tests/variational/problem.c | 2 +- 53 files changed, 1141 insertions(+), 779 deletions(-) delete mode 100755 src/get_mass delete mode 100644 src/get_mass.c delete mode 100644 src/planets.c delete mode 100644 src/planets.h create mode 100644 unit_tests/spk_init/Makefile create mode 100644 unit_tests/spk_init/problem.c create mode 100644 unit_tests/spk_join_masses/Makefile create mode 100644 unit_tests/spk_join_masses/problem.c create mode 100644 unit_tests/spk_load_constants/Makefile create mode 100644 unit_tests/spk_load_constants/problem.c diff --git a/.github/workflows/c.yml b/.github/workflows/c.yml index 0e67ca1..b9b9b0c 100644 --- a/.github/workflows/c.yml +++ b/.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/ @@ -96,3 +88,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 \ No newline at end of file diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml index 3f2324f..3a48b56 100644 --- a/.github/workflows/python.yml +++ b/.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 . diff --git a/MANIFEST.in b/MANIFEST.in index 38ddfd2..92bcc21 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,12 +1,10 @@ include src/assist.c include src/forces.c include src/tools.c -include src/planets.c include src/spk.c include src/assist.h include src/forces.h include src/tools.h -include src/planets.h include src/spk.h include README.md include LICENSE diff --git a/Makefile b/Makefile index de22331..5ee8242 100644 --- a/Makefile +++ b/Makefile @@ -7,14 +7,36 @@ all: libassist clean: $(MAKE) -C src clean - $(MAKE) -C doc clean + @if [ -d doc ] && [ -f doc/Makefile ]; then $(MAKE) -C doc clean; fi @-rm -f *.so - @pip uninstall -y assist - @python setup.py clean --all + @pip uninstall -y assist 2>/dev/null || true + @python setup.py clean --all 2>/dev/null || true @rm -rf assist.* .PHONY: doc doc: - cd doc/doxygen && doxygen - $(MAKE) -C doc html + @if [ -d doc/doxygen ]; then cd doc/doxygen && doxygen; fi + @if [ -d doc ] && [ -f doc/Makefile ]; then $(MAKE) -C doc html; fi + +# 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 \ No newline at end of file diff --git a/README.md b/README.md index 10f8bb3..11a9a26 100644 --- 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 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, * David M. Hernandez, Center for Astrophysics | Harvard & Smithsonian +* Alec Koumjian, Asteroid Institute | B612 Foundation, diff --git a/assist/ephem.py b/assist/ephem.py index e2dd5b6..a7c6e2c 100644 --- a/assist/ephem.py +++ b/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,24 @@ 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_planets", c_void_p), + ("spk_asteroids", c_void_p), + ("AU", c_double), + ("EMRAT", c_double), + ("J2E", c_double), + ("J3E", c_double), + ("J4E", c_double), + ("J2SUN", c_double), + ("RE", c_double), + ("CLIGHT", c_double), + ("ASUN", c_double), + # Derived constants calculated from base constants + ("Re_eq", c_double), + ("Rs_eq", c_double), + ("c_AU_per_day", c_double), + ("c_squared", c_double), + ("over_c_squared", c_double), + ] diff --git a/assist/extras.py b/assist/extras.py index 5842010..2ad3b5d 100644 --- 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 index be2d582..3bd1a4d 100644 --- a/assist/test/test_apophis.py +++ b/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 diff --git a/assist/test/test_basic.py b/assist/test/test_basic.py index fd99808..6362ddf 100644 --- 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): @@ -12,10 +15,10 @@ 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 @@ -23,7 +26,7 @@ def test_ephem(self): 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) diff --git a/assist/test/test_forces.py b/assist/test/test_forces.py index b4c97fc..9b5da61 100644 --- a/assist/test/test_forces.py +++ b/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) diff --git a/assist/test/test_interpolate.py b/assist/test/test_interpolate.py index e2a3875..83d6a4b 100644 --- a/assist/test/test_interpolate.py +++ b/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 diff --git a/data/readme.txt b/data/readme.txt index 73102e7..4104eb9 100644 --- a/data/readme.txt +++ b/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 diff --git a/docs/installation.md b/docs/installation.md index 94c1273..d453750 100644 --- a/docs/installation.md +++ b/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) ``` diff --git a/examples/asteroid/problem.c b/examples/asteroid/problem.c index d6a3c23..90e8868 100644 --- a/examples/asteroid/problem.c +++ b/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"); diff --git a/examples/ephemeris/problem.c b/examples/ephemeris/problem.c index edf6b4c..aad3aa4 100644 --- a/examples/ephemeris/problem.c +++ b/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"); diff --git a/examples/interpolation/problem.c b/examples/interpolation/problem.c index b436232..a9e0884 100644 --- a/examples/interpolation/problem.c +++ b/examples/interpolation/problem.c @@ -15,7 +15,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"); diff --git a/jupyter_examples/5303_Ceres.ipynb b/jupyter_examples/5303_Ceres.ipynb index d462bdf..a5c3fa2 100644 --- a/jupyter_examples/5303_Ceres.ipynb +++ b/jupyter_examples/5303_Ceres.ipynb @@ -55,7 +55,7 @@ "metadata": {}, "outputs": [], "source": [ - "ephem = assist.Ephem(\"../data/linux_p1550p2650.440\", \"../data/sb441-n16.bsp\")" + "ephem = assist.Ephem(\"../data/de440.bsp\", \"../data/sb441-n16.bsp\")" ] }, { diff --git a/jupyter_examples/Apophis.ipynb b/jupyter_examples/Apophis.ipynb index 24e7040..b3f8ac6 100644 --- a/jupyter_examples/Apophis.ipynb +++ b/jupyter_examples/Apophis.ipynb @@ -43,7 +43,7 @@ "metadata": {}, "source": [ "We begin by initializing the ASSIST ephemeris object. This object allows you (and ASSIST) to query the positions and velocities of Solar System objects at any time. For that you need to have downloaded the ephemeris data files from NASA JPL and store them on your computer. \n", - "- https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de440/linux_p1550p2650.440\n", + "- https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de440.bsp\n", "- https://ssd.jpl.nasa.gov/ftp/eph/small_bodies/asteroids_de441/sb441-n16.bsp" ] }, @@ -53,7 +53,7 @@ "metadata": {}, "outputs": [], "source": [ - "ephem = assist.Ephem(\"../data/linux_p1550p2650.440\", \"../data/sb441-n16.bsp\")" + "ephem = assist.Ephem(\"../data/de440.bsp\", \"../data/sb441-n16.bsp\")" ] }, { diff --git a/jupyter_examples/GettingStarted.ipynb b/jupyter_examples/GettingStarted.ipynb index 7955f16..65aa4c7 100644 --- a/jupyter_examples/GettingStarted.ipynb +++ b/jupyter_examples/GettingStarted.ipynb @@ -15,12 +15,12 @@ "\n", "The positions of the massive bodies come from two binary files, both provided by NASA JPL. The first is for the Sun, planets, and Moon, with the latest DE441 ephemeris. The other is for the asteroids, corresponding to DE441. To use assist, you need to download these files and store them on your computer. They can be found at these URLs:\n", "\n", - "- https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de441/linux_m13000p17000.441\n", + "- https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de441.bsp\n", "- https://ssd.jpl.nasa.gov/ftp/eph/small_bodies/asteroids_de441/sb441-n16.bsp\n", "\n", "Note that these are large files, more than 3 GB in total. Instead of the DE441 file, you can also download the DE440 file which is only 100 MB and covers a shorter timespan:\n", "\n", - "- https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de440/linux_p1550p2650.440\n", + "- https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de440.bsp\n", "\n", "Once you have installed ASSIST and downloaded these files, you're ready to go ahead with this notebook. " ] @@ -66,7 +66,7 @@ "metadata": {}, "outputs": [], "source": [ - "ephem = assist.Ephem(\"../data/linux_p1550p2650.440\", \"../data/sb441-n16.bsp\")" + "ephem = assist.Ephem(\"../data/de440.bsp\", \"../data/sb441-n16.bsp\")" ] }, { diff --git a/jupyter_examples/VariationalEquations.ipynb b/jupyter_examples/VariationalEquations.ipynb index ecadf25..b1480d0 100644 --- a/jupyter_examples/VariationalEquations.ipynb +++ b/jupyter_examples/VariationalEquations.ipynb @@ -39,7 +39,7 @@ "metadata": {}, "outputs": [], "source": [ - "ephem = assist.Ephem(\"../data/linux_p1550p2650.440\", \"../data/sb441-n16.bsp\")" + "ephem = assist.Ephem(\"../data/de440.bsp\", \"../data/sb441-n16.bsp\")" ] }, { diff --git a/jupyter_examples/publication_figures.ipynb b/jupyter_examples/publication_figures.ipynb index bbcc099..0a83193 100644 --- a/jupyter_examples/publication_figures.ipynb +++ b/jupyter_examples/publication_figures.ipynb @@ -50,8 +50,8 @@ "metadata": {}, "outputs": [], "source": [ - "#ephem = assist.Ephem(\"../data/linux_p1550p2650.440\", \"../data/sb441-n16.bsp\")\n", - "ephem = assist.Ephem(\"../data/linux_m13000p17000.441\", \"../data/sb441-n16.bsp\")" + "#ephem = assist.Ephem(\"../data/de440.bsp\", \"../data/sb441-n16.bsp\")\n", + "ephem = assist.Ephem(\"../data/de441.bsp\", \"../data/sb441-n16.bsp\")" ] }, { diff --git a/setup.py b/setup.py index 26078d3..bc58fa6 100644 --- a/setup.py +++ b/setup.py @@ -1,12 +1,12 @@ -from codecs import open -import os import inspect -import sys +import os +import sys +from codecs import open from distutils import sysconfig -from distutils.sysconfig import get_python_lib +from distutils.sysconfig import get_python_lib try: - from setuptools import setup, Extension + from setuptools import Extension, setup from setuptools.command.build_ext import build_ext as _build_ext except ImportError: print("Installing ASSIST requires setuptools. Do 'pip install setuptools'.") @@ -43,7 +43,7 @@ def finalize_options(self): # get site-packages dir to add to paths in case REBOUND & ASSIST installed simul in tmp dir rebdirsp = get_python_lib()+'/'#[p for p in sys.path if p.endswith('site-packages')][0]+'/' self.include_dirs.append(rebdir) - sources = [ 'src/assist.c', 'src/spk.c', 'src/planets.c', 'src/forces.c', 'src/tools.c'], + sources = [ 'src/assist.c', 'src/spk.c', 'src/forces.c', 'src/tools.c'], if not "CONDA_BUILD_CROSS_COMPILATION" in os.environ: self.library_dirs.append(rebdir+'/../') @@ -74,7 +74,7 @@ def finalize_options(self): extra_link_args.append('-Wl,-install_name,@rpath/libassist'+suffix) libassistmodule = Extension('libassist', - sources = [ 'src/assist.c','src/spk.c', 'src/planets.c', 'src/forces.c', 'src/tools.c'], + sources = [ 'src/assist.c','src/spk.c', 'src/forces.c', 'src/tools.c'], include_dirs = ['src'], library_dirs = [], runtime_library_dirs = ["."], diff --git a/src/Makefile b/src/Makefile index 4b68938..52e75ad 100644 --- a/src/Makefile +++ b/src/Makefile @@ -22,7 +22,7 @@ ifndef ASSISTGITHASH PREDEF+= -DASSISTGITHASH=$(ASSISTGITHASH) endif -SOURCES=assist.c spk.c planets.c forces.c tools.c +SOURCES=assist.c spk.c forces.c tools.c OBJECTS=$(SOURCES:.c=.o) HEADERS=$(SOURCES:.c=.h) diff --git a/src/assist.c b/src/assist.c index 9002c44..7f8e7ed 100644 --- a/src/assist.c +++ b/src/assist.c @@ -34,7 +34,6 @@ #include "rebound.h" #include "spk.h" -#include "planets.h" #include "forces.h" #define STRINGIFY(s) str(s) @@ -91,7 +90,7 @@ static struct reb_dpconst7 dpcast(struct reb_dp7 dp){ int assist_ephem_init(struct assist_ephem* ephem, char *user_planets_path, char *user_asteroids_path){ - char default_planets_path[] = "/data/linux_m13000p17000.441"; + char default_planets_path[] = "/data/de441.bsp"; char default_asteroids_path[] = "/data/sb441-n16.bsp"; ephem->jd_ref = 2451545.0; // Default jd_ref @@ -114,13 +113,11 @@ 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) { - return ASSIST_ERROR_EPHEM_FILE; - } + // Load both SPK files first + ephem->spk_planets = assist_spk_init(planets_path); int asteroids_path_not_found = 0; - if(user_asteroids_path == NULL){ if(getenv("ASSIST_DIR")==NULL){ asteroids_path_not_found = 1; @@ -132,88 +129,27 @@ 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; } } + + // Now load constants and masses from the planets file and apply them + struct spk_constants_and_masses data = assist_load_spk_constants_and_masses(planets_path); + assist_apply_spk_constants(ephem, &data); + + // Join masses to both planet and asteroid targets + if (ephem->spk_planets) { + assist_spk_join_masses(ephem->spk_planets, &data.masses, ephem->EMRAT); + } + if (ephem->spk_asteroids) { + assist_spk_join_masses(ephem->spk_asteroids, &data.masses, ephem->EMRAT); + } + + // Clean up temporary data + assist_free_spk_constants_and_masses(&data); - if (asteroids_path_not_found != 1){ - // Try to find masses of bodies in spk file in ephemeris constants - for(int n=0; nspl->num; n++){ // loop over all asteroids - int found = 0; - for(int c=0; cjpl->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; - } - } - } - // Use lookup table for new KBO objects in DE440/441 - // Source: https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/README.txt - int massmap[] = { - // ID, SPK_ID - 8001, 2136199, - 8002, 2136108, - 8003, 2090377, - 8004, 2136472, - 8005, 2050000, - 8006, 2084522, - 8007, 2090482, - 8008, 2020000, - 8009, 2055637, - 8010, 2028978, - 8011, 2307261, - 8012, 2174567, - 8013, 3361580, - 8014, 3308265, - 8015, 2055565, - 8016, 2145452, - 8017, 2090568, - 8018, 2208996, - 8019, 2225088, - 8020, 2019521, - 8021, 2120347, - 8022, 2278361, - 8023, 3525142, - 8024, 2230965, - 8025, 2042301, - 8026, 2455502, - 8027, 3545742, - 8028, 2523639, - 8029, 2528381, - 8030, 3515022, - }; - if (found==0){ - int mapped = -1; - for (int m=0; mspl->targets[n].code){ - mapped = massmap[m]; - break; - } - } - if (mapped != -1){ - for(int c=0; cjpl->num; c++){ // loop over all constants (again) - if (strncmp(ephem->jpl->str[c], "MA", 2) == 0) { - int cid = atoi(ephem->jpl->str[c]+2); - if (cid==mapped){ - ephem->spl->targets[n].mass = ephem->jpl->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 ); - } - - } - }else{ + if (asteroids_path_not_found == 1){ fprintf(stderr, "(ASSIST) %s\n", assist_error_messages[ASSIST_ERROR_AST_FILE]); } @@ -233,13 +169,14 @@ 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->spk_planets != NULL){ + assist_spk_free(ephem->spk_planets); } - if (ephem->spl){ - assist_spk_free(ephem->spl); + if (ephem->spk_asteroids != NULL){ + assist_spk_free(ephem->spk_asteroids); } } + void assist_ephem_free(struct assist_ephem* ephem){ assist_ephem_free_pointers(ephem); free(ephem); @@ -278,8 +215,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)); @@ -305,6 +242,7 @@ void assist_init(struct assist_extras* assist, struct reb_simulation* sim, struc assist->nm = 2.0; assist->nn = 5.093; assist->r0 = 1.0; + sim->integrator = REB_INTEGRATOR_IAS15; sim->gravity = REB_GRAVITY_NONE; sim->extras = assist; @@ -369,7 +307,7 @@ void assist_error(struct assist_extras* assist, const char* const msg){ } -struct reb_particle assist_get_particle_with_error(struct assist_ephem* ephem, const int particle_id, const double t, int* error){ +struct reb_particle assist_get_particle_with_error(const struct assist_ephem* ephem, const int particle_id, const double t, int* error){ struct reb_particle p = {0}; double GM = 0; int flag = assist_all_ephem(ephem, NULL, particle_id, t, &GM, &p.x, &p.y, &p.z, &p.vx, &p.vy, &p.vz, &p.ax, &p.ay, &p.az); @@ -379,7 +317,7 @@ struct reb_particle assist_get_particle_with_error(struct assist_ephem* ephem, c } -struct reb_particle assist_get_particle(struct assist_ephem* ephem, const int particle_id, const double t){ +struct reb_particle assist_get_particle(const struct assist_ephem* ephem, const int particle_id, const double t){ int error = 0; struct reb_particle p = assist_get_particle_with_error(ephem, particle_id, t, &error); if (error != ASSIST_SUCCESS){ @@ -494,8 +432,11 @@ void assist_integrate_or_interpolate(struct assist_extras* ax, double t){ double dts = copysign(1., sim->dt_last_done); if ( !(dts*(sim->t-sim->dt_last_done) < dts*t && dts*t < dts*sim->t) ){ // Integrate if requested time not in interval of last timestep + reb_simulation_integrate(sim, t); + } + double h = 1.0-(sim->t -t) / sim->dt_last_done; if (sim->t - t==0.){ memcpy(ax->current_state, sim->particles, sizeof(struct reb_particle)*sim->N); @@ -587,3 +528,4 @@ static void assist_pre_timestep_modifications(struct reb_simulation* sim){ memcpy(assist->last_state, sim->particles, sizeof(struct reb_particle)*sim->N); } + diff --git a/src/assist.h b/src/assist.h index 3f3e43c..5b46f34 100644 --- a/src/assist.h +++ b/src/assist.h @@ -88,8 +88,24 @@ enum ASSIST_BODY { struct assist_ephem { double jd_ref; - struct jpl_s* jpl; - struct spk_s* spl; + struct spk_s* spk_planets; + struct spk_s* spk_asteroids; + // SPK constants moved directly into assist_ephem + double AU; // definition of AU + double EMRAT; // Earth/Moon mass ratio + double J2E; // Other constant names follow JPL + double J3E; + double J4E; + double J2SUN; + double RE; + double CLIGHT; + double ASUN; + // Derived constants calculated from base constants + double Re_eq; // Earth radius in AU + double Rs_eq; // Sun radius in AU + double c_AU_per_day; // Speed of light in AU/day + double c_squared; // Speed of light squared in (AU/day)^2 + double over_c_squared; // 1/c^2 in (day/AU)^2 }; struct assist_cache_item { @@ -172,7 +188,7 @@ struct reb_simulation* assist_create_interpolated_simulation(struct reb_simulati void assist_integrate_or_interpolate(struct assist_extras* ax, double t); // Find particle position and velocity based on ephemeris data -struct reb_particle assist_get_particle(struct assist_ephem* ephem, const int particle_id, const double t); +struct reb_particle assist_get_particle(const struct assist_ephem* ephem, const int particle_id, const double t); /** * @brief Find particle position and velocity based on ephemeris data. @@ -182,11 +198,12 @@ struct reb_particle assist_get_particle(struct assist_ephem* ephem, const int pa * @param t The time at which to find the particle. * @param error Pointer to an integer to store the error code. */ -struct reb_particle assist_get_particle_with_error(struct assist_ephem* ephem, const int particle_id, const double t, int* error); +struct reb_particle assist_get_particle_with_error(const struct assist_ephem* ephem, const int particle_id, const double t, int* error); // Functions called from python: void assist_init(struct assist_extras* assist, struct reb_simulation* sim, struct assist_ephem* ephem); void assist_free_pointers(struct assist_extras* assist); +void assist_ephem_free_pointers(struct assist_ephem* ephem); void test_vary(struct reb_simulation* sim, FILE *vfile); @@ -200,7 +217,7 @@ struct assist_ephem* assist_ephem_create(char *planets_file_name, char *asteroid * @details This function prepares an ASSIST ephemeris structure for use. * @param ephem The assist_ephem pointer to initialize. * @param user_planets_path The path to the planets file. If NULL, the - * default path (/data/linux_m13000p17000.441) will be used. + * default path (/data/de441.bsp) will be used. * @param user_asteroids_path The path to the asteroids file. If NULL, the * default path (/data/sb441-n16.bsp) will be used. * @return ASSIST_SUCCESS if successful, otherwise an error code. @@ -208,7 +225,6 @@ struct assist_ephem* assist_ephem_create(char *planets_file_name, char *asteroid /// int assist_ephem_init(struct assist_ephem* ephem, char *user_planets_path, char *user_asteroids_path); - /** * @brief Converts an ASSIST supported simulation to a pure REBOUND simulation. * @details This function returns a new REBOUND simulation that contains the sun, the planets and all particles @@ -221,5 +237,5 @@ int assist_ephem_init(struct assist_ephem* ephem, char *user_planets_path, char * @return A new REBOUND simulation. Needs to be freed by the caller. */ /// -struct reb_simulation* assist_simulation_convert_to_rebound(struct reb_simulation* r, struct assist_ephem* ephem, int merge_moon); +struct reb_simulation* assist_simulation_convert_to_rebound(const struct reb_simulation* r, const struct assist_ephem* ephem, int merge_moon); #endif diff --git a/src/forces.c b/src/forces.c index 016b8ea..9749d5d 100644 --- a/src/forces.c +++ b/src/forces.c @@ -33,7 +33,6 @@ #include "assist.h" #include "rebound.h" #include "spk.h" -#include "planets.h" #include "forces.h" @@ -49,7 +48,6 @@ static void assist_additional_force_eih_GR(struct reb_simulation* sim, double xo void assist_additional_forces(struct reb_simulation* sim){ // implement additional_forces here - struct assist_extras* assist = (struct assist_extras*) sim->extras; struct assist_ephem* ephem = assist->ephem; if (assist->ephem_cache){ @@ -131,7 +129,6 @@ void assist_additional_forces(struct reb_simulation* sim){ FILE *eih_file = NULL; // Uncomment this line and recompile for testing. //eih_file = fopen("eih_acc.out", "w"); - if (assist->forces & ASSIST_FORCE_GR_EIH){ assist_additional_force_eih_GR(sim, xo, yo, zo, vxo, vyo, vzo, axo, ayo, azo, @@ -174,7 +171,7 @@ void assist_additional_forces(struct reb_simulation* sim){ } } -int assist_all_ephem(struct assist_ephem* ephem, struct assist_ephem_cache* ephem_cache, const int i, const double t, double* const GM, +int assist_all_ephem(const struct assist_ephem* ephem, struct assist_ephem_cache* ephem_cache, const int i, const double t, double* const GM, double* const x, double* const y, double* const z, double* const vx, double* const vy, double* const vz, double* const ax, double* const ay, double* const az @@ -204,19 +201,56 @@ int assist_all_ephem(struct assist_ephem* ephem, struct assist_ephem_cache* ephe // Get position and mass of massive body i. if(i < ASSIST_BODY_NPLANETS){ - int flag = assist_jpl_calc(ephem->jpl, jd_ref, t, i, GM, x, y, z, vx, vy, vz, ax, ay, az); - if(flag != ASSIST_SUCCESS) return(flag); + if (ephem->spk_planets){ + // Get position and mass of planet i + // Map the assist body enum to the SPK target code + struct { + int code; + int assist_code; + } spk_assist_planet_map[] = { + {10, ASSIST_BODY_SUN}, // Sun + {1, ASSIST_BODY_MERCURY}, // Mercury + {2, ASSIST_BODY_VENUS}, // Venus + {399, ASSIST_BODY_EARTH}, // Earth + {301, ASSIST_BODY_MOON}, // Moon + {4, ASSIST_BODY_MARS}, // Mars + {5, ASSIST_BODY_JUPITER}, // Jupiter + {6, ASSIST_BODY_SATURN}, // Saturn + {7, ASSIST_BODY_URANUS}, // Uranus + {8, ASSIST_BODY_NEPTUNE}, // Neptune + {9, ASSIST_BODY_PLUTO} // Pluto + }; + + int code = -1; + int array_size = sizeof(spk_assist_planet_map) / sizeof(spk_assist_planet_map[0]); + for (int j = 0; j < array_size; j++) { + if (spk_assist_planet_map[j].assist_code == i) { + code = spk_assist_planet_map[i].code; + break; + } + } + if (code == -1) { + return(ASSIST_ERROR_NEPHEM); + } + + int flag = assist_spk_calc_planets(ephem, jd_ref, t, code, GM, x, y, z, vx, vy, vz, ax, ay, az); + if(flag != ASSIST_SUCCESS) return(flag); + + }else{ + return(ASSIST_ERROR_NEPHEM); + } + }else{ // Get position and mass of asteroid i-ASSIST_BODY_NPLANETS. - int flag = assist_spk_calc(ephem->spl, jd_ref, t, i-ASSIST_BODY_NPLANETS, GM, x, y, z); + int flag = assist_spk_calc(ephem->spk_asteroids, jd_ref, t, i-ASSIST_BODY_NPLANETS, GM, x, y, z); if(flag != ASSIST_SUCCESS) return(flag); + // Translate SPK positions from heliocentric to barycentric. double GMs, xs, ys, zs; double vxs, vys, vzs, axs, ays, azs; // Not needed flag = assist_all_ephem(ephem, ephem_cache, ASSIST_BODY_SUN, t, &GMs, &xs, &ys, &zs, &vxs, &vys, &vzs, &axs, &ays, &azs); if(flag != ASSIST_SUCCESS) return(flag); - // Translate massive asteroids from heliocentric to barycentric. *x += xs; *y += ys; *z += zs; // velocities and accelerations are not needed for these // bodies @@ -292,16 +326,16 @@ static void assist_additional_force_direct(struct reb_simulation* sim, double xo ASSIST_BODY_JUPITER, ASSIST_BODY_SUN }; - int spl_num = 0; - if (ephem->spl){ - spl_num += ephem->spl->num; + int ast_num = 0; + if (ephem->spk_asteroids){ + ast_num += ephem->spk_asteroids->num; } // Direct forces from massives bodies - for (int k=0; k < ASSIST_BODY_NPLANETS + spl_num; k++){ + for (int k=0; k < ASSIST_BODY_NPLANETS + ast_num; k++){ int i; // ordered index - if (k>=spl_num){ - i = order[k-spl_num]; // planets and sun last + if (k>=ast_num){ + i = order[k-ast_num]; // planets and sun last }else{ i = k + ASSIST_BODY_NPLANETS; // asteroids first } @@ -321,6 +355,7 @@ static void assist_additional_force_direct(struct reb_simulation* sim, double xo return; } + // Loop over test particles for (int j=0; j=spl_num){ - i = order[k-spl_num]; // planets and sun last + if (k>=ast_num){ + i = order[k-ast_num]; // planets and sun last }else{ i = k + ASSIST_BODY_NPLANETS; // asteroids first } @@ -462,13 +497,15 @@ static void assist_additional_force_earth_J2J4(struct reb_simulation* sim, doubl assist_all_ephem(ephem, assist->ephem_cache, ASSIST_BODY_EARTH, t, &GM, &xe, &ye, &ze, &vxe, &vye, &vze, &axe, &aye, &aze); const double GMearth = GM; + double xr, yr, zr; //, vxr, vyr, vzr, axr, ayr, azr; xr = xe; yr = ye; zr = ze; - const double J2e = ephem->jpl->J2E; - const double J3e = ephem->jpl->J3E; - const double J4e = ephem->jpl->J4E; - const double Re_eq = ephem->jpl->RE/ephem->jpl->AU; + double J2e = ephem->J2E; + double J3e = ephem->J3E; + double J4e = ephem->J4E; + double Re_eq = ephem->Re_eq; + // Unit vector to equatorial pole at the epoch // Note also that the pole orientation is not changing during @@ -597,7 +634,6 @@ static void assist_additional_force_earth_J2J4(struct reb_simulation* sim, doubl int tp = vc.testparticle; struct reb_particle* const particles_var1 = particles + vc.index; if(tp == j){ - // Variational particle coords const double ddx = particles_var1[0].x; const double ddy = particles_var1[0].y; @@ -637,6 +673,7 @@ static void assist_additional_force_earth_J2J4(struct reb_simulation* sim, doubl } } } + } static void assist_additional_force_solar_J2(struct reb_simulation* sim, double xo, double yo, double zo, FILE *outfile){ @@ -661,9 +698,8 @@ static void assist_additional_force_solar_J2(struct reb_simulation* sim, double assist_all_ephem(ephem, assist->ephem_cache, ASSIST_BODY_SUN, t, &GM, &xr, &yr, &zr, &vxr, &vyr, &vzr, &axr, &ayr, &azr); const double GMsun = GM; - const double au = ephem->jpl->AU; - const double Rs_eq = ephem->jpl->ASUN/au; - const double J2s = ephem->jpl->J2SUN; + double Rs_eq = ephem->Rs_eq; + double J2s = ephem->J2SUN; // Hard-coded constants. BEWARE! double RAs = 286.13*M_PI/180.; @@ -737,7 +773,6 @@ static void assist_additional_force_solar_J2(struct reb_simulation* sim, double int tp = vc.testparticle; struct reb_particle* const particles_var1 = particles + vc.index; if(tp == j){ - // Variational particle coords double ddx = particles_var1[0].x; double ddy = particles_var1[0].y; @@ -767,11 +802,8 @@ static void assist_additional_force_solar_J2(struct reb_simulation* sim, double particles_var1[0].az += daz; } - } - } - } static void assist_additional_force_non_gravitational(struct reb_simulation* sim, @@ -1069,10 +1101,7 @@ static void assist_additional_force_potential_GR(struct reb_simulation* sim, const double jd_ref = ephem->jd_ref; // Nobili and Roxburgh GR treatment - - const double au = ephem->jpl->AU; - const double c = (ephem->jpl->CLIGHT/au)*86400; - const double C2 = c*c; + const double C2 = ephem->c_squared; const unsigned int N = sim->N; // N includes real+variational particles const unsigned int N_real = N - sim->N_var; @@ -1176,11 +1205,7 @@ static void assist_additional_force_simple_GR(struct reb_simulation* sim, const double jd_ref = ephem->jd_ref; // Damour and Deruelle solar GR treatment - - const double au = ephem->jpl->AU; - const double c = (ephem->jpl->CLIGHT/au)*86400; - - const double C2 = c*c; + const double C2 = ephem->c_squared; const unsigned int N = sim->N; // N includes real+variational particles const unsigned int N_real = N - sim->N_var; @@ -1311,10 +1336,8 @@ static void assist_additional_force_eih_GR(struct reb_simulation* sim, // This is one of three options for GR. // This version is only rarely needed. - const double au = ephem->jpl->AU; - const double c = (ephem->jpl->CLIGHT/au)*86400; - const double C2 = c*c; - const double over_C2 = 1./(c*c); + const double C2 = ephem->c_squared; + const double over_C2 = ephem->over_c_squared; const unsigned int N = sim->N; // N includes real+variational particles const unsigned int N_real = N - sim->N_var; @@ -2011,10 +2034,8 @@ static void assist_additional_force_eih_GR_orig(struct reb_simulation* sim, // This is one of three options for GR. // This version is only rarely needed. - const double au = ephem->jpl->AU; - const double c = (ephem->jpl->CLIGHT/au)*86400; - const double C2 = c*c; // This could be stored as C2. - const double over_C2 = 1./(c*c); + const double C2 = ephem->c_squared; + const double over_C2 = ephem->over_c_squared; const unsigned int N = sim->N; // N includes real+variational particles const unsigned int N_real = N - sim->N_var; diff --git a/src/forces.h b/src/forces.h index 5516e71..59d0240 100644 --- a/src/forces.h +++ b/src/forces.h @@ -29,7 +29,7 @@ // It includes all forces, including gravity. void assist_additional_forces(struct reb_simulation* sim); -int assist_all_ephem(struct assist_ephem* ephem, struct assist_ephem_cache* cache, const int i, const double t, double* const GM, +int assist_all_ephem(const struct assist_ephem* ephem, struct assist_ephem_cache* cache, const int i, const double t, double* const GM, double* const x, double* const y, double* const z, double* const vx, double* const vy, double* const vz, double* const ax, double* const ay, double* const az diff --git a/src/get_mass b/src/get_mass deleted file mode 100755 index 07ee30ecd9a4462796b2c9b3412116ddfcaaf2c4..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 53464 zcmeHQ4Rlo1wZ0P)jR@Wef~KvOQ4SD4CeiJ2=7k_zM~`gvktv=C7F<43=sz znaXuC5?c_Sg*LUe(mt)Ug^E^T5`GdWXb4IXtwuyU5w#HurlOhm?Q_qa+zimNy4tng zdh6WPbI#uT?7h!9-}%Wt_saCOFFzSLSCUM_BuPp@G9&dTNm4+vfRdzdASp;rXQB0W z+ugSRx|5_sAJb6wv>f@#Km$(aEL+j6p(@Nicc@I_Jh~!LI)_&~oj!NHk25IZ{k41x zgMs|lmuWi0Cl0C;(_mzsPOs0u&>IuQ`@7|4&EKP3L+5paWbSXUf8@_u?yjtL*Fh=X zUnlo>f;%9daF0Ew^v_jO?HwF&yua`*TK|f;1LB8_H|AI)oKBb5S65NJh|A;srB2iQ z74mxGjn9z2qQdF)SKnVzUFNK)_S7(b==_#*e=Xbz@$o*4{%M^=n4f&>EUWX@J7>() zR!gqdFMVCam$ZJ3YuM>rR5OSg_4D-pQenrQT+N)RxW#+racbGy7oO}V7G{<$sVoCe4~`$!j~6a}W}W}( z8y5bwyJY;%f+t^}2)P9bU6v$yxO5H@(QIy#OB1dJ&7mL-DFfrD@z8l@e5l^g8s~0;`O(zbY0&3D1>CYTqVqM)XQKe6^_qX){ekINI#r25AK zb>}fTCDibBG-_4%Sk*VI>KBUoNBrtPo@ zDPJ|dN6;ntM{g$RUP?r!1a|o!?0uRjb+MiXmy(A+`|eK0re7}|m2rX6cEwWb}^6(!_L01jcKfdv20lu%&;TK8Q8b1VpR z(;FJv^O8Oo%T^+qh6~vPv$E@D! zmZJ6SZd&kP95~Tt@?CPKE-H|d3C!In%b%4ybE%hJv8}MC1)qd4oA)xD8fq4JGqR* z2v)acciWl{$pIS%WwUh1fvF(!Ew9OeNg&kSjfcNdLJxeUr~^vdrwK}Qr&aBmi zRWthYEUnt26a!lPhN@= zdWGaK(TS|6JC9zdgsy9VhEX0dVPKBn&AI5jJva>-erVXTk0_ya%&B@PI~;kKT&nM8 zp8)+^(1EbLd^{GKRo$r>-2o%D9koTM-8#jJZq$LYtJ}8_qo^H`(c~kvDHlD;#~;v} zxSMW#J?Cw`MBlYhr{Hfe`7!cJ96ts8AEsPJlOkauU-A&-;XhWdJp?(-iK&3t8OH0;;kv{xK1SB!XB=0~{5L%j+!46(#I1t&z)ly?>FK~$^ zCiyv&yqtD~;0|4PfO6=tN%v8vkQ$XW2u=X_W<`XVcb0qBk3v z4Ktk~NNfh)p?(^=Qbkmrp#MD0G5?@SmkV!VC zZ7>>c<7!i5YQN!<`k3S=+;F4O@KLUo8&j+2lBF?;hZ}A-8s5d#@?&b#xLSiz?P@N3 z*(e;tZIzhqFs??&YE1Yk{u32@Xk%K~D14I%eQA!ES_fBKZdCg%7j8ES*Kymzm~D`& zH5=6yb77}ZSjuf@#cc27YAr^!8@X_gQJBtci(-Ly-{aUGODEB;pnd#@ z7O?Wz2B=E@C{8Q%$OD}qA3_DlE+SD0 z@*t5=utpte3n(k2tQ6_+9vx}*q(D1ZSLlmV zSV{_;4GMo7qR^by0G@PcN9s=cHF`I?fmWK+Xd9AF*r4-OoDOX_(xDxu`=Nsqbj%NJ zE7GBzr0XQc>9lYi+A+FLa-0rr6Y|q#&>0b@L)(FLXb0$i=v+7UhZZ^M(7M)jM#t&2 zavfTMx=w1G4lO(K(__$~zZ+Ux`JHCFs;!Ld=t~ewFH&Yx`)t{}BGQGhMRgYGBN(W* z2RE`+gyo;z9oaRCiRdMRUR$7kU|+Ul3UssI^0erG(P)2lk>b$(MH|QsF?nR` zGx#!VhV>!hN6pa-NyiGxq5dUu2TBSkf=T;F!x@jZ#yHd-1Y8t583L}zm$-W)g7n=; zLDtMdMLnqmU!Yc$;8q$-AK8 zXx%a;7OGweuXY3-ffI7mH>n}@B!YbzF915_N6v!))UC<#BLlFqs;?@+poZx#x#=hq zy^afjEsR?bUBd1#gDdFLuA3$MjydG z=aNv1CyNnwLUKEfos57MQWhirB$1g!x$x+bixHQS$O5-4AjME@8|X@OdoHy^SGJVf z72vyUJ3(idF75qE3}3wMX#CiWiT((aq`k~Ih~S^*)2n{AY(TzBk{+tCtGGCfv5-pW zdtF5UxoEyj=iy=OIlN!)Wfa;Hx8=#-i$29yPd)1<6-MZ z>B!4PM{(84o5IU1US{yp!pl@%qNhxs1_(Qf{G|EY5U&zm1e9IvQe3we;W*qGe1Q&{ z;A%R8f-Aa^1Nl46_E3G6qE1aK2q~!)LqGKix{FcFH!N?hfAqeK-Zjw#*;QO+Q*1r7 zswGR#k|t9z{(|bA;Hch)kiDT({1$ho*J} ze@7=}@P{pAAEXl$og?Fd^-9PV*34`zv{rf|(T{0OWvpUr3k&BWCCFBzUA+xk!#r{| zz(skF1$S$D9g*}Hp^Y~_*bsQfGU{!QJa2xB2k%vAQGOe@aCfUmk*l#+ie9IcI0sAD zyo$O_3UgB7zZ?wwXT<#K{-8O3J)C|Fr+A%I^bwgCDK27u@vFHF{0ZNnM+f8`~o6tsswAhR!!8gXWt}##~=; z1)q#lCwtS;Qzx6fuhjS3zjMA7YV2Q`ZyAQbE~L-*1w-e1CeSk{JvEZeF@G^tg10H7 z4(ap!CZS35Je$okfw`N3N%PE>Cxbg~o@ZLnXw08J&xhgkeaz5%Wa6QDp3UZY&?(Jx zNT25!hIuwqPmAzUt^_M>R688e?Y*dtJl#ffJ)6ySm?@DgZmtIp4DGy#?e8pYzK_u6 znOhlLv0mFhr_T3|-hVLPgXhI$&0mE1s|~Jr_Duf#XXkIK=C6b{#G|ozJh|v)4F{4G%or19h(nD^GAVZa5aUX z!4=^%`I~qa{v4V=3QB{kmotA9te&=ihyQZwkmLJTsQIG+Hn_T(`RhEBzaPZ>>E9lx zWPfkjAu(sBT^(rL+CcyM=fP~mcI^DH2VZ`N{v-LwZ;~mB1Ct2mX)r9(`#!1DP3lkf z;Fjf}?ZMx-poshy8e7(Yq7NW8FjQ8dqz@q_vB2xpcnj`XAeYoV9F+EvgS|`gz?j}8 zJ&fMqU6P`%+i)@(-AZjlChD!>(?lnv^!Kley)4>f0ZS3)7huW!8kpnNQWwN}>>9Zk z=1OCKvK$zR+Rgda+cGBwNaMxKnT2;u3hf7ZUYi^^g1^hK-AuOVC%YcLr8V=iQPu%j zFhC-;E!x)BJ2sS`Y}(g9EI^ix{Yi4*NqE{qP--Ez2L_JN&snFqlwLtS1abDz_n{3h zCy|y^e-o&EjLHVp&8VwDQD3F{OYu@ZQMs9X`S(y8C$xzBYrOh8x)$p#a1eU*mKzZ( zMm8M%ksnP&p094&jw)OB2{!-7>lx+53E9eCcjRyU19dLPp**^NW?kBCA$slGW%FU__}@VJF0DDpm9RMc(w zJ40QYuCDE;#cWV$CxsyE_;eNSVKv>=yj6lYbeAblUCGG2;L7ymuP~jXgG7#C7nR%&)E;&u|aY3Yite5B_PlVY>74p4s*gc zxbZcbGi8n9n*9A^HhFyUu|6dTJ+x^JeRn`{MU{$o=ggjs1NU8}?XT06wxbr)NxEJ* zj`dJ6>sW-={IB5!>mTddoW7gUMSIXZ>WX_IM@&?LubJbIJ6t1{-5JW?-LJS_wY#?A zvQM$4c?8R6btY>*iuS+`Q7N-)SGLe#B#Ax(Z_>UE&?#`y$Z@Qplt4e$7=85ZqYr$p zla2l7%gyw+MH_z+LRd7cPdAd2fQpSH7YnNc83O>Ek`@b6?8J!&RjN<0s5^b!ya&JE+>#PJ3MN8?layp=mNUx*jwl;H*`Yl>pIl)O|W ziqBg8O96Yy$a~nxN|LmaOrQX8OY)NeD=CK(|L!n#o=tnxVG-+wp+O685Pd2c`4L-8 zB}rPDLN+OR7DI~`1>ZT_dJ+{AUIC!pd_GS4XtqoL-b4-c(Rev_Y`oLW2C7^KrP z*r9R=O>2irpMN+G724@|Qxh2hhtc$x^C;iK@-WNWS-zd+J6Ybz^4%=o!}3>I-o^56 zmcP#O11vwt@4&AF@1x+!4C;GZfaD(e#`I{C5cDG~V71u?Wj? zN0oZtCJG-6rQU9$t|zLKs9TBJLeyPEZ6vCMs9zCvFHx;TEg|X&q5?!U6ZIrf4->VX zsCuF{5mimpcB0CN>L%)LqCO&uUhtHkAnG=vlCWcXrxG=asEI^fMpOn-6NvgYQPYST zP1K!44F?4uL!Q2A8r!#(m8r*bHQvli?cY~oC(#MW5{icL8}qClk7xe;H-2h;Z{)dq z_mBKhZhF~y>F(ju^g@TT&|#fno8`0>Dz*ZtvZC5w@2tw6l$DiLn|%ZF36m#hOEtCb zYRgsW-m8X7o-#{%*>I`Anw~YeOrqEgF&iDx(j`bZmZdBtdZ@%)(PP?Ilw>kVCR3?t zm?`D_L=*ixPX!f}rrk^`8cqEkAzdKlGmlbY18l>xlwdYXrrH!Ku@G;A`bJ5^TJXBE zXS6h|0iO`?kNA@szJYptC*5$VlyI4pK;1d3rw89JNI7eQRyOd`g3gn?mY3#<8oz*- z8+a*A(&UA_Y~^JiFLNhr^=H)?p(GFx2nYlO0s;YnfIvVXAP^7;2m}NI0s(=5KtLcM z5D*9m1Ox&C0fB%(Kp-Fx5C{ka1Ofs9fq+0jARrJB2nYlO0s;YnfIvVXAP^7;2m}NI z0s(=5KtLcM5D*9m1Ox&C0fB%(Kp-Fx5C{ka1Oftq|ECBX#BU4GzvDm5%a3_^f|t@Q zTK#Zdj^^dJd6~h>iM*W3%iDNa$jiHVS;ot1Ue@#SVO}=#@(EtH^72=_+{nu_O$UvI zeq(`tQk{NV0ly(5;a7p!|7ZAlDDfj+?EeFte!W9H0s(=5KtLcM5D*9m1Ox&C0fB%( zKp-Fx5C{ka1Ofs9fq+0jARrJB2nYlO0s;YnfIvVXAP^7;2m}NI0s(=5KtLcM5D*9m z1Ox&C0fB%(Kp-Fx5C{ka1Ofs9fq+0jARrJB2nYlO0s;YnfIvVXAQ1S!fp!v=OTY;U72)Ooz$Yq%WC|-S}3Gd?{yn7?l_XDvymVnrosSp`7*}a!WsQrt4B{$lyl!H$ zz<8=^%7~sqbZuQlwa-IrF0plPrlt^$QTWpAGw#f{T|b!=3L#wN_PMJU6ITST((88L zN6Z2+Ri$N%E3{V18LM}omN<8n*X?uGmim0|y6Wuf$*7*>Ri%JK@-~==@mF(pIk7dh z?rN%QA*MD4!5U)g+@)*`t;Be%2S>4<34C?ct}5nuLyWIo0%!OU1}A>kAjw%(?sV6? zU4Eb2Sq@*Afyq^+71h#|duuD5(#_u5`<&9^g*9bMq-Ra^`v=3MmK*Sk2=~|2-6zej zsIKry3mI680JKZ0tSa?-rSCv#3N-s}MbLj~byod&lb)KpS{~jcIbo(C0X*}8M>w7o1 z*LmIE*k5Yj$dkRkAG)}`&Y!KlzJCsId!5(qTZzC!1+{PFseOIFz0K`)p7d$_#{QnJ z|N8zs%P$Ty7u+`9O3pluiO8#_VxY$DYw`8v$faHi@$Muoj+T9{e1b7+w1(< z+Uw^LPFlAAb^dJa_4Df-Zm;uaYpT1y1!k+=u51yx2a$Y}=vpBDx?>BN@KQHq+umAm+!Fm1f z+bquO{V&n&`S{9peLg>RoY(PQ3Lf2~uCphA+{ix`$G67uzl-Cy#PK`g_`Pxb>v8;{ zIR5=OzBi8VkK=jgW8*nj8%&Hh$MLCg{MB)MRvbSij=v?2zb%e;#PM_DcxN169>>#s zd>`o%q(_m0NXIs*(tjg8iS#3+wMgrbmLsW1 zKSugbwgdN)MRJ}dCU>$3+(|Jxi<%k5;%-K=5&~X_-EbpW!Ub7?0LyYm%AD6)0m#8+Qc zf1P(>*5dmY)lQk@&019Z!2MZx7FB1?%G7PIy}oo(_H~md)JSlJd9ST3t#8h5r#i*+#&M{a0{;*|n2%JU5hP(SX=ipsOdN Qv;i?$Y -#include -#include -#include -#include "planets.h" -#include "spk.h" - -int main(int argc, char **argv) -{ - struct jpl_s *pl; - double m; - int n; - - if ((pl = jpl_init()) == NULL) - abort(); - - for (n = 0; n < 1000; n++) { - m = jpl_mass(pl, n); - - if (m > 0.0) - fprintf(stdout, "%d\t%e\n", n, m); - } - - jpl_free(pl); - return 0; -} diff --git a/src/planets.c b/src/planets.c deleted file mode 100644 index 4bd492f..0000000 --- a/src/planets.c +++ /dev/null @@ -1,363 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include "spk.h" -#include "planets.h" -#include "assist.h" - -/* - * assist_jpl_work - * - * Interpolate the appropriate Chebyshev polynomial coefficients. - * - * ncf - number of coefficients per component - * ncm - number of components (ie: 3 for most) - * niv - number of intervals / sets of coefficients - * - */ - -void assist_jpl_work(double *P, int ncm, int ncf, int niv, double t0, double t1, double *u, double *v, double *w) -{ - double T[24], S[24]; - double U[24]; - double t, c; - int p, m, n, b; - - // adjust to correct interval - t = t0 * (double)niv; - t0 = 2.0 * fmod(t, 1.0) - 1.0; - c = (double)(niv * 2) / t1 / 86400.0; - - b = (int)t; - - // set up Chebyshev polynomials and derivatives - T[0] = 1.0; T[1] = t0; - S[0] = 0.0; S[1] = 1.0; - U[0] = 0.0; U[1] = 0.0; U[2] = 4.0; - - for (p = 2; p < ncf; p++) { - T[p] = 2.0 * t0 * T[p-1] - T[p-2]; - S[p] = 2.0 * t0 * S[p-1] + 2.0 * T[p-1] - S[p-2]; - } - for (p = 3; p < ncf; p++) { - U[p] = 2.0 * t0 * U[p-1] + 4.0 * S[p-1] - U[p-2]; - } - - // compute the position/velocity - for (m = 0; m < ncm; m++) { - u[m] = v[m] = w[m] = 0.0; - n = ncf * (m + b * ncm); - - for (p = 0; p < ncf; p++) { - u[m] += T[p] * P[n+p]; - v[m] += S[p] * P[n+p] * c; - w[m] += U[p] * P[n+p] * c * c; - } - } -} - -/* - * assist_jpl_init - * - * Initialise everything needed ... probaly not be compatible with a non-430 file. - * - */ - -static double getConstant(struct jpl_s* jpl, char* name){ - for (int p = 0; p < jpl->num; p++) { - if (strncmp(name,jpl->str[p],6)==0){ - return jpl->con[p]; - } - } - fprintf(stderr,"WARNING: Constant [%s] not found in ephemeris file.\n",name); - return 0; -} - -struct jpl_s * assist_jpl_init(char *str) -{ - struct stat sb; - ssize_t ret; - int fd; - - if ((fd = open(str, O_RDONLY)) < 0){ - return NULL; - } - - if (fstat(fd, &sb) < 0){ - close(fd); - fprintf(stderr, "Error while trying to determine filesize.\n"); - return NULL; - } - - - // skip the header and constant names for now - if (lseek(fd, 0x0A5C, SEEK_SET) < 0){ - close(fd); - fprintf(stderr, "Error while seeking to header.\n"); - return NULL; - } - - struct jpl_s* jpl = calloc(1, sizeof(struct jpl_s)); - - // read header - ret = read(fd, &jpl->beg, sizeof(double)); // Start JD - ret += read(fd, &jpl->end, sizeof(double)); // End JD - ret += read(fd, &jpl->inc, sizeof(double)); // Days per block - ret += read(fd, &jpl->num, sizeof(int32_t)); // Number of constants - ret += read(fd, &jpl->cau, sizeof(double)); // AU to km - ret += read(fd, &jpl->cem, sizeof(double)); // Ratio between Earth/Moon - - // number of coefficients for all components - for (int p = 0; p < JPL_N; p++){ - jpl->ncm[p] = 3; - } - // exceptions: - jpl->ncm[JPL_NUT] = 2; // nutations - jpl->ncm[JPL_TDB] = 1; // TT-TDB - - for (int p = 0; p < 12; p++) { // Columns 1-12 of Group 1050 - ret += read(fd, &jpl->off[p], sizeof(int32_t)); - ret += read(fd, &jpl->ncf[p], sizeof(int32_t)); - ret += read(fd, &jpl->niv[p], sizeof(int32_t)); - } - - ret += read(fd, &jpl->ver, sizeof(int32_t)); // Version. e.g. 440 - ret += read(fd, &jpl->off[12], sizeof(int32_t)); // Columns 13 of Group 1050 - ret += read(fd, &jpl->ncf[12], sizeof(int32_t)); - ret += read(fd, &jpl->niv[12], sizeof(int32_t)); - - // Get all the constant names - jpl->str = calloc(jpl->num, sizeof(char *)); - - // retrieve the names of the first 400 constants - lseek(fd, 0x00FC, SEEK_SET); - for (int p = 0; p < 400; p++) { // Group 1040 - jpl->str[p] = calloc(8, sizeof(char)); - read(fd, jpl->str[p], 6); - } - - // read the remaining constant names - lseek(fd, 0x0B28, SEEK_SET); - for (int p = 400; p < jpl->num; p++) { - jpl->str[p] = calloc(8, sizeof(char)); - read(fd, jpl->str[p], 6); - } - - for (int p = 13; p < 15; p++) { // Columns 14 and 15 of Group 1050 - ret += read(fd, &jpl->off[p], sizeof(int32_t)); - ret += read(fd, &jpl->ncf[p], sizeof(int32_t)); - ret += read(fd, &jpl->niv[p], sizeof(int32_t)); - } - - // adjust for correct indexing (ie: zero based) - for (int p = 0; p < JPL_N; p++){ - jpl->off[p] -= 1; - } - - // save file size, and determine 'kernel size' or 'block size' (=8144 bytes for DE440/441) - jpl->len = sb.st_size; - jpl->rec = sizeof(double) * 2; - - for (int p = 0; p < JPL_N; p++){ - jpl->rec += sizeof(double) * jpl->ncf[p] * jpl->niv[p] * jpl->ncm[p]; - } - - // memory map the file, which makes us thread-safe with kernel caching - jpl->map = mmap(NULL, jpl->len, PROT_READ, MAP_SHARED, fd, 0); - - if (jpl->map == NULL){ - close(fd); - free(jpl); // note constants leak - fprintf(stderr, "Error while calling mmap().\n"); - return NULL; - } - - // Read constants - jpl->con = calloc(jpl->num, sizeof(double)); - lseek(fd, jpl->rec, SEEK_SET); // Starts at offset of 1 block size - for (int p = 0; p < jpl->num; p++){ - read(fd, &jpl->con[p], sizeof(double)); - //printf("%6d %s %.5e\n",p,jpl->str[p],jpl->con[p]); - } - - - // Find masses - jpl->mass[0] = getConstant(jpl, "GMS "); // Sun - jpl->mass[1] = getConstant(jpl, "GM1 "); // Mercury - jpl->mass[2] = getConstant(jpl, "GM2 "); - double emrat = getConstant(jpl, "EMRAT "); // Earth Moon Ratio - double gmb = getConstant(jpl, "GMB "); // Earth Moon combined - jpl->mass[3] = (emrat/(1.+emrat)) * gmb; // Earth - jpl->mass[4] = 1./(1+emrat) * gmb; // Moon - jpl->mass[5] = getConstant(jpl, "GM4 "); // Mars - jpl->mass[6] = getConstant(jpl, "GM5 "); // Jupiter - jpl->mass[7] = getConstant(jpl, "GM6 "); - jpl->mass[8] = getConstant(jpl, "GM7 "); - jpl->mass[9] = getConstant(jpl, "GM8 "); - jpl->mass[10] = getConstant(jpl, "GM9 "); // Pluto - - - // Other constants - jpl->J2E = getConstant(jpl, "J2E "); - jpl->J3E = getConstant(jpl, "J3E "); - jpl->J4E = getConstant(jpl, "J4E "); - jpl->J2SUN = getConstant(jpl, "J2SUN "); - jpl->AU = getConstant(jpl, "AU "); - jpl->RE = getConstant(jpl, "RE "); - jpl->CLIGHT = getConstant(jpl, "CLIGHT"); - jpl->ASUN = getConstant(jpl, "ASUN "); - - // this file descriptor is no longer needed since we are memory mapped - if (close(fd) < 0) { - fprintf(stderr, "Error while closing file.\n"); - } -#if defined(MADV_RANDOM) - if (madvise(jpl->map, jpl->len, MADV_RANDOM) < 0){ - fprintf(stderr, "Error during madvise.\n"); - } -#endif - - return jpl; - -} - -void assist_jpl_free(struct jpl_s *jpl) { - if (jpl == NULL){ - return; - } - if (munmap(jpl->map, jpl->len) < 0){ - fprintf(stderr, "Error during munmap().\n"); - } - for (int p = 0; p < jpl->num; p++){ - free(jpl->str[p]); - } - free(jpl->str); - free(jpl->con); - free(jpl); -} - -/* - * jpl_calc - * - * Caculate the position+velocity in _equatorial_ coordinates. - * Assumes pos is initially zero. - */ -enum ASSIST_STATUS assist_jpl_calc(struct jpl_s *jpl, double jd_ref, double jd_rel, int body, - double* const GM, - double* const out_x, double* const out_y, double* const out_z, - double* const out_vx, double* const out_vy, double* const out_vz, - double* const out_ax, double* const out_ay, double* const out_az){ - double t, *z; - u_int32_t blk; - - if (jpl == NULL || jpl->map == NULL) - return ASSIST_ERROR_EPHEM_FILE; - if(body<0 || body >= ASSIST_BODY_NPLANETS) - return(ASSIST_ERROR_NEPHEM); - - struct mpos_s pos; - - // Get mass, position, velocity, and mass of body i in barycentric coords. - *GM = jpl->mass[body]; - - // check if covered by this file - if (jd_ref + jd_rel < jpl->beg || jd_ref + jd_rel > jpl->end) - return ASSIST_ERROR_COVERAGE; - - // compute record number and 'offset' into record - blk = (u_int32_t)((jd_ref + jd_rel - jpl->beg) / jpl->inc); - z = (double*)jpl->map + (blk + 2) * jpl->rec/sizeof(double); - t = ((jd_ref - jpl->beg - (double)blk * jpl->inc) + jd_rel) / jpl->inc; - - switch (body) { // The indices in the jpl-> arrays match the JPL component index for the body - case ASSIST_BODY_SUN: - assist_jpl_work(&z[jpl->off[10]], jpl->ncm[10], jpl->ncf[10], jpl->niv[10], t, jpl->inc, pos.u, pos.v, pos.w); - break; - case ASSIST_BODY_MERCURY: - assist_jpl_work(&z[jpl->off[JPL_MER]], jpl->ncm[JPL_MER], jpl->ncf[JPL_MER], jpl->niv[JPL_MER], t, jpl->inc, pos.u, pos.v, pos.w); - break; - case ASSIST_BODY_VENUS: - assist_jpl_work(&z[jpl->off[JPL_VEN]], jpl->ncm[JPL_VEN], jpl->ncf[JPL_VEN], jpl->niv[JPL_VEN], t, jpl->inc, pos.u, pos.v, pos.w); - break; - case ASSIST_BODY_EARTH: - { - struct mpos_s emb, lun; - assist_jpl_work(&z[jpl->off[JPL_EMB]], jpl->ncm[JPL_EMB], jpl->ncf[JPL_EMB], jpl->niv[JPL_EMB], t, jpl->inc, emb.u, emb.v, emb.w); // earth moon barycenter - assist_jpl_work(&z[jpl->off[JPL_LUN]], jpl->ncm[JPL_LUN], jpl->ncf[JPL_LUN], jpl->niv[JPL_LUN], t, jpl->inc, lun.u, lun.v, lun.w); - - vecpos_set(pos.u, emb.u); - vecpos_off(pos.u, lun.u, -1.0 / (1.0 + jpl->cem)); - - vecpos_set(pos.v, emb.v); - vecpos_off(pos.v, lun.v, -1.0 / (1.0 + jpl->cem)); - - vecpos_set(pos.w, emb.w); - vecpos_off(pos.w, lun.w, -1.0 / (1.0 + jpl->cem)); - } - break; - case ASSIST_BODY_MOON: - { - struct mpos_s emb, lun; - assist_jpl_work(&z[jpl->off[JPL_EMB]], jpl->ncm[JPL_EMB], jpl->ncf[JPL_EMB], jpl->niv[JPL_EMB], t, jpl->inc, emb.u, emb.v, emb.w); - assist_jpl_work(&z[jpl->off[JPL_LUN]], jpl->ncm[JPL_LUN], jpl->ncf[JPL_LUN], jpl->niv[JPL_LUN], t, jpl->inc, lun.u, lun.v, lun.w); - - vecpos_set(pos.u, emb.u); - vecpos_off(pos.u, lun.u, jpl->cem / (1.0 + jpl->cem)); - - vecpos_set(pos.v, emb.v); - vecpos_off(pos.v, lun.v, jpl->cem / (1.0 + jpl->cem)); - - vecpos_set(pos.w, emb.w); - vecpos_off(pos.w, lun.w, jpl->cem / (1.0 + jpl->cem)); - } - break; - case ASSIST_BODY_MARS: - assist_jpl_work(&z[jpl->off[JPL_MAR]], jpl->ncm[JPL_MAR], jpl->ncf[JPL_MAR], jpl->niv[JPL_MAR], t, jpl->inc, pos.u, pos.v, pos.w); - break; - case ASSIST_BODY_JUPITER: - assist_jpl_work(&z[jpl->off[JPL_JUP]], jpl->ncm[JPL_JUP], jpl->ncf[JPL_JUP], jpl->niv[JPL_JUP], t, jpl->inc, pos.u, pos.v, pos.w); - break; - case ASSIST_BODY_SATURN: - assist_jpl_work(&z[jpl->off[JPL_SAT]], jpl->ncm[JPL_SAT], jpl->ncf[JPL_SAT], jpl->niv[JPL_SAT], t, jpl->inc, pos.u, pos.v, pos.w); - break; - case ASSIST_BODY_URANUS: - assist_jpl_work(&z[jpl->off[JPL_URA]], jpl->ncm[JPL_URA], jpl->ncf[JPL_URA], jpl->niv[JPL_URA], t, jpl->inc, pos.u, pos.v, pos.w); - break; - case ASSIST_BODY_NEPTUNE: - assist_jpl_work(&z[jpl->off[JPL_NEP]], jpl->ncm[JPL_NEP], jpl->ncf[JPL_NEP], jpl->niv[JPL_NEP], t, jpl->inc, pos.u, pos.v, pos.w); - break; - case ASSIST_BODY_PLUTO: - assist_jpl_work(&z[jpl->off[JPL_PLU]], jpl->ncm[JPL_PLU], jpl->ncf[JPL_PLU], jpl->niv[JPL_PLU], t, jpl->inc, pos.u, pos.v, pos.w); - break; - default: - return ASSIST_ERROR_NEPHEM; // body not found - break; - } - - // Convert to au/day and au/day^2 - vecpos_div(pos.u, jpl->cau); - vecpos_div(pos.v, jpl->cau/86400.); - vecpos_div(pos.w, jpl->cau/(86400.*86400.)); - - *out_x = pos.u[0]; - *out_y = pos.u[1]; - *out_z = pos.u[2]; - *out_vx = pos.v[0]; - *out_vy = pos.v[1]; - *out_vz = pos.v[2]; - *out_ax = pos.w[0]; - *out_ay = pos.w[1]; - *out_az = pos.w[2]; - - return(ASSIST_SUCCESS); - -} - diff --git a/src/planets.h b/src/planets.h deleted file mode 100644 index a33c5b8..0000000 --- a/src/planets.h +++ /dev/null @@ -1,74 +0,0 @@ -#ifndef _JPL_EPHEM_H -#define _JPL_EPHEM_H -#include "assist.h" - -struct jpl_s * assist_jpl_init(char *str); -void assist_jpl_free(struct jpl_s *jpl); -void assist_jpl_work(double *P, int ncm, int ncf, int niv, double t0, double t1, double *u, double *v, double *w); -enum ASSIST_STATUS assist_jpl_calc(struct jpl_s *pl, double jd_ref, double jd_rel, int body, - double* const GM, - double* const x, double* const y, double* const z, - double* const vx, double* const vy, double* const vz, - double* const ax, double* const ay, double* const az); - -// Order of columns in JPL Ephemeris file -enum { - JPL_MER, // Mercury - JPL_VEN, // Venus - JPL_EMB, // Earth - JPL_MAR, // Mars - JPL_JUP, // Jupiter - JPL_SAT, // Saturn - JPL_URA, // Uranus - JPL_NEP, // Neptune - JPL_PLU, // Pluto - JPL_LUN, // Moon (geocentric) - JPL_SUN, // the Sun - JPL_NUT, // nutations - JPL_LIB, // lunar librations - JPL_MAN, // lunar mantle - JPL_TDB, // TT-TDB (< 2 ms) - - JPL_N, // Number of columns - }; - -struct jpl_s { - double beg, end; // begin and end times - double inc; // time step size - double cau; // definition of AU - double cem; // Earth/Moon mass ratio - int32_t num; // number of constants - int32_t ver; // ephemeris version - int32_t off[JPL_N]; // indexing offset - int32_t ncf[JPL_N]; // number of chebyshev coefficients - int32_t niv[JPL_N]; // number of interpolation intervals - int32_t ncm[JPL_N]; // number of components / dimension - double mass[JPL_N]; // G*mass for all bodies - double J2E; // Other constant names follow JPL - double J3E; - double J4E; - double J2SUN; - double AU; - double RE; - double CLIGHT; - double ASUN; - size_t len, rec; // file and record sizes - void *map; // memory mapped location - double *con; // constant values - char **str; // constant names -}; - -// From Weryk's code -/////// private interface : - -static inline void vecpos_off(double *u, const double *v, const double w) - { u[0] += v[0] * w; u[1] += v[1] * w; u[2] += v[2] * w; } -static inline void vecpos_set(double *u, const double *v) - { u[0] = v[0]; u[1] = v[1]; u[2] = v[2]; } -static inline void vecpos_nul(double *u) - { u[0] = u[1] = u[2] = 0.0; } -static inline void vecpos_div(double *u, double v) - { u[0] /= v; u[1] /= v; u[2] /= v; } - -#endif // _JPL_EPHEM_H - diff --git a/src/spk.c b/src/spk.c index c415ad8..e2fddfd 100644 --- a/src/spk.c +++ b/src/spk.c @@ -1,4 +1,3 @@ - // spk.c - code to handle spice kernel position files // https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/daf.html @@ -9,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -20,8 +20,6 @@ #include "spk.h" - - /* * spk_free * @@ -47,6 +45,85 @@ int assist_spk_free(struct spk_s *pl) } + +void parse_comments(int fd, int first_summary_record, char **comments) { + // Calculate the number of records in the comment section + int num_records = first_summary_record - 2; // Number of comment records + if (num_records <= 0) { + *comments = NULL; + return; + } + + // Allocate a buffer for the comments + int comment_length = num_records * RECORD_LENGTH; + char *buffer = malloc(comment_length + 1); // +1 for null-terminator + if (!buffer) { + perror("malloc"); + exit(EXIT_FAILURE); + } + + // Read the comment records + ssize_t bytes_read; + size_t total_bytes_read = 0; + // Seek to beginning of file + if (lseek(fd, 0, SEEK_SET) == -1) { + perror("lseek"); + free(buffer); + exit(EXIT_FAILURE); + } + // Read in each comment record + for (int i = 0; i < num_records; i++) { + if (lseek(fd, (1 + i) * RECORD_LENGTH, SEEK_SET) == -1) { + perror("lseek"); + free(buffer); + exit(EXIT_FAILURE); + } + bytes_read = read(fd, buffer + total_bytes_read, RECORD_LENGTH); + if (bytes_read == -1) { + perror("read"); + free(buffer); + exit(EXIT_FAILURE); + } + + // Remove all null character and end of comment markers from the end of the buffer + // by deleting the bytes. These pad the ends and create extra newlines. + while (buffer[total_bytes_read + bytes_read - 1] == '\0' || buffer[total_bytes_read + bytes_read - 1] == '\4') { + bytes_read--; + } + + total_bytes_read += bytes_read; + } + + + // DAF comments use the null character to indicate end of a line + // replace with newline character + for (char *p = buffer; p < buffer + total_bytes_read; p++) { + if (*p == '\0') { + *p = '\n'; + } + } + + // Assign the buffer to the comments pointer + *comments = buffer; + + +} + +// sscanf doesn't know the 'D' scientific notation, so replace +// it with 'E' when surrounded by digit and sign characters +void replace_d_with_e(char *str) { + for (int i = 0; str[i] != '\0'; i++) { + if (str[i] == 'D' || str[i] == 'd') { + str[i] = 'e'; + } + } +} + + + + + + /* * assist_spk_init * @@ -56,74 +133,155 @@ int assist_spk_free(struct spk_s *pl) static double inline _jul(double eph) { return 2451545.0 + eph / 86400.0; } +// Populate mass data for spk targets from mass data structure +void assist_spk_join_masses(struct spk_s *sp, const struct mass_data* masses, double emrat) { + if (sp == NULL || masses == NULL || masses->names == NULL) { + return; + } -#define record_length 1024 + // Create an array based mapping of the GMX and target code formats + struct { + const char *name; + int code; + } planet_codes[] = { + {"GMS", 10}, // Sun + {"GM1", 1}, // Mercury + {"GM2", 2}, // Venus + {"GMB", 399}, // Earth + {"GMB", 3}, // Earth-Moon Barycenter + {"GMB", 301}, // Moon + {"GM4", 4}, // Mars + {"GM5", 5}, // Jupiter + {"GM6", 6}, // Saturn + {"GM7", 7}, // Uranus + {"GM8", 8}, // Neptune + {"GM9", 9} // Pluto + }; -struct spk_s * assist_spk_init(const char *path) { - // For file format information, see https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/daf.html + // Join the mass data by iterating through the targets + for (int m = 0; m < sp->num; m++) { + // Skip if mass is already set + if (sp->targets[m].mass != 0) { + continue; + } - // Format for one summary - struct sum_s { - double beg; // begin epoch, seconds since J2000.0 - double end; // ending epoch - int tar; // target code - int cen; // centre code (10 = sun) - int ref; // reference frame (1 = J2000.0) - int ver; // type of ephemeris (2 = chebyshev) - int one; // initial array address - int two; // final array address - }; + // Determine the label we are matching for in the masses array + // If it is a planet, use the lookup value from planet_codes + // If it is an asteroid, we format it as "MA" + code + // Allocate the label and make sure it starts as a null string + char *mass_label = calloc(64, sizeof(char)); + + for (int i = 0; i < sizeof(planet_codes) / sizeof(planet_codes[0]); i++) { + if (sp->targets[m].code == planet_codes[i].code) { + strncpy(mass_label, planet_codes[i].name, 63); + mass_label[63] = '\0'; // Ensure null termination + break; + } + } + // If mass label is still empty, it is an asteroid + if (strlen(mass_label) == 0) { + // Format the asteroid code to be MAdddd where dddd is 4 digit + // 0 masked target code - 2000000 + sprintf(mass_label, "MA%04d", sp->targets[m].code - 2000000); + } - // File is split into records. We read one record at a time. - union { - char buf[record_length]; - struct { - double next; // The record number of the next summary record in the file. Zero if this is the final summary record. - double prev; // The record number of the previous summary record in the file. Zero if this is the initial summary record. - double nsum; // Number of summaries in this record - struct sum_s s[25]; // Summaries (25 is the maximum) - } summary; // Summary record - // See: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/daf.html#The%20File%20Record - struct { - char locidw[8]; // An identification word - int nd; // The number of double precision components in each array summary. - int ni; // The number of integer components in each array summary. - char locifn[60];// The internal name or description of the array file. - int fward; // The record number of the initial summary record in the file. - int bward; // The record number of the final summary record in the file. - } file; // File record - } record; + // Find the mass in the mass arrays + for (int i = 0; i < masses->count; i++) { + if (strcmp(masses->names[i], mass_label) == 0) { + // Earth and moon mass is stored as one value + // Use the ratio to split it using emrat + if (sp->targets[m].code == 399) { + sp->targets[m].mass = masses->values[i] * (emrat / (1. + emrat)); + } else if (sp->targets[m].code == 301) { + sp->targets[m].mass = masses->values[i] * (1./(1.+emrat)); + } else { + sp->targets[m].mass = masses->values[i]; + } + break; + } + } - // Try opening file. - int fd = open(path, O_RDONLY); - if (fd < 0){ + if (sp->targets[m].mass == 0 && sp->targets[m].code != 199 && sp->targets[m].code != 299) { + printf("Mass not found for target code: %d\n", sp->targets[m].code); + } + + free(mass_label); + } +} + +// Load the file record of an spk file + // For file format information, see https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/daf.html +union record_t * assist_load_spk_file_record(int fd) { + // Allocate memory for the record + union record_t *record = calloc(1, sizeof(union record_t)); + if (!record) { + fprintf(stderr, "Memory allocation failed.\n"); + close(fd); + return NULL; + } + + // Seek to the beginning of the file + if (lseek(fd, 0, SEEK_SET) == -1) { + perror("lseek"); + free(record); + close(fd); + return NULL; + } + + // Read the file record + ssize_t bytesRead = read(fd, record, RECORD_LENGTH); + if (bytesRead != RECORD_LENGTH) { + if (bytesRead == -1) { + perror("read"); + } else { + fprintf(stderr, "Incomplete read. Expected %d bytes, got %zd bytes.\n", RECORD_LENGTH, bytesRead); + } + free(record); + close(fd); return NULL; } - // Read the file record. - read(fd, &record, record_length); // Check if the file is a valid Double Precision Array File - if (strncmp(record.file.locidw, "DAF/SPK", 7) != 0) { + if (strncmp(record->file.locidw, "DAF/SPK", 7) != 0) { fprintf(stderr,"Error parsing DAF/SPK file. Incorrect header.\n"); + free(record); close(fd); return NULL; } - // Check that the size of a summary record is equal to the size of our struct. - int nc = 8 * ( record.file.nd + (record.file.ni + 1) / 2 ); + // Check that the size of a summary record is equal to the size of our struct + int nc = 8 * (record->file.nd + (record->file.ni + 1) / 2); if (nc != sizeof(struct sum_s)) { fprintf(stderr,"Error parsing DAF/SPK file. Wrong size of summary record.\n"); + free(record); close(fd); return NULL; } + return record; +} + + +// Initialize the targets of a single spk file +// Note that target masses will not be populated until assist_spk_join_masses +// is called. +struct spk_s * assist_spk_init(const char *path) { + // Try opening file. + int fd = open(path, O_RDONLY); + if (fd < 0){ + return NULL; + } + + // Load the file record + union record_t * record = assist_load_spk_file_record(fd); + // Seek until the first summary record using the file record's fward pointer. // Record numbers start from 1 not 0 so we subtract 1 to get to the correct record. - lseek(fd, (record.file.fward - 1) * record_length, SEEK_SET); - read(fd, record.buf, record_length); + lseek(fd, (record->file.fward - 1) * RECORD_LENGTH, SEEK_SET); + read(fd, record->buf, RECORD_LENGTH); // We are at the first summary block, validate - if ((int64_t)record.buf[8] != 0) { + if ((int64_t)record->buf[8] != 0) { fprintf(stderr, "Error parsing DAF/SPL file. Cannot find summary block.\n"); close(fd); return NULL; @@ -132,45 +290,55 @@ struct spk_s * assist_spk_init(const char *path) { // okay, let's go struct spk_s* pl = calloc(1, sizeof(struct spk_s)); while (1){ // Loop over records - for (int b = 0; b < (int)record.summary.nsum; b++) { // Loop over summaries - struct sum_s* sum = &record.summary.s[b]; // get current summary + for (int b = 0; b < (int)record->summary.nsum; b++) { // Loop over summaries + struct sum_s* sum = &record->summary.s[b]; // get current summary // Index in our arrays for current target - int m = pl->num - 1; + int m = -1; + + // Check to see if we are adding to an existing target + for (int i = 0; i < pl->num; i++) { + if (pl->targets[i].code == sum->tar) { + m = i; + break; + } + } - // New target? - if (pl->num==0 || sum->tar != pl->targets[m].code) { + // If this is a new target, add it to the list + if (m == -1) { + m = pl->num; if (pl->num <= pl->allocated_num){ pl->allocated_num += 32; // increase space in batches of 32 pl->targets = realloc(pl->targets, pl->allocated_num*sizeof(struct spk_target)); } - m++; + pl->num++; pl->targets[m].code = sum->tar; pl->targets[m].cen = sum->cen; pl->targets[m].beg = _jul(sum->beg); pl->targets[m].res = _jul(sum->end) - pl->targets[m].beg; pl->targets[m].one = calloc(32768, sizeof(int)); pl->targets[m].two = calloc(32768, sizeof(int)); - pl->targets[m].ind = 0; - pl->num++; + pl->targets[m].ind = -1; + // Set default of mass to 0 + pl->targets[m].mass = 0; } // add index for target + pl->targets[m].ind++; pl->targets[m].one[pl->targets[m].ind] = sum->one; pl->targets[m].two[pl->targets[m].ind] = sum->two; pl->targets[m].end = _jul(sum->end); - pl->targets[m].ind++; } // Location of next record - long n = (long)record.summary.next - 1; + long n = (long)record->summary.next - 1; if (n<0){ // this is already the last record. break; }else{ // Find and read next record - lseek(fd, n * record_length, SEEK_SET); - read(fd, record.buf, record_length); + lseek(fd, n * RECORD_LENGTH, SEEK_SET); + read(fd, record->buf, RECORD_LENGTH); } } @@ -201,8 +369,9 @@ struct spk_s * assist_spk_init(const char *path) { } -enum ASSIST_STATUS assist_spk_calc(struct spk_s *pl, double jde, double rel, int m, double* GM, double* out_x, double* out_y, double* out_z) +enum ASSIST_STATUS assist_spk_calc(const struct spk_s *pl, double jde, double rel, int m, double* GM, double* out_x, double* out_y, double* out_z) { + if(pl==NULL){ return(ASSIST_ERROR_AST_FILE); } @@ -211,7 +380,6 @@ enum ASSIST_STATUS assist_spk_calc(struct spk_s *pl, double jde, double rel, int return(ASSIST_ERROR_NAST); } struct spk_target* target = &(pl->targets[m]); - if (jde + rel < target->beg || jde + rel > target->end){ return ASSIST_ERROR_COVERAGE; } @@ -279,3 +447,236 @@ enum ASSIST_STATUS assist_spk_calc(struct spk_s *pl, double jde, double rel, int return ASSIST_SUCCESS; } +struct spk_target* assist_spk_find_target(const struct spk_s *pl, int code) { + for (int i = 0; i < pl->num; i++) { + if (pl->targets[i].code == code) { + return &(pl->targets[i]); + } + } + return NULL; +} + +struct mpos_s assist_spk_target_pos(const struct spk_s *pl, const struct spk_target* target, double jde, double rel) +{ + int n, b, p, P, R; + double T[32]; + double S[32]; + double U[32]; + double *val, z; + struct mpos_s pos = {0}; + + // find location of 'directory' describing the data records + n = (int)((jde + rel - target->beg) / target->res); + val = (double *)pl->map + target->two[n] - 1; + + // record size and number of coefficients per coordinate + R = (int)val[-1]; + P = (R - 2) / 3; // must be < 32 !! + + // pick out the precise record + b = (int)(((jde - _jul(val[-3])) + rel) / (val[-2] / 86400.0)); + val = (double *)pl->map + (target->one[n] - 1) + b * R; + + // scale to interpolation units + z = ((jde - _jul(val[0])) + rel) / (val[1] / 86400.0); + + // Calculate the scaling factor 'c' + double c = 1.0 / val[1]; + + // set up Chebyshev polynomials + T[0] = 1.0; T[1] = z; + S[0] = 0.0; S[1] = 1.0; + U[0] = 0.0; U[1] = 0.0; U[2] = 4.0; + + for (p = 2; p < P; p++) { + T[p] = 2.0 * z * T[p-1] - T[p-2]; + S[p] = 2.0 * z * S[p-1] + 2.0 * T[p-1] - S[p-2]; + } + for (p = 3; p < P; p++) { + U[p] = 2.0 * z * U[p-1] + 4.0 * S[p-1] - U[p-2]; + } + + + for (n = 0; n < 3; n++) { + b = 2 + n * P; + pos.u[n] = pos.v[n] = pos.w[n] = 0.0; + + // sum interpolation stuff + for (p = 0; p < P; p++) { + double coeff = val[b + p]; + pos.u[n] += coeff * T[p]; + pos.v[n] += coeff * S[p] * c; + pos.w[n] += coeff * U[p] * c * c; + } + } + + return pos; +} + + +// Calculate the position and velocity of planets from DE440 SPK file +enum ASSIST_STATUS assist_spk_calc_planets(const struct assist_ephem* ephem, double jde, double rel, int code, double* GM, double* out_x, double* out_y, double* out_z, double* out_vx, double* out_vy, double* out_vz, double* out_ax, double* out_ay, double* out_az) +{ + + struct spk_s* pl = ephem->spk_planets; + + + struct spk_target* target = assist_spk_find_target(pl, code); + + if (target == NULL) { + return(ASSIST_ERROR_NEPHEM); + } + + if (jde + rel < target->beg || jde + rel > target->end){ + return ASSIST_ERROR_COVERAGE; + } + + *GM = target->mass; // Note mass constants defined in DE440/441 ephemeris files. If not found mass of 0 is used. + + struct mpos_s pos = assist_spk_target_pos(pl, target, jde, rel); + + // Earth and Moon must be translated from EMB to SSB frame + if (code == 301 || code == 399) { + struct spk_target* emb = assist_spk_find_target(pl, 3); + struct mpos_s emb_pos = assist_spk_target_pos(pl, emb, jde, rel); + + for (int i = 0; i < 3; i++) { + pos.u[i] += emb_pos.u[i]; + pos.v[i] += emb_pos.v[i]; + pos.w[i] += emb_pos.w[i]; + } + + } + + // Convert to AU and AU/day + const double au = ephem->AU; + const double seconds_per_day = 86400.; + + for (int i = 0; i < 3; i++) { + pos.u[i] /= au; + pos.v[i] /= au / seconds_per_day; + pos.w[i] /= au / (seconds_per_day * seconds_per_day); + } + + *out_x = pos.u[0]; + *out_y = pos.u[1]; + *out_z = pos.u[2]; + *out_vx = pos.v[0]; + *out_vy = pos.v[1]; + *out_vz = pos.v[2]; + *out_ax = pos.w[0]; + *out_ay = pos.w[1]; + *out_az = pos.w[2]; + + return ASSIST_SUCCESS; +} + +// Load both constants and masses from SPK file comments in one pass +struct spk_constants_and_masses assist_load_spk_constants_and_masses(const char *path) { + struct spk_constants_and_masses data = {0}; + + // Try opening file. + int fd = open(path, O_RDONLY); + if (fd < 0){ + return data; + } + + // Load the file record + union record_t * record = assist_load_spk_file_record(fd); + if (!record) { + close(fd); + return data; + } + + char *comments = NULL; + parse_comments(fd, record->file.fward, &comments); + free(record); + + if (comments == NULL) { + close(fd); + return data; + } + + char *line = strtok(comments, "\n"); + char key[64]; + char value_str[64]; + double value; + int in_constants_section = 0; + + while (line) { + if (strstr(line, "Initial conditions and constants used for integration:")) { + in_constants_section = 1; + } + + if (in_constants_section) { + replace_d_with_e(line); + if (sscanf(line, "%63s %63s", key, value_str) == 2) { + // Convert the value string to double + value = strtod(value_str, NULL); + + // Parse constants + if (strcmp(key, "cau") == 0) data.AU = value; + else if (strcmp(key, "EMRAT") == 0) data.EMRAT = value; + else if (strcmp(key, "J2E") == 0) data.J2E = value; + else if (strcmp(key, "J3E") == 0) data.J3E = value; + else if (strcmp(key, "J4E") == 0) data.J4E = value; + else if (strcmp(key, "J2SUN") == 0) data.J2SUN = value; + else if (strcmp(key, "AU") == 0) data.AU = value; + else if (strcmp(key, "RE") == 0) data.RE = value; + else if (strcmp(key, "CLIGHT") == 0) data.CLIGHT = value; + else if (strcmp(key, "ASUN") == 0) data.ASUN = value; + // Parse masses + else if (strncmp(key, "GM", 2) == 0 || strncmp(key, "MA", 2) == 0) { + data.masses.names = realloc(data.masses.names, (data.masses.count + 1) * sizeof(char *)); + data.masses.values = realloc(data.masses.values, (data.masses.count + 1) * sizeof(double)); + data.masses.names[data.masses.count] = strdup(key); + data.masses.values[data.masses.count] = value; + data.masses.count++; + } + } + } + + line = strtok(NULL, "\n"); + } + + free(comments); + close(fd); + return data; +} + +// Apply constants from parsed data to ephem structure +void assist_apply_spk_constants(struct assist_ephem* ephem, const struct spk_constants_and_masses* data) { + if (!ephem || !data) return; + + // Apply base constants + ephem->AU = data->AU; + ephem->EMRAT = data->EMRAT; + ephem->J2E = data->J2E; + ephem->J3E = data->J3E; + ephem->J4E = data->J4E; + ephem->J2SUN = data->J2SUN; + ephem->RE = data->RE; + ephem->CLIGHT = data->CLIGHT; + ephem->ASUN = data->ASUN; + + // Calculate derived constants + ephem->Re_eq = ephem->RE / ephem->AU; // Earth radius in AU + ephem->Rs_eq = ephem->ASUN / ephem->AU; // Sun radius in AU + ephem->c_AU_per_day = (ephem->CLIGHT / ephem->AU) * 86400; // Speed of light in AU/day + ephem->c_squared = ephem->c_AU_per_day * ephem->c_AU_per_day; // c^2 in (AU/day)^2 + ephem->over_c_squared = 1.0 / ephem->c_squared; // 1/c^2 in (day/AU)^2 +} + +// Free constants and masses data structure +void assist_free_spk_constants_and_masses(struct spk_constants_and_masses* data) { + if (data && data->masses.names) { + for (int i = 0; i < data->masses.count; i++) { + free(data->masses.names[i]); + } + free(data->masses.names); + free(data->masses.values); + data->masses.names = NULL; + data->masses.values = NULL; + data->masses.count = 0; + } +} diff --git a/src/spk.h b/src/spk.h index 4de1a5f..0d0df4e 100644 --- a/src/spk.h +++ b/src/spk.h @@ -1,4 +1,3 @@ - // spk.h // https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/daf.html @@ -10,6 +9,9 @@ #include "assist.h" +#define RECORD_LENGTH 1024 +#define MAX_COMMENT_LENGTH 1000000 + struct mpos_s { double u[3]; double v[3]; @@ -17,7 +19,29 @@ struct mpos_s { double jde; }; +struct mass_data { + char **names; + double *values; + int count; +}; + +// Combined structure for constants and masses from SPK file +struct spk_constants_and_masses { + // Constants + double AU; + double EMRAT; + double J2E; + double J3E; + double J4E; + double J2SUN; + double RE; + double CLIGHT; + double ASUN; + // Mass data + struct mass_data masses; +}; +// Represents available data for a target in a SPK file struct spk_target { int code; // Target code int cen; // Centre target @@ -31,6 +55,7 @@ struct spk_target { }; +// Represents a collection of targets in an SPK file struct spk_s { struct spk_target* targets; int num; // number of targets @@ -39,10 +64,55 @@ struct spk_s { size_t len; // map length }; + // Format for one summary +struct sum_s { + double beg; // begin epoch, seconds since J2000.0 + double end; // ending epoch + int tar; // target code + int cen; // centre code (10 = sun) + int ref; // reference frame (1 = J2000.0) + int ver; // type of ephemeris (2 = chebyshev) + int one; // initial array address + int two; // final array address +}; + +// File is split into records. We read one record at a time. +union record_t { + char buf[RECORD_LENGTH]; + struct { + double next; // The record number of the next summary record in the file. Zero if this is the final summary record. + double prev; // The record number of the previous summary record in the file. Zero if this is the initial summary record. + double nsum; // Number of summaries in this record + struct sum_s s[25]; // Summaries (25 is the maximum) + } summary; // Summary record + // See: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/daf.html#The%20File%20Record + struct { + char locidw[8]; // An identification word + int nd; // The number of double precision components in each array summary. + int ni; // The number of integer components in each array summary. + char locifn[60];// The internal name or description of the array file. + int fward; // The record number of the initial summary record in the file. + int bward; // The record number of the final summary record in the file. + int free; // Next available free record + } file; // File record +}; + + + + int assist_spk_free(struct spk_s *pl); struct spk_s * assist_spk_init(const char *path); -enum ASSIST_STATUS assist_spk_calc(struct spk_s *pl, double jde, double rel, int m, double* GM, double* x, double* y, double* z); +union record_t * assist_load_spk_file_record(int fd); +struct spk_constants_and_masses assist_load_spk_constants_and_masses(const char *path); +void assist_apply_spk_constants(struct assist_ephem* ephem, const struct spk_constants_and_masses* data); +void assist_free_spk_constants_and_masses(struct spk_constants_and_masses* data); +void parse_comments(int fd, int first_summary_record, char **comments); +void assist_spk_join_masses(struct spk_s *sp, const struct mass_data* masses, double emrat); +enum ASSIST_STATUS assist_spk_calc(const struct spk_s *pl, double jde, double rel, int m, double* GM, double* x, double* y, double* z); +enum ASSIST_STATUS assist_spk_calc_planets(const struct assist_ephem* ephem, double jde, double rel, int m, double* GM, double* x, double* y, double* z, double* vx, double* vy, double* vz, double* ax, double* ay, double* az); +struct spk_target* assist_spk_find_target(const struct spk_s *pl, int code); +struct mpos_s assist_spk_target_pos(const struct spk_s *pl, const struct spk_target* target, double jde, double rel); #endif // _SPK_H diff --git a/src/tools.c b/src/tools.c index d5db407..bc8c6d9 100644 --- a/src/tools.c +++ b/src/tools.c @@ -33,7 +33,7 @@ // This function returns a new rebound simulation with all particles added. // Return simulation must be freed by caller. // if merge_moon = 1, then the Earth and Moon are added as a single particle. -struct reb_simulation* assist_simulation_convert_to_rebound(struct reb_simulation* r, struct assist_ephem* ephem, int merge_moon){ +struct reb_simulation* assist_simulation_convert_to_rebound(const struct reb_simulation* r, const struct assist_ephem* ephem, int merge_moon){ struct reb_simulation* r2 = reb_simulation_create(); r2->t = r->t; r2->dt = r->dt; diff --git a/unit_tests/5303_Ceres/problem.c b/unit_tests/5303_Ceres/problem.c index 1497b92..3f5b480 100644 --- a/unit_tests/5303_Ceres/problem.c +++ b/unit_tests/5303_Ceres/problem.c @@ -10,7 +10,7 @@ int main(int argc, char* argv[]){ struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); diff --git a/unit_tests/apophis/problem.c b/unit_tests/apophis/problem.c index 94d7b69..3e798de 100644 --- a/unit_tests/apophis/problem.c +++ b/unit_tests/apophis/problem.c @@ -11,7 +11,7 @@ int main(int argc, char* argv[]){ struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); @@ -70,8 +70,8 @@ int main(int argc, char* argv[]){ double diff = reb_particle_distance(&p_final, &r->particles[0]) * au2meter; // in meter // Ensure accuracy is better than 250m - printf("Difference: %.2f m\n",diff); - assert(diff < 250); + fprintf(stderr,"diff to JPL Small Body code %fm\n",diff); + assert(diff < 200); assist_free(ax); assist_ephem_free(ephem); diff --git a/unit_tests/benchmark/problem.c b/unit_tests/benchmark/problem.c index 44cad88..ec208b0 100644 --- a/unit_tests/benchmark/problem.c +++ b/unit_tests/benchmark/problem.c @@ -10,7 +10,7 @@ double runtime(){ struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); diff --git a/unit_tests/convert_simulation/problem.c b/unit_tests/convert_simulation/problem.c index 3479eec..ba0c85e 100644 --- a/unit_tests/convert_simulation/problem.c +++ b/unit_tests/convert_simulation/problem.c @@ -12,7 +12,7 @@ const double au2meter = 149597870700; int main(int argc, char* argv[]){ struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); diff --git a/unit_tests/ephem_cache/problem.c b/unit_tests/ephem_cache/problem.c index 1f55c33..933001d 100644 --- a/unit_tests/ephem_cache/problem.c +++ b/unit_tests/ephem_cache/problem.c @@ -34,7 +34,7 @@ struct reb_particle integrate(struct assist_ephem* ephem, int cache_on, int dire int main(int argc, char* argv[]){ struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); diff --git a/unit_tests/holman/problem.c b/unit_tests/holman/problem.c index 4c5957d..f01ba33 100644 --- a/unit_tests/holman/problem.c +++ b/unit_tests/holman/problem.c @@ -31,7 +31,7 @@ int main(int argc, char* argv[]){ struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); diff --git a/unit_tests/holman_reverse/problem.c b/unit_tests/holman_reverse/problem.c index d8b3a10..799e99d 100644 --- a/unit_tests/holman_reverse/problem.c +++ b/unit_tests/holman_reverse/problem.c @@ -10,7 +10,7 @@ int main(int argc, char* argv[]){ struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); diff --git a/unit_tests/interpolation/problem.c b/unit_tests/interpolation/problem.c index 722b76a..052ceff 100644 --- a/unit_tests/interpolation/problem.c +++ b/unit_tests/interpolation/problem.c @@ -13,7 +13,7 @@ int main(int argc, char* argv[]){ remove("out.bin"); struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); diff --git a/unit_tests/onthefly_backwards_interpolation/problem.c b/unit_tests/onthefly_backwards_interpolation/problem.c index ddbec7d..462d3c7 100644 --- a/unit_tests/onthefly_backwards_interpolation/problem.c +++ b/unit_tests/onthefly_backwards_interpolation/problem.c @@ -6,7 +6,7 @@ int main(int argc, char* argv[]){ struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); diff --git a/unit_tests/onthefly_interpolation/problem.c b/unit_tests/onthefly_interpolation/problem.c index 86b1135..076e042 100644 --- a/unit_tests/onthefly_interpolation/problem.c +++ b/unit_tests/onthefly_interpolation/problem.c @@ -6,7 +6,7 @@ int main(int argc, char* argv[]){ struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); @@ -51,7 +51,7 @@ int main(int argc, char* argv[]){ // Compare results for (int i=0; i +#include +#include +#include "spk.h" + +int main(int argc, char *argv[]) { + char *planet_path = "../../data/de440.bsp"; + char *asteroid_path = "../../data/sb441-n16.bsp"; + + struct spk_s *planet_spk = assist_spk_init(planet_path); + if (!planet_spk) { + fprintf(stderr, "Failed to initialize SPK data.\n"); + return EXIT_FAILURE; + } + + // Assert that there are 14 targets in the planet SPK data + // 1 : Mercury + // 2 : Venus + // 3 : Earth-Moon barycenter + // 4 : Mars + // 5 : Jupiter + // 6 : Saturn + // 7 : Uranus + // 8 : Neptune + // 9 : Pluto + // 10 : Sun + // 301 : Moon + // 399 : Earth + // 199 : Mercury (Mercury barycenter) + // 299 : Venus (Venus barycenter) + + assert(planet_spk->num == 14); + int counted_targets = 0; + + for (int i = 0; i < 16; i++) { + // Assuming a valid target has a non-zero code, adjust the condition as necessary + if (planet_spk->targets[i].code != 0) { + counted_targets++; + } + } + assert(counted_targets == 14); + + assist_spk_free(planet_spk); + + struct spk_s *asteroid_spk = assist_spk_init(asteroid_path); + if (!asteroid_spk) { + fprintf(stderr, "Failed to initialize SPK data.\n"); + return EXIT_FAILURE; + } + + // Assert that there are 16 targets in the asteroid SPK data + assert(asteroid_spk->num == 16); + counted_targets = 0; + + for (int i = 0; i < 20; i++) { + // Assuming a valid target has a non-zero code, adjust the condition as necessary + if (asteroid_spk->targets[i].code != 0) { + counted_targets++; + } + } + + assist_spk_free(asteroid_spk); + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/unit_tests/spk_join_masses/Makefile b/unit_tests/spk_join_masses/Makefile new file mode 100644 index 0000000..e9701f7 --- /dev/null +++ b/unit_tests/spk_join_masses/Makefile @@ -0,0 +1,46 @@ +ifndef REB_DIR +ifneq ($(wildcard ../../../rebound/.*),) # Check for REBOUND in default location +REB_DIR=../../../rebound +endif +ifneq ($(wildcard ../../../../rebound/.*),) # Check for ASSIST being inside REBOUND directory +REB_DIR=../../../ +endif +endif +ifndef REB_DIR # REBOUND is not in default location and REB_DIR is not set + $(error ASSIST not in the same directory as REBOUND. To use a custom location, you Must set the REB_DIR environment variable for the path to your rebound directory, e.g., export REB_DIR=/Users/dtamayo/rebound.) +endif +PROBLEMDIR=$(shell basename `dirname \`pwd\``)"/"$(shell basename `pwd`) + +include $(REB_DIR)/src/Makefile.defs + +ASSIST_DIR=../../ + +all: librebound.so libassist.so + @echo "" + @echo "Compiling problem file ..." + $(CC) -I$(ASSIST_DIR)/src/ -I$(REB_DIR)/src/ -Wl,-rpath,./ $(OPT) $(PREDEF) problem.c -L. -lassist -lrebound $(LIB) -o rebound + @echo "" + @echo "Problem file compiled successfully." + +librebound.so: + @echo "Compiling shared library librebound.so ..." + $(MAKE) -C $(REB_DIR)/src/ + @echo "Creating link for shared library librebound.so ..." + @-rm -f librebound.so + @ln -s $(REB_DIR)/src/librebound.so . + +libassist.so: librebound.so $(ASSIST_DIR)/src/*.h $(ASSIST_DIR)/src/*.c + @echo "Compiling shared library libassist.so ..." + $(MAKE) -C $(ASSIST_DIR)/src/ + @-rm -f libassist.so + @ln -s $(ASSIST_DIR)/src/libassist.so . + +clean: + @echo "Cleaning up shared library librebound.so ..." + @-rm -f librebound.so + $(MAKE) -C $(REB_DIR)/src/ clean + @echo "Cleaning up shared library libassist.so ..." + @-rm -f libassist.so + $(MAKE) -C $(ASSIST_DIR)/src/ clean + @echo "Cleaning up local directory ..." + @-rm -vf rebound diff --git a/unit_tests/spk_join_masses/problem.c b/unit_tests/spk_join_masses/problem.c new file mode 100644 index 0000000..31b4625 --- /dev/null +++ b/unit_tests/spk_join_masses/problem.c @@ -0,0 +1,64 @@ +#include +#include +#include +#include "assist.h" +#include "spk.h" + + +void assert_populated_masses(struct spk_s *spk_data) { + if (!spk_data) { + printf("SPK data is NULL.\n"); + return; + } + + // Ensure we have targets loaded + assert(spk_data->num > 0); + printf("Found %d targets in SPK data.\n", spk_data->num); + + for (int i = 0; i < spk_data->num; i++) { + struct spk_target *target = &spk_data->targets[i]; + // We do not calculate mass for mercury and venus barycenters + if (target->code != 199 && target->code != 299) { + // Assert that the masses are not zero + assert(target->mass != 0.0); + } + } +} + + +int main(int argc, char *argv[]) { + + const char *planets_path = "../../data/de440.bsp"; + const char *asteroids_path = "../../data/sb441-n16.bsp"; + + struct assist_ephem ephem = {0}; + + // Load planets + ephem.spk_planets = assist_spk_init(planets_path); + if (!ephem.spk_planets) { + fprintf(stderr, "Failed to initialize planet SPK data.\n"); + return EXIT_FAILURE; + } + + // Load asteroids + ephem.spk_asteroids = assist_spk_init(asteroids_path); + if (!ephem.spk_asteroids) { + fprintf(stderr, "Failed to initialize asteroid SPK data.\n"); + return EXIT_FAILURE; + } + + // Load constants and masses, then apply them + struct spk_constants_and_masses data = assist_load_spk_constants_and_masses(planets_path); + assist_apply_spk_constants(&ephem, &data); + assist_spk_join_masses(ephem.spk_planets, &data.masses, ephem.EMRAT); + assist_spk_join_masses(ephem.spk_asteroids, &data.masses, ephem.EMRAT); + assist_free_spk_constants_and_masses(&data); + + assert_populated_masses(ephem.spk_planets); + assert_populated_masses(ephem.spk_asteroids); + + assist_spk_free(ephem.spk_planets); + assist_spk_free(ephem.spk_asteroids); + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/unit_tests/spk_load_constants/Makefile b/unit_tests/spk_load_constants/Makefile new file mode 100644 index 0000000..e9701f7 --- /dev/null +++ b/unit_tests/spk_load_constants/Makefile @@ -0,0 +1,46 @@ +ifndef REB_DIR +ifneq ($(wildcard ../../../rebound/.*),) # Check for REBOUND in default location +REB_DIR=../../../rebound +endif +ifneq ($(wildcard ../../../../rebound/.*),) # Check for ASSIST being inside REBOUND directory +REB_DIR=../../../ +endif +endif +ifndef REB_DIR # REBOUND is not in default location and REB_DIR is not set + $(error ASSIST not in the same directory as REBOUND. To use a custom location, you Must set the REB_DIR environment variable for the path to your rebound directory, e.g., export REB_DIR=/Users/dtamayo/rebound.) +endif +PROBLEMDIR=$(shell basename `dirname \`pwd\``)"/"$(shell basename `pwd`) + +include $(REB_DIR)/src/Makefile.defs + +ASSIST_DIR=../../ + +all: librebound.so libassist.so + @echo "" + @echo "Compiling problem file ..." + $(CC) -I$(ASSIST_DIR)/src/ -I$(REB_DIR)/src/ -Wl,-rpath,./ $(OPT) $(PREDEF) problem.c -L. -lassist -lrebound $(LIB) -o rebound + @echo "" + @echo "Problem file compiled successfully." + +librebound.so: + @echo "Compiling shared library librebound.so ..." + $(MAKE) -C $(REB_DIR)/src/ + @echo "Creating link for shared library librebound.so ..." + @-rm -f librebound.so + @ln -s $(REB_DIR)/src/librebound.so . + +libassist.so: librebound.so $(ASSIST_DIR)/src/*.h $(ASSIST_DIR)/src/*.c + @echo "Compiling shared library libassist.so ..." + $(MAKE) -C $(ASSIST_DIR)/src/ + @-rm -f libassist.so + @ln -s $(ASSIST_DIR)/src/libassist.so . + +clean: + @echo "Cleaning up shared library librebound.so ..." + @-rm -f librebound.so + $(MAKE) -C $(REB_DIR)/src/ clean + @echo "Cleaning up shared library libassist.so ..." + @-rm -f libassist.so + $(MAKE) -C $(ASSIST_DIR)/src/ clean + @echo "Cleaning up local directory ..." + @-rm -vf rebound diff --git a/unit_tests/spk_load_constants/problem.c b/unit_tests/spk_load_constants/problem.c new file mode 100644 index 0000000..9b9bee2 --- /dev/null +++ b/unit_tests/spk_load_constants/problem.c @@ -0,0 +1,33 @@ +#include +#include +#include +#include "assist.h" +#include "spk.h" + + +void assert_constants_exist(struct assist_ephem *ephem) { + assert(ephem->AU != 0.0); + assert(ephem->EMRAT != 0.0); + assert(ephem->J2E != 0.0); + assert(ephem->J3E != 0.0); + assert(ephem->J4E != 0.0); + assert(ephem->J2SUN != 0.0); + assert(ephem->RE != 0.0); + assert(ephem->CLIGHT != 0.0); + assert(ephem->ASUN != 0.0); + printf("All constants exist.\n"); +} + +int main(int argc, char *argv[]) { + + const char *planets_path = "../../data/de440.bsp"; + + struct assist_ephem ephem = {0}; + struct spk_constants_and_masses data = assist_load_spk_constants_and_masses(planets_path); + assist_apply_spk_constants(&ephem, &data); + assist_free_spk_constants_and_masses(&data); + + assert_constants_exist(&ephem); + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/unit_tests/variational/problem.c b/unit_tests/variational/problem.c index ecaf4bd..813823f 100644 --- a/unit_tests/variational/problem.c +++ b/unit_tests/variational/problem.c @@ -12,7 +12,7 @@ const double au2meter = 149597870700; int main(int argc, char* argv[]){ struct assist_ephem* ephem = assist_ephem_create( - "../../data/linux_p1550p2650.440", + "../../data/de440.bsp", "../../data/sb441-n16.bsp"); if (ephem == NULL){ fprintf(stderr,"Error initializing assist_ephem.\n"); From d832ef36866f19b059bc03e30877ec22e6b85083 Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Tue, 27 May 2025 14:10:00 -0400 Subject: [PATCH 2/6] Add note to run all tests --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 5ee8242..9ea076c 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ libassist: $(MAKE) -C src @ln -f -s src/libassist.so . - @if [ "$(MAKELEVEL)" -eq "0" ]; then echo "To compile the example problems, go to a subdirectory of examples/ and execute make there."; fi + @if [ "$(MAKELEVEL)" -eq "0" ]; then echo "To run tests, run 'make test' or go to a subdirectory of examples/ and execute make there."; fi all: libassist From 31e318807f5665740f55a2ff9e879507cf2345e0 Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Tue, 24 Jun 2025 10:05:48 -0400 Subject: [PATCH 3/6] check file type --- src/spk.c | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/spk.c b/src/spk.c index e2fddfd..77f4616 100644 --- a/src/spk.c +++ b/src/spk.c @@ -261,6 +261,28 @@ union record_t * assist_load_spk_file_record(int fd) { return record; } +// Check if file is an ASCII source file that needs conversion +static int is_ascii_source_file(int fd) { + char buf[1024]; + ssize_t bytes_read; + + // Read first chunk of file + lseek(fd, 0, SEEK_SET); + bytes_read = read(fd, buf, sizeof(buf)); + if (bytes_read <= 0) { + return 0; + } + + // Check for ASCII source file markers + // These files typically contain header information about the ephemeris + // and start with comment blocks + if (strstr(buf, "JPL") && strstr(buf, "EPHEMERIS") && + (strstr(buf, "ASCII") || strstr(buf, "ascii"))) { + return 1; + } + + return 0; +} // Initialize the targets of a single spk file // Note that target masses will not be populated until assist_spk_join_masses @@ -272,8 +294,23 @@ struct spk_s * assist_spk_init(const char *path) { return NULL; } + // Check for ASCII source file + if (is_ascii_source_file(fd)) { + fprintf(stderr, "(ASSIST) Error: ASCII source file format detected. This file needs to be converted to binary SPK (.bsp) format.\n"); + fprintf(stderr, "(ASSIST) Please download the pre-converted binary SPK files from:\n"); + fprintf(stderr, "(ASSIST) https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/\n"); + fprintf(stderr, "(ASSIST) and https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/asteroids/\n"); + fprintf(stderr, "(ASSIST) Look for files with .bsp extension (e.g., de440.bsp, sb441-n16.bsp)\n"); + close(fd); + return NULL; + } + // Load the file record union record_t * record = assist_load_spk_file_record(fd); + if (!record) { + close(fd); + return NULL; + } // Seek until the first summary record using the file record's fward pointer. // Record numbers start from 1 not 0 so we subtract 1 to get to the correct record. From 011dd44161d7629ca402bce36350578fcee0c662 Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Tue, 24 Jun 2025 10:22:16 -0400 Subject: [PATCH 4/6] Try with loader_path and fix --- setup.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/setup.py b/setup.py index bc58fa6..0b06ca3 100644 --- a/setup.py +++ b/setup.py @@ -53,6 +53,13 @@ def finalize_options(self): ext.extra_link_args.append('-Wl,-rpath,'+rebdir+'/../') ext.runtime_library_dirs.append(rebdirsp) ext.extra_link_args.append('-Wl,-rpath,'+rebdirsp) + + # Add loader-relative paths for virtual environments where site-packages gets moved + if sys.platform == 'darwin': + ext.extra_link_args.append('-Wl,-rpath,@loader_path') + elif sys.platform.startswith('linux'): + ext.extra_link_args.append('-Wl,-rpath,$ORIGIN') + print(extra_link_args) print(rebdir+'/../') print(rebdirsp) @@ -72,6 +79,11 @@ def finalize_options(self): vars = sysconfig.get_config_vars() vars['LDSHARED'] = vars['LDSHARED'].replace('-bundle', '-shared') extra_link_args.append('-Wl,-install_name,@rpath/libassist'+suffix) + # Add additional linker flags for better compatibility with virtual environments + extra_link_args.append('-Wl,-rpath,@loader_path') +elif sys.platform.startswith('linux'): + # Add Linux equivalent for virtual environment compatibility + extra_link_args.append('-Wl,-rpath,$ORIGIN') libassistmodule = Extension('libassist', sources = [ 'src/assist.c','src/spk.c', 'src/forces.c', 'src/tools.c'], From 3ea74e8431db940c6b9c8724ceba1f0c18c69bfe Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Tue, 1 Jul 2025 09:59:39 -0400 Subject: [PATCH 5/6] try different --- setup.py | 32 ++++++++++--------- .../benchmark/benchmark_MacBook-Pro.local.txt | 2 ++ 2 files changed, 19 insertions(+), 15 deletions(-) create mode 100644 unit_tests/benchmark/benchmark_MacBook-Pro.local.txt diff --git a/setup.py b/setup.py index 0b06ca3..8aca202 100644 --- a/setup.py +++ b/setup.py @@ -46,30 +46,34 @@ def finalize_options(self): sources = [ 'src/assist.c', 'src/spk.c', 'src/forces.c', 'src/tools.c'], if not "CONDA_BUILD_CROSS_COMPILATION" in os.environ: + # Add library directories for build-time linking self.library_dirs.append(rebdir+'/../') self.library_dirs.append(rebdirsp) + for ext in self.extensions: - ext.runtime_library_dirs.append(rebdir+'/../') - ext.extra_link_args.append('-Wl,-rpath,'+rebdir+'/../') - ext.runtime_library_dirs.append(rebdirsp) - ext.extra_link_args.append('-Wl,-rpath,'+rebdirsp) + # Clear any existing runtime_library_dirs to avoid conflicts + ext.runtime_library_dirs = [] - # Add loader-relative paths for virtual environments where site-packages gets moved + # Use loader-relative paths for runtime linking if sys.platform == 'darwin': + # Primary path: same directory as the loading library ext.extra_link_args.append('-Wl,-rpath,@loader_path') + # Backup path: current directory (for compatibility) + ext.extra_link_args.append('-Wl,-rpath,@loader_path/.') elif sys.platform.startswith('linux'): + # Linux equivalent of @loader_path ext.extra_link_args.append('-Wl,-rpath,$ORIGIN') + ext.extra_link_args.append('-Wl,-rpath,$ORIGIN/.') - print(extra_link_args) - print(rebdir+'/../') - print(rebdirsp) + print("Library directories:", self.library_dirs) + print("Extra link args:", ext.extra_link_args) else: # For conda-forge cross-compile builds rebdir=get_python_lib(prefix=os.environ["PREFIX"]) self.library_dirs.append(rebdir) for ext in self.extensions: ext.extra_link_args.append('-Wl,-rpath,'+rebdir) - print(ext.extra_link_args) + print("Cross-compile extra link args:", ext.extra_link_args) from distutils.version import LooseVersion @@ -79,17 +83,15 @@ def finalize_options(self): vars = sysconfig.get_config_vars() vars['LDSHARED'] = vars['LDSHARED'].replace('-bundle', '-shared') extra_link_args.append('-Wl,-install_name,@rpath/libassist'+suffix) - # Add additional linker flags for better compatibility with virtual environments - extra_link_args.append('-Wl,-rpath,@loader_path') elif sys.platform.startswith('linux'): - # Add Linux equivalent for virtual environment compatibility - extra_link_args.append('-Wl,-rpath,$ORIGIN') + # Linux doesn't need install_name, but we can add other flags if needed + pass libassistmodule = Extension('libassist', sources = [ 'src/assist.c','src/spk.c', 'src/forces.c', 'src/tools.c'], include_dirs = ['src'], library_dirs = [], - runtime_library_dirs = ["."], + runtime_library_dirs = [], # Will be set by build_ext.finalize_options libraries=['rebound'+suffix[:suffix.rfind('.')]], define_macros=[ ('LIBASSIST', None) ], extra_compile_args=['-fstrict-aliasing', '-O3','-std=c99', '-fPIC', '-D_GNU_SOURCE', '-Wpointer-arith', ghash_arg], @@ -138,4 +140,4 @@ def finalize_options(self): tests_require=["numpy","matplotlib","rebound"], test_suite="assist.test", ext_modules = [libassistmodule], - zip_safe=False) + zip_safe=False) \ No newline at end of file diff --git a/unit_tests/benchmark/benchmark_MacBook-Pro.local.txt b/unit_tests/benchmark/benchmark_MacBook-Pro.local.txt new file mode 100644 index 0000000..fbc9114 --- /dev/null +++ b/unit_tests/benchmark/benchmark_MacBook-Pro.local.txt @@ -0,0 +1,2 @@ +0.021243s 0.006105s 1.1.9 May 27 2025 13:55:45 575dc570608a93deb448144d66ede3176576b33a +0.020310s 0.002809s 1.1.9 May 27 2025 13:55:45 575dc570608a93deb448144d66ede3176576b33a From b1efd344b12d4eaaad05e00b77ba0d1803bc5437 Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Wed, 2 Jul 2025 09:42:23 -0400 Subject: [PATCH 6/6] Now detect legacy ephemeris files --- src/spk.c | 100 +++++++++++++++++--- src/spk.h | 12 +++ unit_tests/spk_detection/Makefile | 54 +++++++++++ unit_tests/spk_detection/problem.c | 142 +++++++++++++++++++++++++++++ 4 files changed, 296 insertions(+), 12 deletions(-) create mode 100644 unit_tests/spk_detection/Makefile create mode 100644 unit_tests/spk_detection/problem.c diff --git a/src/spk.c b/src/spk.c index 77f4616..5147e9a 100644 --- a/src/spk.c +++ b/src/spk.c @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -261,8 +262,57 @@ union record_t * assist_load_spk_file_record(int fd) { return record; } -// Check if file is an ASCII source file that needs conversion -static int is_ascii_source_file(int fd) { +// Detect legacy ephemeris file formats and provide helpful error messages +// (enum ephemeris_file_format_t is now defined in spk.h) + +// Simple signature check for legacy binary format +int assist_detect_legacy_binary_signature(int fd) { + // Legacy binary files have constant names at offset 0x00FC (252 bytes) + // Check for typical ephemeris constant names like "GMS ", "GM1 ", etc. + + char const_names[6 * 3]; // Read first 3 constant names (6 chars each) + if (lseek(fd, 0x00FC, SEEK_SET) != 0x00FC) { + return 0; // Can't seek to expected position + } + + if (read(fd, const_names, sizeof(const_names)) != sizeof(const_names)) { + return 0; // Can't read constant names + } + + // Check if names look like typical ephemeris constants + // They should be printable ASCII with trailing spaces, like "GMS ", "GM1 " + int valid_names = 0; + for (int i = 0; i < 3; i++) { + char* name = &const_names[i * 6]; + + // Check if it looks like a constant name: starts with letters, may have trailing spaces + if ((name[0] >= 'A' && name[0] <= 'Z') || (name[0] >= 'a' && name[0] <= 'z')) { + int has_letters = 0; + int has_spaces = 0; + for (int j = 0; j < 6; j++) { + if ((name[j] >= 'A' && name[j] <= 'Z') || (name[j] >= 'a' && name[j] <= 'z') || + (name[j] >= '0' && name[j] <= '9')) { + has_letters++; + } else if (name[j] == ' ') { + has_spaces++; + } else if (name[j] == '\0') { + // Null termination is ok + break; + } + } + + // Valid constant name should have at least 2 letters and some spaces for padding + if (has_letters >= 2 && (has_spaces > 0 || strlen(name) < 6)) { + valid_names++; + } + } + } + + // If at least 2 out of 3 names look like ephemeris constants, it's likely a legacy file + return (valid_names >= 2); +} + +ephemeris_file_format_t assist_detect_ephemeris_file_format(int fd) { char buf[1024]; ssize_t bytes_read; @@ -270,7 +320,19 @@ static int is_ascii_source_file(int fd) { lseek(fd, 0, SEEK_SET); bytes_read = read(fd, buf, sizeof(buf)); if (bytes_read <= 0) { - return 0; + return FILE_FORMAT_UNKNOWN; + } + + // Check for valid .bsp SPK file (DAF/SPK header) first + if (bytes_read >= 8 && strncmp(buf, "DAF/SPK ", 8) == 0) { + return FILE_FORMAT_VALID_BSP; + } + + // Check for legacy binary format signature + if (assist_detect_legacy_binary_signature(fd)) { + // Reset file position after signature check + lseek(fd, 0, SEEK_SET); + return FILE_FORMAT_BINARY_LEGACY; } // Check for ASCII source file markers @@ -278,10 +340,27 @@ static int is_ascii_source_file(int fd) { // and start with comment blocks if (strstr(buf, "JPL") && strstr(buf, "EPHEMERIS") && (strstr(buf, "ASCII") || strstr(buf, "ascii"))) { - return 1; + return FILE_FORMAT_ASCII_SOURCE; } - return 0; + return FILE_FORMAT_UNKNOWN; +} + +static void print_legacy_format_error(ephemeris_file_format_t format, const char* filepath) { + if (format == FILE_FORMAT_ASCII_SOURCE) { + fprintf(stderr, "(ASSIST) Error: ASCII ephemeris source file detected. ASSIST requires binary SPK (.bsp) format.\n"); + } else if (format == FILE_FORMAT_BINARY_LEGACY) { + fprintf(stderr, "(ASSIST) Error: Legacy binary ephemeris file detected. ASSIST requires binary SPK (.bsp) format.\n"); + } else { + fprintf(stderr, "(ASSIST) Error: Unsupported ephemeris file format. ASSIST requires binary SPK (.bsp) format.\n"); + } + + fprintf(stderr, "(ASSIST) Download from: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440.bsp\n"); + + // Check if we're running in Python environment + if (getenv("PYTHONPATH") || getenv("VIRTUAL_ENV") || getenv("CONDA_DEFAULT_ENV")) { + fprintf(stderr, "(ASSIST) Or install via: pip install naif-de440 jpl-small-bodies-de441-n16\n"); + } } // Initialize the targets of a single spk file @@ -294,13 +373,10 @@ struct spk_s * assist_spk_init(const char *path) { return NULL; } - // Check for ASCII source file - if (is_ascii_source_file(fd)) { - fprintf(stderr, "(ASSIST) Error: ASCII source file format detected. This file needs to be converted to binary SPK (.bsp) format.\n"); - fprintf(stderr, "(ASSIST) Please download the pre-converted binary SPK files from:\n"); - fprintf(stderr, "(ASSIST) https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/\n"); - fprintf(stderr, "(ASSIST) and https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/asteroids/\n"); - fprintf(stderr, "(ASSIST) Look for files with .bsp extension (e.g., de440.bsp, sb441-n16.bsp)\n"); + // Check for legacy ephemeris file formats + ephemeris_file_format_t format = assist_detect_ephemeris_file_format(fd); + if (format != FILE_FORMAT_VALID_BSP) { + print_legacy_format_error(format, path); close(fd); return NULL; } diff --git a/src/spk.h b/src/spk.h index 0d0df4e..10262d7 100644 --- a/src/spk.h +++ b/src/spk.h @@ -41,6 +41,14 @@ struct spk_constants_and_masses { struct mass_data masses; }; +// Ephemeris file format detection +typedef enum { + FILE_FORMAT_VALID_BSP = 0, // Valid .bsp SPK file + FILE_FORMAT_ASCII_SOURCE = 1, // ASCII source file (needs conversion) + FILE_FORMAT_BINARY_LEGACY = 2, // Binary legacy format (.440, .441, etc.) + FILE_FORMAT_UNKNOWN = 3 // Unknown/unsupported format +} ephemeris_file_format_t; + // Represents available data for a target in a SPK file struct spk_target { int code; // Target code @@ -114,5 +122,9 @@ enum ASSIST_STATUS assist_spk_calc_planets(const struct assist_ephem* ephem, dou struct spk_target* assist_spk_find_target(const struct spk_s *pl, int code); struct mpos_s assist_spk_target_pos(const struct spk_s *pl, const struct spk_target* target, double jde, double rel); +// File format detection functions +int assist_detect_legacy_binary_signature(int fd); +ephemeris_file_format_t assist_detect_ephemeris_file_format(int fd); + #endif // _SPK_H diff --git a/unit_tests/spk_detection/Makefile b/unit_tests/spk_detection/Makefile new file mode 100644 index 0000000..dace851 --- /dev/null +++ b/unit_tests/spk_detection/Makefile @@ -0,0 +1,54 @@ +ifndef REB_DIR +ifneq ($(wildcard ../../../rebound/.*),) # Check for REBOUND in default location +REB_DIR=../../../rebound +endif +ifneq ($(wildcard ../../../../rebound/.*),) # Check for ASSIST being inside REBOUND directory +REB_DIR=../../../ +endif +endif +ifndef REB_DIR # REBOUND is not in default location and REB_DIR is not set + $(error ASSIST not in the same directory as REBOUND. To use a custom location, you Must set the REB_DIR environment variable for the path to your rebound directory, e.g., export REB_DIR=/Users/dtamayo/rebound.) +endif +PROBLEMDIR=$(shell basename `dirname \`pwd\``)"/"$(shell basename `pwd`) + +include $(REB_DIR)/src/Makefile.defs + +ASSIST_DIR=../../ + +all: librebound.so libassist.so + @echo "" + @echo "Compiling problem file ..." + $(CC) -I$(ASSIST_DIR)/src/ -I$(REB_DIR)/src/ -Wl,-rpath,./ $(OPT) $(PREDEF) problem.c -L. -lassist -lrebound $(LIB) -o rebound + @echo "" + @echo "Problem file compiled successfully." + +librebound.so: + @echo "Compiling shared library librebound.so ..." + $(MAKE) -C $(REB_DIR)/src/ + @echo "Creating link for shared library librebound.so ..." + @-rm -f librebound.so + @ln -s $(REB_DIR)/src/librebound.so . + +libassist.so: librebound.so $(ASSIST_DIR)/src/*.h $(ASSIST_DIR)/src/*.c + @echo "Compiling shared library libassist.so ..." + $(MAKE) -C $(ASSIST_DIR)/src/ + @-rm -f libassist.so + @ln -s $(ASSIST_DIR)/src/libassist.so . + +clean: + @echo "Cleaning up shared library librebound.so ..." + @-rm -f librebound.so + $(MAKE) -C $(REB_DIR)/src/ clean + @echo "Cleaning up shared library libassist.so ..." + @-rm -f libassist.so + $(MAKE) -C $(ASSIST_DIR)/src/ clean + @echo "Cleaning up local directory ..." + @-rm -vf rebound + @-rm -vf test_*.txt test_*.bin test_*.bsp + +test: all + @echo "Running SPK detection tests..." + ./rebound + @echo "All tests completed!" + +.PHONY: all clean test \ No newline at end of file diff --git a/unit_tests/spk_detection/problem.c b/unit_tests/spk_detection/problem.c new file mode 100644 index 0000000..f208516 --- /dev/null +++ b/unit_tests/spk_detection/problem.c @@ -0,0 +1,142 @@ +#include +#include +#include +#include +#include +#include +#include +#include "../../src/spk.h" + +// Embedded file samples from real ephemeris files + +// Valid BSP file header (first 64 bytes from de440.bsp) +static const unsigned char valid_bsp_sample[] = { + 0x44, 0x41, 0x46, 0x2f, 0x53, 0x50, 0x4b, 0x20, 0x02, 0x00, 0x00, 0x00, 0x06, 0x00, 0x00, 0x00, + 0x4e, 0x49, 0x4f, 0x32, 0x53, 0x50, 0x4b, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, + 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, + 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20 +}; + +// Legacy binary file sample (301 bytes from linux_p1550p2650.440 including signature at 0xFC) +static const unsigned char legacy_binary_sample[] = { + 0x4a, 0x50, 0x4c, 0x20, 0x50, 0x6c, 0x61, 0x6e, 0x65, 0x74, 0x61, 0x72, 0x79, 0x20, 0x45, 0x70, + 0x68, 0x65, 0x6d, 0x65, 0x72, 0x69, 0x73, 0x20, 0x44, 0x45, 0x34, 0x34, 0x30, 0x2f, 0x4c, 0x45, + 0x34, 0x34, 0x30, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, + 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, + 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, + 0x20, 0x20, 0x20, 0x20, 0x53, 0x74, 0x61, 0x72, 0x74, 0x20, 0x45, 0x70, 0x6f, 0x63, 0x68, 0x3a, + 0x20, 0x4a, 0x45, 0x44, 0x3d, 0x20, 0x20, 0x32, 0x32, 0x38, 0x37, 0x31, 0x38, 0x34, 0x2e, 0x35, + 0x20, 0x20, 0x31, 0x35, 0x34, 0x39, 0x2d, 0x44, 0x45, 0x43, 0x2d, 0x32, 0x31, 0x20, 0x30, 0x30, + 0x3a, 0x30, 0x30, 0x3a, 0x30, 0x30, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, + 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, + 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x46, 0x69, 0x6e, 0x61, 0x6c, 0x20, 0x45, 0x70, + 0x6f, 0x63, 0x68, 0x3a, 0x20, 0x4a, 0x45, 0x44, 0x3d, 0x20, 0x20, 0x32, 0x36, 0x38, 0x38, 0x39, + 0x37, 0x36, 0x2e, 0x35, 0x20, 0x20, 0x32, 0x36, 0x35, 0x30, 0x2d, 0x4a, 0x41, 0x4e, 0x2d, 0x32, + 0x35, 0x20, 0x30, 0x30, 0x3a, 0x30, 0x30, 0x3a, 0x30, 0x30, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, + 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, + 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x44, 0x45, 0x4e, 0x55, + 0x4d, 0x20, 0x4c, 0x45, 0x4e, 0x55, 0x4d, 0x20, 0x54, 0x44, 0x41, 0x54, 0x45, 0x46, 0x54, 0x44, + 0x41, 0x54, 0x45, 0x42, 0x4a, 0x44, 0x45, 0x50, 0x4f, 0x43, 0x43, 0x45, 0x4e, 0x54, 0x45, 0x52, + 0x43, 0x4c, 0x49, 0x47, 0x48, 0x54, 0x42, 0x45, 0x54, 0x41, 0x20, 0x20, 0x47, 0x4d, 0x53, 0x20, + 0x20, 0x20, 0x47, 0x4d, 0x31, 0x20, 0x20, 0x20, 0x47, 0x4d, 0x32, 0x20, 0x20, 0x20 +}; + +// ASCII source file sample +static const char ascii_source_sample[] = "JPL Planetary EPHEMERIS DE440 ASCII source file\n"; + +static int create_temp_file(const char* filename, const void* data, size_t size) { + int fd = open(filename, O_CREAT | O_WRONLY | O_TRUNC, 0644); + if (fd < 0) return -1; + + if (write(fd, data, size) != (ssize_t)size) { + close(fd); + return -1; + } + + close(fd); + return 0; +} + +int main() { + int fd; + ephemeris_file_format_t format; + struct spk_s* spk; + + // Test valid BSP file detection + assert(create_temp_file("test_valid.bsp", valid_bsp_sample, sizeof(valid_bsp_sample)) == 0); + fd = open("test_valid.bsp", O_RDONLY); + assert(fd >= 0); + format = assist_detect_ephemeris_file_format(fd); + close(fd); + assert(format == FILE_FORMAT_VALID_BSP); + unlink("test_valid.bsp"); + + // Test legacy binary file detection + assert(create_temp_file("test_legacy.440", legacy_binary_sample, sizeof(legacy_binary_sample)) == 0); + fd = open("test_legacy.440", O_RDONLY); + assert(fd >= 0); + format = assist_detect_ephemeris_file_format(fd); + close(fd); + assert(format == FILE_FORMAT_BINARY_LEGACY); + + // Test legacy binary file rejection through ASSIST API + spk = assist_spk_init("test_legacy.440"); + assert(spk == NULL); + unlink("test_legacy.440"); + + // Test ASCII source file detection + assert(create_temp_file("test_ascii.440", ascii_source_sample, strlen(ascii_source_sample)) == 0); + fd = open("test_ascii.440", O_RDONLY); + assert(fd >= 0); + format = assist_detect_ephemeris_file_format(fd); + close(fd); + assert(format == FILE_FORMAT_ASCII_SOURCE); + + // Test ASCII source file rejection through ASSIST API + spk = assist_spk_init("test_ascii.440"); + assert(spk == NULL); + unlink("test_ascii.440"); + + // Test unknown file format + const char unknown_sample[] = "Not an ephemeris file"; + assert(create_temp_file("test_unknown.dat", unknown_sample, strlen(unknown_sample)) == 0); + fd = open("test_unknown.dat", O_RDONLY); + assert(fd >= 0); + format = assist_detect_ephemeris_file_format(fd); + close(fd); + assert(format == FILE_FORMAT_UNKNOWN); + unlink("test_unknown.dat"); + + // Test signature detection function directly + assert(create_temp_file("test_sig_valid.bsp", valid_bsp_sample, sizeof(valid_bsp_sample)) == 0); + fd = open("test_sig_valid.bsp", O_RDONLY); + assert(assist_detect_legacy_binary_signature(fd) == 0); + close(fd); + unlink("test_sig_valid.bsp"); + + assert(create_temp_file("test_sig_legacy.440", legacy_binary_sample, sizeof(legacy_binary_sample)) == 0); + fd = open("test_sig_legacy.440", O_RDONLY); + assert(assist_detect_legacy_binary_signature(fd) == 1); + close(fd); + unlink("test_sig_legacy.440"); + + // Test with real BSP file if available + const char* test_files[] = {"../../data/de440.bsp", "/tmp/de440.bsp", "./de440.bsp"}; + for (int i = 0; i < 3; i++) { + fd = open(test_files[i], O_RDONLY); + if (fd >= 0) { + format = assist_detect_ephemeris_file_format(fd); + close(fd); + assert(format == FILE_FORMAT_VALID_BSP); + + spk = assist_spk_init(test_files[i]); + if (spk) { + assert(spk->num == 14); // Expected target count + assist_spk_free(spk); + } + break; + } + } + + return 0; +} \ No newline at end of file