From ce03e944918d0818e854cfed4dc55650a293aa86 Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Thu, 12 Sep 2024 18:05:54 -0400 Subject: [PATCH] Implement SPICE (bsp) for planet ephemeris --- .github/workflows/c.yml | 29 +- .github/workflows/python.yml | 4 +- .python-version | 1 + MANIFEST.in | 2 - Makefile | 22 + README.md | 31 +- assist/ephem.py | 20 +- 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 | 121 ++--- src/assist.h | 10 +- src/forces.c | 103 ++-- src/get_mass | Bin 53464 -> 0 bytes src/get_mass.c | 33 -- src/planets.c | 363 ------------- src/planets.h | 74 --- src/spk.c | 486 ++++++++++++++++-- src/spk.h | 72 +++ unit_tests/5303_Ceres/problem.c | 2 +- unit_tests/apophis/problem.c | 6 +- unit_tests/benchmark/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 | 2 +- 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 | 58 +++ unit_tests/spk_load_constants/Makefile | 46 ++ unit_tests/spk_load_constants/problem.c | 39 ++ unit_tests/variational/problem.c | 2 +- 51 files changed, 1074 insertions(+), 741 deletions(-) create mode 100644 .python-version 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 f85d353..c4286d5 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/ @@ -91,3 +83,18 @@ jobs: run: | make LD_LIBRARY_PATH=../../../rebound/src/ ./rebound + - name: Test Loading Spice Kernels + working-directory: ./unit_tests/spk_init/ + run: | + make + LD_LIBRARY_PATH=../../../rebound/src/ ./rebound + - name: Test Loading Constants from Planet SPK + working-directory: ./unit_tests/spk_load_constants/ + run: | + make + LD_LIBRARY_PATH=../../../rebound/src/ ./rebound + - name: Test Joining Mass Data from SPK + working-directory: ./unit_tests/spk_join_masses/ + run: | + make + LD_LIBRARY_PATH=../../../rebound/src/ ./rebound \ 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/.python-version b/.python-version new file mode 100644 index 0000000..686c879 --- /dev/null +++ b/.python-version @@ -0,0 +1 @@ +assist diff --git a/MANIFEST.in b/MANIFEST.in index 8bc4723..b5fb7ec 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,10 +1,8 @@ include src/assist.c include src/forces.c -include src/planets.c include src/spk.c include src/assist.h include src/forces.h -include src/planets.h include src/spk.h include README.md include LICENSE diff --git a/Makefile b/Makefile index de22331..ea1a995 100644 --- a/Makefile +++ b/Makefile @@ -18,3 +18,25 @@ doc: cd doc/doxygen && doxygen $(MAKE) -C doc html + +# Iterate through each directory in the unit_tests directory and compile the test files +# then run them +test: + @echo "Running tests ..." + @set -e; \ + for dir in $(wildcard unit_tests/*); do \ + if [ -d $$dir ]; then \ + echo "Entering directory $$dir"; \ + if $(MAKE) -C $$dir; then \ + if [ -f "$$dir/rebound" ]; then \ + (cd $$dir && LD_LIBRARY_PATH=$(REB_DIR):. ./rebound); \ + else \ + echo "No rebound executable found in $$dir"; \ + exit 1; \ + fi \ + else \ + echo "Make failed in $$dir"; \ + exit 1; \ + fi \ + fi \ + done \ 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..4229373 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,10 @@ def get_particle(self, body: Union[int, str], t: float) -> rebound.Particle: def __del__(self) -> None: clibassist.assist_ephem_free_pointers(byref(self)) - _fields_ = [("jd_ref", c_double), - ("_pl", c_void_p), - ("_spl", c_void_p), - ] + _fields_ = [ + ("jd_ref", c_double), + ("spk_global", c_void_p), + ("spk_planets", c_void_p), + ("spk_asteroids", c_void_p), + ] diff --git a/assist/extras.py b/assist/extras.py index ea80243..747483b 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 3ae4524..c86e7c5 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'], + sources = [ 'src/assist.c', 'src/spk.c', 'src/forces.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'], + sources = [ 'src/assist.c','src/spk.c', 'src/forces.c'], include_dirs = ['src'], library_dirs = [], runtime_library_dirs = ["."], diff --git a/src/Makefile b/src/Makefile index 0551f13..d3079bb 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 +SOURCES=assist.c spk.c forces.c OBJECTS=$(SOURCES:.c=.o) HEADERS=$(SOURCES:.c=.h) diff --git a/src/assist.c b/src/assist.c index 1ba197b..23237a4 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,9 +113,9 @@ 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; - } + ephem->spk_planets = assist_spk_init(planets_path); + ephem->spk_global = assist_load_spk_constants(planets_path); + assist_spk_join_masses(ephem->spk_planets, ephem->spk_global); int asteroids_path_not_found = 0; @@ -132,87 +131,14 @@ int assist_ephem_init(struct assist_ephem* ephem, char *user_planets_path, char } if (asteroids_path_not_found != 1){ - if ((ephem->spl = assist_spk_init(asteroids_path)) == NULL) { + if ((ephem->spk_asteroids = assist_spk_init(asteroids_path)) == NULL) { asteroids_path_not_found = 1; } } + // Join the mass data from whichever planets file we are using. if (asteroids_path_not_found != 1){ - // Try to find masses of bodies in spk file in ephemeris constants - for(int n=0; 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 ); - } - - } + assist_spk_join_masses(ephem->spk_asteroids, ephem->spk_global); }else{ fprintf(stderr, "(ASSIST) %s\n", assist_error_messages[ASSIST_ERROR_AST_FILE]); } @@ -233,13 +159,17 @@ 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->spk_asteroids != NULL){ + assist_spk_free(ephem->spk_asteroids); } - if (ephem->spl){ - assist_spk_free(ephem->spl); + if (ephem->spk_global != NULL){ + assist_free_spk_constants(ephem->spk_global); } } + void assist_ephem_free(struct assist_ephem* ephem){ assist_ephem_free_pointers(ephem); free(ephem); @@ -278,8 +208,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)); @@ -485,8 +415,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); @@ -578,3 +511,17 @@ static void assist_pre_timestep_modifications(struct reb_simulation* sim){ memcpy(assist->last_state, sim->particles, sizeof(struct reb_particle)*sim->N); } +double assist_get_constant(struct assist_ephem* ephem, const char* constant_name) { + struct spk_constants* con = &ephem->spk_global->con; + if (strcmp(constant_name, "AU") == 0) return con->AU; + if (strcmp(constant_name, "EMRAT") == 0) return con->EMRAT; + if (strcmp(constant_name, "J2E") == 0) return con->J2E; + if (strcmp(constant_name, "J3E") == 0) return con->J3E; + if (strcmp(constant_name, "J4E") == 0) return con->J4E; + if (strcmp(constant_name, "J2SUN") == 0) return con->J2SUN; + if (strcmp(constant_name, "RE") == 0) return con->RE; + if (strcmp(constant_name, "CLIGHT") == 0) return con->CLIGHT; + if (strcmp(constant_name, "ASUN") == 0) return con->ASUN; + // If the constant is not found, return some error value (NaN) + return NAN; +} diff --git a/src/assist.h b/src/assist.h index 00470fc..3bcb17f 100644 --- a/src/assist.h +++ b/src/assist.h @@ -88,8 +88,9 @@ enum ASSIST_BODY { struct assist_ephem { double jd_ref; - struct jpl_s* jpl; - struct spk_s* spl; + struct spk_global* spk_global; + struct spk_s* spk_planets; + struct spk_s* spk_asteroids; }; struct assist_cache_item { @@ -181,6 +182,7 @@ struct reb_particle assist_get_particle_with_error(struct assist_ephem* ephem, c // 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); @@ -194,12 +196,12 @@ 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. */ /// int assist_ephem_init(struct assist_ephem* ephem, char *user_planets_path, char *user_asteroids_path); - +double assist_get_constant(struct assist_ephem* ephem, const char* name); #endif diff --git a/src/forces.c b/src/forces.c index 1d621d5..305ddbd 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){ @@ -129,7 +127,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, @@ -202,19 +199,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 @@ -290,16 +324,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 } @@ -317,6 +351,7 @@ static void assist_additional_force_direct(struct reb_simulation* sim, double xo reb_simulation_error(sim, assist_error_messages[flag]); } + // 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 } @@ -456,13 +491,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 = assist_get_constant(assist->ephem, "J2E"); + double J3e = assist_get_constant(assist->ephem, "J3E"); + double J4e = assist_get_constant(assist->ephem, "J4E"); + double Re_eq = assist_get_constant(assist->ephem, "RE")/assist_get_constant(assist->ephem, "AU"); + // Unit vector to equatorial pole at the epoch // Note also that the pole orientation is not changing during @@ -591,7 +628,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; @@ -631,6 +667,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){ @@ -655,9 +692,9 @@ 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 au = assist_get_constant(assist->ephem, "AU"); + double Rs_eq = assist_get_constant(assist->ephem, "ASUN")/au; + double J2s = assist_get_constant(assist->ephem, "J2SUN"); // Hard-coded constants. BEWARE! double RAs = 286.13*M_PI/180.; @@ -731,7 +768,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; @@ -761,11 +797,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, @@ -1068,9 +1101,8 @@ 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; + double au = assist_get_constant(assist->ephem, "AU"); + double c = (assist_get_constant(assist->ephem, "CLIGHT") / au) * 86400; const double C2 = c*c; const unsigned int N = sim->N; // N includes real+variational particles @@ -1175,9 +1207,8 @@ 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; + double au = assist_get_constant(assist->ephem, "AU"); + double c = (assist_get_constant(assist->ephem, "CLIGHT") / au) * 86400; const double C2 = c*c; @@ -1310,8 +1341,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; + double au = assist_get_constant(assist->ephem, "AU"); + double c = (assist_get_constant(assist->ephem, "CLIGHT") / au) * 86400; const double C2 = c*c; const double over_C2 = 1./(c*c); 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..e6e8330 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 * @@ -46,6 +44,190 @@ int assist_spk_free(struct spk_s *pl) return 0; } +/* + * assist_free_spk_constants + * + */ +int assist_free_spk_constants(struct spk_global *sg) { + + if (sg->masses.names) { + for (int i = 0; i < sg->masses.count; i++) { + free(sg->masses.names[i]); + } + free(sg->masses.names); + } + + free(sg->masses.values); + free(sg); + + return ASSIST_SUCCESS; +} + +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) { + char *d_char = strchr(str, 'D'); + if (d_char != NULL) { + // Ensure it's a scientific notation before replacing + if ((d_char > str && isdigit(*(d_char - 1))) && (isdigit(*(d_char + 1)) || (*(d_char + 1) == '+' || *(d_char + 1) == '-'))) { + *d_char = 'E'; + } + } +} + + +// Reads in the constants and masses from the comments of the DE440.bsp file +struct spk_global * assist_load_spk_constants(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); + + char *comments = NULL; + parse_comments(fd, record->file.fward, &comments); + + if (comments == NULL) { + printf("No comments found.\n"); + return NULL; + } + + struct spk_global* sg = malloc(sizeof(struct spk_global)); + if (!sg) { + // Handle memory allocation failure + close(fd); + return NULL; + } + + // Initialize the spk_global struct using designated initializers + *sg = (struct spk_global){ + .con = { + .AU = 0, + .EMRAT = 0, + .J2E = 0, + .J3E = 0, + .J4E = 0, + .J2SUN = 0, + .RE = 0, + .CLIGHT = 0, + .ASUN = 0 + }, + .masses = { + .names = NULL, + .values = NULL, + .count = 0 + } + }; + + 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); + if (strcmp(key, "cau") == 0) sg->con.AU = value; + else if (strcmp(key, "EMRAT") == 0) sg->con.EMRAT = value; + else if (strcmp(key, "J2E") == 0) sg->con.J2E = value; + else if (strcmp(key, "J3E") == 0) sg->con.J3E = value; + else if (strcmp(key, "J4E") == 0) sg->con.J4E = value; + else if (strcmp(key, "J2SUN") == 0) sg->con.J2SUN = value; + else if (strcmp(key, "AU") == 0) sg->con.AU = value; + else if (strcmp(key, "RE") == 0) sg->con.RE = value; + else if (strcmp(key, "CLIGHT") == 0) sg->con.CLIGHT = value; + else if (strcmp(key, "ASUN") == 0) sg->con.ASUN = value; + else if (strncmp(key, "GM", 2) == 0 || strncmp(key, "MA", 2) == 0) { + sg->masses.names = realloc(sg->masses.names, (sg->masses.count + 1) * sizeof(char *)); + sg->masses.values = realloc(sg->masses.values, (sg->masses.count + 1) * sizeof(double)); + sg->masses.names[sg->masses.count] = strdup(key); + sg->masses.values[sg->masses.count] = value; + sg->masses.count++; + } + } + } + + line = strtok(NULL, "\n"); + } + + free(comments); + close(fd); + return sg; +} /* * assist_spk_init @@ -56,74 +238,148 @@ 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 the global masses array +void assist_spk_join_masses(struct spk_s *sp, struct spk_global *sg) { + if (sp == NULL || sg == 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++) { + // 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); + } - // 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 - }; + // Find the mass in the masses array + for (int i = 0; i < sg->masses.count; i++) { + if (strcmp(sg->masses.names[i], mass_label) == 0) { + // Earth and moon mass is stored as one value + // Use the ratio to split it using pl->con.EMRAT + if (sp->targets[m].code == 399) { + sp->targets[m].mass = sg->masses.values[i] * (sg->con.EMRAT / (1. + sg->con.EMRAT)); + } else if (sp->targets[m].code == 301) { + sp->targets[m].mass = sg->masses.values[i] * (1./(1.+sg->con.EMRAT)); + } else { + sp->targets[m].mass = sg->masses.values[i]; + } + break; + } + } - // 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; + 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); + } + } +} - // Try opening file. - int fd = open(path, O_RDONLY); - if (fd < 0){ +// 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,8 +388,8 @@ 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; @@ -152,6 +408,8 @@ struct spk_s * assist_spk_init(const char *path) { pl->targets[m].one = calloc(32768, sizeof(int)); pl->targets[m].two = calloc(32768, sizeof(int)); pl->targets[m].ind = 0; + // Set default of mass to 0 + pl->targets[m].mass = 0; pl->num++; } @@ -163,14 +421,14 @@ struct spk_s * assist_spk_init(const char *path) { } // 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); } } @@ -279,3 +537,127 @@ 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(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(struct spk_s *pl, 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(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 = assist_get_constant(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; +} diff --git a/src/spk.h b/src/spk.h index 4de1a5f..51bf31e 100644 --- a/src/spk.h +++ b/src/spk.h @@ -10,6 +10,9 @@ #include "assist.h" +#define RECORD_LENGTH 1024 +#define MAX_COMMENT_LENGTH 1000000 + struct mpos_s { double u[3]; double v[3]; @@ -17,7 +20,33 @@ struct mpos_s { double jde; }; +struct spk_constants { + 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; +}; + +struct mass_data { + char **names; + double *values; + int count; +}; + +// Represents global values of constants and masses +// fetched from the comments of an SPK file (used with DE440) +struct spk_global { + struct spk_constants con; + 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 +60,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 +69,52 @@ 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); +int assist_free_spk_constants(struct spk_global *sg); struct spk_s * assist_spk_init(const char *path); +union record_t * assist_load_spk_file_record(int fd); +struct spk_global * assist_load_spk_constants(const char *path); +void parse_comments(int fd, int first_summary_record, char **comments); +void assist_spk_join_masses(struct spk_s *sp, struct spk_global *sg); enum ASSIST_STATUS assist_spk_calc(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(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); #endif // _SPK_H 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/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..d02e4a2 --- /dev/null +++ b/unit_tests/spk_join_masses/problem.c @@ -0,0 +1,58 @@ +#include +#include +#include +#include "spk.h" + + +void assert_populated_masses(struct spk_s *spk_data) { + if (!spk_data) { + printf("SPK data is NULL.\n"); + return; + } + + 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 spk_global *sg = assist_load_spk_constants(planets_path); + if (!sg) { + fprintf(stderr, "Failed to initialize SPK data.\n"); + return EXIT_FAILURE; + } + struct spk_s *planet_data = assist_spk_init(planets_path); + if (!planet_data) { + fprintf(stderr, "Failed to initialize SPK data.\n"); + return EXIT_FAILURE; + } + + assist_spk_join_masses(planet_data, sg); + assert_populated_masses(planet_data); + + struct spk_s *asteroid_data = assist_spk_init(asteroids_path); + if (!asteroid_data) { + fprintf(stderr, "Failed to initialize SPK data.\n"); + return EXIT_FAILURE; + } + + assist_spk_join_masses(asteroid_data, sg); + assert_populated_masses(asteroid_data); + + assist_spk_free(planet_data); + assist_spk_free(asteroid_data); + + assist_free_spk_constants(sg); + + 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..026b0f3 --- /dev/null +++ b/unit_tests/spk_load_constants/problem.c @@ -0,0 +1,39 @@ +#include +#include +#include +#include "spk.h" + + +void assert_constants_exist(struct spk_global *sg) { + struct spk_constants *sc = &sg->con; + assert(sc->AU != 0.0); + assert(sc->EMRAT != 0.0); + assert(sc->J2E != 0.0); + assert(sc->J3E != 0.0); + assert(sc->J4E != 0.0); + assert(sc->J2SUN != 0.0); + assert(sc->RE != 0.0); + assert(sc->CLIGHT != 0.0); + assert(sc->ASUN != 0.0); + assert(sg->masses.count = 420); + for (int i = 0; i < sg->masses.count; i++) { + assert(sg->masses.names[i] != NULL); + assert(sg->masses.values[i] != 0.0); + } + printf("All constants exist.\n"); +} + +int main(int argc, char *argv[]) { + + const char *planets_path = "../../data/de440.bsp"; + + struct spk_global *sg = assist_load_spk_constants(planets_path); + if (!sg) { + fprintf(stderr, "Failed to initialize SPK data.\n"); + return EXIT_FAILURE; + } + + assert_constants_exist(sg); + + 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");