From 9fabfcd9dd432b097ce83a897d27a9dcc7c13fea Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Sun, 10 Nov 2024 22:19:31 +0100 Subject: [PATCH] CSPICE integration (#88) * Adjust solsys-calcepth tests * Initial CSPICE support * Tweaks to CSPICE support and test * Documentation for CSPICE support --- CHANGELOG.md | 8 + Makefile | 12 ++ README.md | 86 +++++++-- config.mk | 8 + include/solarsystem.h | 10 + src/solsys-calceph.c | 2 +- src/solsys-cspice.c | 395 ++++++++++++++++++++++++++++++++++++++++ test/Makefile | 13 ++ test/src/test-calceph.c | 35 ++-- test/src/test-cspice.c | 208 +++++++++++++++++++++ 10 files changed, 745 insertions(+), 32 deletions(-) create mode 100644 src/solsys-cspice.c create mode 100644 test/src/test-cspice.c diff --git a/CHANGELOG.md b/CHANGELOG.md index 7fb8d0df..dff89881 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,14 @@ Changes expected for the next feature release, expected around 1 February 2025. the CALCEPH development files (C headers and unversioned static or shared libraries depending on the needs of the build). + - #86: NAIF CSPICE integration: `novas_use_cspice()`, `novas_use_cspice_planets()`, `novas_use_cspice_ephem()` + to use the NAIF CSPICE library for all Solar-system sources, major planets only, or for other bodies only. + `NOVAS_EPHEM_OBJECTS` should use NAIF IDs with CSPICE (or else -1 for name-based lookup). These functions are + provided by the `libsolsys-cspice.so[.1]` and/or `.a` plugin libraries, which are built contingent on the + `CSPICE_SUPPORT` variable being set to 1 prior to the build. The build of the plugin module requires an accessible + installation of the CSPICE development files (C headers and unversioned static or shared libraries depending on the + needs of the build). + - #87: Added `novas_planet_for_name()` function to return the NOVAS planet ID for a given (case insensitive) name. - NOVAS-NAIF conversions for major planets (and Sun, Moon, SSB): `novas_to_naif_planet()` (planet centers), diff --git a/Makefile b/Makefile index a5b07eda..75281682 100644 --- a/Makefile +++ b/Makefile @@ -61,6 +61,12 @@ ifeq ($(CALCEPH_SUPPORT),1) SHARED_TARGETS += $(LIB)/libsolsys-calceph.so endif +ifeq ($(CSPICE_SUPPORT),1) + CPPFLAGS += -DUSE_CSPICE=1 + SOLSYS_TARGETS += $(OBJ)/solsys-cspice.o + SHARED_TARGETS += $(LIB)/libsolsys-cspice.so +endif + # Default target for packaging with Linux distributions .PHONY: distro distro: $(SHARED_TARGETS) $(DOC_TARGETS) @@ -118,6 +124,8 @@ $(LIB)/libsolsys-ephem.so: $(LIB)/libsolsys-ephem.so.$(SO_VERSION) $(LIB)/libsolsys-calceph.so: $(LIB)/libsolsys-calceph.so.$(SO_VERSION) +$(LIB)/libsolsys-cspice.so: $(LIB)/libsolsys-cspice.so.$(SO_VERSION) + $(LIB)/libnovas.so: $(LIB)/libsupernovas.so $(LIB)/libsolsys%.so.$(SO_VERSION): LDFLAGS += -L$(LIB) -lsupernovas @@ -145,6 +153,10 @@ $(LIB)/libsolsys-ephem.so.$(SO_VERSION): $(SRC)/solsys-ephem.c | $(LIB)/libsuper $(LIB)/libsolsys-calceph.so.$(SO_VERSION): LDFLAGS += -lcalceph $(LIB)/libsolsys-calceph.so.$(SO_VERSION): $(SRC)/solsys-calceph.c | $(LIB)/libsupernovas.so +# Shared library: libsolsys-cspice.so.1 (standalone solsys2.c functionality) +$(LIB)/libsolsys-cspice.so.$(SO_VERSION): LDFLAGS += -lcspice +$(LIB)/libsolsys-cspice.so.$(SO_VERSION): $(SRC)/solsys-cspice.c | $(LIB)/libsupernovas.so + # Static library: libnovas.a $(LIB)/libnovas.a: $(OBJECTS) | $(LIB) Makefile diff --git a/README.md b/README.md index 2a00e449..e7654b8a 100644 --- a/README.md +++ b/README.md @@ -87,6 +87,8 @@ Outside contributions are very welcome. See from JPL and/or in INPOP 2.0/3.0 format. - [NAIF SPICE toolkit](https://naif.jpl.nasa.gov/naif/toolkit.html) for integrating Solar-system ephemeris via from JPL. + - [Smithsonian/cspice-sharedlib](https://github.com/Smithsonian/cspice-sharedlib) for building CSPICE as a shared + library for dynamic linking. - [IAU Minor Planet Center](https://www.minorplanetcenter.net/iau/mpc.html) provides another source of ephemeris data. @@ -201,8 +203,15 @@ the necessary variables in the shell prior to invoking `make`. For example: setting `CALCEPH_SUPPORT = 1` in `config.mk` or in the shell prior to the build. When enabled it will build `libsolsys-calceph.so[.1]` and/or `.a` supplemental libraries, depending on the build target. The build of the modules requires an accessible installation of the CALCEPH development files (C headers and unversioned static or - shared libraries depending on the needs of the build). You might want to set `LD_LIBRARY_PATH`, and/or `CPPFLAGS` - to include an appropriate `-I` option, to help locate these prior to calling `make`. + shared libraries depending on the needs of the build). + + - You can enable integration with the [NAIF CSPICE Toolkit](https://naif.jpl.nasa.gov/naif/toolkit.html), by setting + `CSPICE_SUPPORT = 1` in `config.mk` or in the shell prior to the build. When enabled it will build + `libsolsys-cspice.so[.1]` and/or `.a` supplemental libraries, depending on the build target. The build of the + modules requires an accessible installation of the CSPICE development files (C headers, under a `cspice/` + sub-folder in the header search path, and unversioned static or shared libraries depending on the needs of the + build). You might want to check out the [Smithsonian/cspice-sharedlib](https://github.com/Smithsonian/cspice-sharedlib) + repository for building CSPICE as a shared library. - If you want to use the CIO locator binary file for `cio_location()`, you can specify the path to the CIO locator file (e.g. `/usr/local/share/supernovas/CIO_RA.TXT`) on your system e.g. by setting the `CIO_LOCATOR_FILE` shell @@ -814,6 +823,10 @@ before that level of accuracy is reached. to specify whether `object.number` in `NOVAS_EPHEM_OBJECT` type objects is a NAIF ID (default) or else a CALCEPH ID number of the Solar-system body. + - NAIF CSPICE integration: `novas_use_cspice()`, `novas_use_cspice_planets()`, `novas_use_cspice_ephem()` to use the + NAIF CSPICE library for all Solar-system sources, major planets only, or for other bodies only. + `NOVAS_EPHEM_OBJECTS` should use NAIF IDs with CSPICE (or else -1 for name-based lookup). + - NAIF/NOVAS ID conversions for major planets (and Sun, Moon, SSB): `novas_to_naif_planet()`, `novas_to_dexxx_planet()`, and `naif_to_novas_planet()`. @@ -821,6 +834,8 @@ before that level of accuracy is reached. - Added `novas_planet_for_name()` function to return the NOVAS planet ID for a given (case insensitive) name. + + ### Refinements to the NOVAS C API @@ -913,8 +928,9 @@ NASA/JPL also provides [generic ephemerides](https://naif.jpl.nasa.gov/pub/naif/ planets, satellites thereof, the 300 largest asteroids, the Lagrange points, and some Earth orbiting stations. - [Universal ephemeris data / service integration](#universal-ephemerides) - - [Built-in support for CALCEPH integration](#calceph-integration) - - [Built-in support for (old) JPL major planet ephemerides](#builtin-ephem-readers) + - [Optional support for CALCEPH integration](#calceph-integration) + - [Optional support for NAIF CSPICE toolkit integration](#cspice-integration) + - [Alternative support for (older) JPL major planet ephemerides](#builtin-ephem-readers) - [Explicit linking of custom ephemeris functions](#explicit-ephem-linking) @@ -964,18 +980,19 @@ provided you compiled SuperNOVAS with `BUILTIN_SOLSYS_EPHEM = 1` (in `config.mk` -### Built-in support for CALCEPH integration +### Optional support for CALCEPH integration The [CALCEPH](https://www.imcce.fr/recherche/equipes/asd/calceph/) library provides an easy-to-use access to JPL and -INPOP ephemeris files from C/C++. As of version 1.2, we provide built-in support for integrating the CALCEPH C library +INPOP ephemeris files from C/C++. As of version 1.2, we provide optional support for interfacing the CALCEPH C library with SuperNOVAS for handling Solar-system objects. -Prior to building SuperNOVAS simply set `CALCEPH_SUPPORT` to 1 in `config.mk` or in your shell. Depending on the +Prior to building SuperNOVAS simply set `CALCEPH_SUPPORT` to 1 in `config.mk` or in your environment. Depending on the build target, it will build `libsolsys-calceph.so[.1]` (target `shared`) or `libsolsys-calceph.a` (target `static`) libraries, which provide the `novas_use_calceph()` and `novas_use_calceph_planets()` functions. -Of course, you will need access to the CALCEPH C development files (C headers and unversioned libraries) for the build -to succeed. Here is an example on how you'd use CALCEPH with SuperNOVAS in your application code: +Of course, you will need access to the CALCEPH C development files (C headers and unversioned `libcalceph.so` or `.a` +library) for the build to succeed. Here is an example on how you'd use CALCEPH with SuperNOVAS in your application +code: ```c #include @@ -1003,16 +1020,63 @@ You can obtain readily available standard ephemeris files from [NASA/JPL](https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/), such as [DE440](https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440.bsp), which covers the major planets, plus Sun, Moon, and Solar-System Barycenter (SSB) for times between 1550 AD and 2650 AD. Or, you can use the -[JPL HORIZONS](https://ssd.jpl.nasa.gov/horizons/app.html#/) system to generate cutom ephemeris data for pretty much +[JPL HORIZONS](https://ssd.jpl.nasa.gov/horizons/app.html#/) system to generate custom ephemeris data for pretty much all known solar systems bodies, down to the tiniest rocks. All of these should work with the `solsys-calceph` plugin. And, when linking your application, don't forget to add `-lsolsys-calceph` to your link flags. That's all there is to it. + +### Optional support for NAIF CSPICE toolkit integration + +The [NAIF CSPICE Toolkit](https://naif.jpl.nasa.gov/naif/toolkit.html) is the canonical standard library for JPL +ephemeris files from C/C++. As of version 1.2, we provide optional support for interfacing CSPICE with SuperNOVAS for +handling Solar-system objects. + +Prior to building SuperNOVAS simply set `CSPICE_SUPPORT` to 1 in `config.mk` or in your environment. Depending on the +build target, it will build `libsolsys-cspice.so[.1]` (target `shared`) or `libsolsys-cspice.a` (target `static`) +libraries, which provide the `novas_use_cspice()`, `novas_use_cspice_planets()`, and `novas_use_cspice_ephem()` +functions to enable CSPICE for providing data for all Solar-system sources, or for major planets only, or for other +bodies only, respectively. + +Of course, you will need access to the CSPICE development files (C headers, installed under a `cspice/` directory +of an header search location, and the unversioned `libcspice.so` or `.a` library) for the build to succeed. You may +want to check out the [Smithsonian/cspice-sharedlib](https://github.com/Smithsonian/cspice-sharedlib) GitHub +repository to help you build CSPICE with shared libraries and dynamically linked tools. + +Here is an example on how you might use CSPICE with SuperNOVAS in your application code: + +```c + // You can load the desired kernels for CSPICE, using the CSPICE API. + // E.g. load DE440s and the Mars satellites from /data/ephem: + furnsh_c("/data/ephem/de440s.bsp"); + furnsh_c("/data/ephem/mar097.bsp"); + ... + + // check for CSPICE errors + if(return_c()) { + // oops, one of the kernels must not have loaded... + ... + } + + // Then use CSPICE as your SuperNOVAS ephemeris provider + novas_use_cspice(); +``` + +You can obtain readily available standard ephemeris files from +[NASA/JPL](https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/), such as +[DE440](https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440.bsp), which covers the major planets, plus +Sun, Moon, and Solar-System Barycenter (SSB) for times between 1550 AD and 2650 AD. Or, you can use the +[JPL HORIZONS](https://ssd.jpl.nasa.gov/horizons/app.html#/) system to generate custom ephemeris data for pretty much +all known solar systems bodies, down to the tiniest rocks. All of these will work with the `solsys-cspice` plugin. + +And, when linking your application, don't forget to add `-lsolsys-cspice` to your link flags. That's all there is to +it. + -### Built-in support for (old) JPL major planet ephemerides +### Alternative support for (older) JPL major planet ephemerides If you only need support for major planets, you may be able to use one of the modules included in the SuperNOVAS distribution. The modules `solsys1.c` and `solsys2.c` provide built-in support to older JPL ephemerides (DE200 to DE421), diff --git a/config.mk b/config.mk index 8a69c285..8324f941 100644 --- a/config.mk +++ b/config.mk @@ -103,6 +103,14 @@ DEFAULT_SOLSYS ?= 3 CALCEPH_SUPPORT ?= 0 +# Whether or not to build solsys-cspice libraries. You need the CSPICE +# development libraries (libcspice.so and/or libcspice.a) installed in +# LD_LIBRARY_PATH, and CSPICE header files in /usr/include/cspice or some +# other accessible location (you may also set an appropriate -I +# option to CPPFLAGS prior to calling make). +CSPICE_SUPPORT ?= 0 + + # cppcheck options for 'check' target. You can add additional options by # setting the CHECKEXTRA variable (e.g. in shell) prior to invoking 'make'. CHECKOPTS ?= --enable=performance,warning,portability,style --language=c \ diff --git a/include/solarsystem.h b/include/solarsystem.h index 7acea5e7..f6635155 100644 --- a/include/solarsystem.h +++ b/include/solarsystem.h @@ -342,6 +342,16 @@ int novas_calceph_use_ids(enum novas_id_type idtype); #endif /* USE_CALCEPH */ +// in solsys-cspice.c +#if USE_CSPICE + +int novas_use_cspice(); + +int novas_use_cspice_ephem(); + +int novas_use_cspice_planets(); + +#endif /* USE_CALCEPH */ /// \cond PRIVATE diff --git a/src/solsys-calceph.c b/src/solsys-calceph.c index 13f0b9a2..ae3444fa 100644 --- a/src/solsys-calceph.c +++ b/src/solsys-calceph.c @@ -316,7 +316,7 @@ static int novas_calceph(const char *name, long id, double jd_tdb_high, double j id = i; } - // Always return psoitions and velocities w.r.t. the SSB + // Always return positions and velocities w.r.t. the SSB if(origin) *origin = NOVAS_SSB; diff --git a/src/solsys-cspice.c b/src/solsys-cspice.c new file mode 100644 index 00000000..e1aa9ee3 --- /dev/null +++ b/src/solsys-cspice.c @@ -0,0 +1,395 @@ +/** + * @file + * + * @author A. Kovacs + * + * SuperNOVAS major planet ephemeris lookup implementation via the NAIF CSPICE library + * See https://naif.jpl.nasa.gov/naif/toolkit.html + * + */ + +#include +#include +#include + +/// \cond PRIVATE +#define USE_CSIPCE 1 ///< NOVAS cspice integration prototypes +#define __NOVAS_INTERNAL_API__ ///< Use definitions meant for internal use by SuperNOVAS only +#include "novas.h" +#include "naif.h" + +#include "cspice/SpiceUsr.h" +#include "cspice/SpiceZpr.h" // for reset_c + +#define SPICE_FILE_SIZE 128 ///< (bytes) Maximum length of SPICE file names. +#define SPICE_WORD_SIZE 80 ///< (bytes) Maximum length for SPICE strings + +/// Multiplicative normalization for the positions returned by km to AU +#define NORM_POS (1e3 / NOVAS_AU) + +/// Multiplicative normalization for the velocities returned by km/s to AU/day +#define NORM_VEL (NORM_POS * 86400.0) +/// \endcond + +/// Semaphore for thread-safe access of ephemerides +static sem_t *sem; + +static int mutex_lock() { + if(!sem) { + sem = (sem_t *) calloc(1, sizeof(sem_t)); + if(!sem) { + perror("ERROR! solsys-cspice: alloc sem_t"); + exit(errno); + } + sem_init(sem, 0, 1); + } + + if(sem_wait(sem) != 0) + return novas_error(-1, errno, "mutex_lock()", "sem_wait()"); + + return 0; +} + +static int mutex_unlock() { + sem_post(sem); + return 0; +} + +/** + * Provides an interface between the NAIF CSPICE C library and NOVAS-C for regular (reduced) + * precision applications. The user must set the cspice ephemeris binary data to use using the + * novas_use_cspice() or novas_use_cspice_planet() to activate CSPICE as the NOVAS ephemeris + * provider. + * + * This call is generally thread safe (notwithstanding outside access to the ephemeris files), + * even if CSPICE itself may not be. All ephemeris access will be mutexed to ensure sequential + * access under the hood. + * + * The call will use whatever ephemeris (SPK) files were loaded by the CSPICE library prior + * to the call (see furnsh_c() function) + * + * REFERENCES: + *
    + *
  1. NAIF CSPICE: https://naif.jpl.nasa.gov/naif/toolkit.html
    2. + *
  2. Kaplan, G. H. "NOVAS: Naval Observatory Vector Astrometry + * Subroutines"; USNO internal document dated 20 Oct 1988; + * revised 15 Mar 1990.
    3. + *
+ * + * @param jd_tdb [day] Two-element array containing the Julian date, which may be + * split any way (although the first element is usually the + * "integer" part, and the second element is the "fractional" + * part). Julian date is on the TDB or "T_eph" time scale. + * @param body Major planet number (or that for Sun, Moon, or Solar-system barycenter) + * @param origin NOVAS_BARYCENTER (0) or NOVAS_HELIOCENTER (1) + * -- relative to which to report positions and velocities. + * @param[out] position [AU] Position vector of 'body' at jd_tdb; equatorial rectangular + * coordinates in AU referred to the ICRS. + * @param[out] velocity [AU/day] Velocity vector of 'body' at jd_tdb; equatorial rectangular + * system referred to the ICRS, in AU/day. + * @return 0 if successful, or else 1 if the 'body' is invalid, or 2 if the + * 'origin' is invalid, or 3 if there was an error providing ephemeris + * data. + * + * @sa planet_cspice + * @sa novas_use_cspice() + * @sa novas_use_cspice_planet() + * + * @author Attila Kovacs + * @since 1.2 + */ +static short planet_cspice_hp(const double jd_tdb[2], enum novas_planet body, enum novas_origin origin, double *position, + double *velocity) { + static const char *fn = "planet_cspice_hp"; + + SpiceDouble pv[6]; + SpiceDouble lt; + SpiceInt target, center; + double tdb2000; // seconds past J2000 TDB + int i, err; + + if(!jd_tdb) + return novas_error(-1, EINVAL, fn, "jd_tdb input time array is NULL."); + + target = novas_to_naif_planet(body); + if(target < 0) return novas_trace(fn, 1, 0); + + switch(origin) { + case NOVAS_BARYCENTER: + center = NAIF_SSB; + break; + case NOVAS_HELIOCENTER: + center = NAIF_SUN; + break; + default: + return novas_error(2, EINVAL, fn, "Invalid origin type: %d", origin); + } + + tdb2000 = (jd_tdb[0] + jd_tdb[1] - NOVAS_JD_J2000) * 86400.0; + + prop_error(fn, mutex_lock(sem), 0); + + // Try with proper planet center NAIF ID first... + // + // See https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/frames.html#Reference%20Frames + // "J2000" and "ICRF" are treated the same, with "J2000" being the compatibility label. + reset_c(); + spkez_c(target, tdb2000, "J2000", "NONE", center, pv, <); + err = return_c(); + reset_c(); + + if(err) { + SpiceInt alt = novas_to_dexxx_planet(body); + if(alt != target) { + // Try with DExxx ID (barycenter vs planet center) + spkez_c(alt, tdb2000, "J2000", "NONE", center, pv, <); + err = return_c(); + reset_c(); + } + } + + mutex_unlock(sem); + + if(err) + return novas_error(3, EAGAIN, fn, "spkez_c() error (NOVAS ID=%d)", body); + + for(i = 3; --i >= 0;) { + if(position) position[i] = pv[i] * NORM_POS; + if(velocity) velocity[i] = pv[3 + i] * NORM_VEL; + } + + return 0; +} + +/** + * Provides an interface between the NAIF CSPICE library and NOVAS-C for regular (reduced) + * precision applications. + * + * This call is generally thread safe (notwithstanding outside access to the ephemeris files), + * even if CSPICE itself may not be. All ephemeris access will be mutexed to ensure sequential + * access under the hood. + * + * The call will use whatever ephemeris (SPK) files were loaded by the CSPICE library prior + * to the call (see furnsh_c() function) + * + * REFERENCES: + *
    + *
  1. NAIF CSPICE: https://naif.jpl.nasa.gov/naif/toolkit.html
    2. + *
  2. Kaplan, G. H. "NOVAS: Naval Observatory Vector Astrometry + * Subroutines"; USNO internal document dated 20 Oct 1988; + * revised 15 Mar 1990.
    3. + *
+ * + * @param jd_tdb [day] Two-element array containing the Julian date, which may be + * split any way (although the first element is usually the + * "integer" part, and the second element is the "fractional" + * part). Julian date is on the TDB or "T_eph" time scale. + * @param body Major planet number (or that for Sun, Moon, or Solar-system barycenter) + * @param origin NOVAS_BARYCENTER (0) or NOVAS_HELIOCENTER (1), or 2 for Earth geocenter + * -- relative to which to report positions and velocities. + * @param[out] position [AU] Position vector of 'body' at jd_tdb; equatorial rectangular + * coordinates in AU referred to the ICRS. + * @param[out] velocity [AU/day] Velocity vector of 'body' at jd_tdb; equatorial rectangular + * system referred to the ICRS, in AU/day. + * @return 0 if successful, or else an error code of solarsystem(). + * + * @sa planet_eph_manager_hp() + * @sa planet_ephem_provider() + * @sa ephem_open() + * @sa set_planet_provider() + * @sa solarsystem() + * + * @author Attila Kovacs + * @since 1.2 + */ +static short planet_cspice(double jd_tdb, enum novas_planet body, enum novas_origin origin, double *position, + double *velocity) { + const double tjd[2] = { jd_tdb, 0.0 }; + + prop_error("planet_cspice", planet_cspice_hp(tjd, body, origin, position, velocity), 0); + return 0; +} + +/** + * Generic ephemeris handling via the NAIF CSPICE library. This call is generally thread safe + * (notwithstanding outside access to the ephemeris files), even if CSPICE itself may not be. + * The ephemeris access will be mutexed to ensure sequential access under the hood. + * + * The call will use whatever ephemeris (SPK) files were loaded by the CSPICE library prior + * to the call (see furnsh_c() function) + * + * @param name The name of the solar-system body. It is important only if the 'id' is + * -1. + * @param id The NAIF ID number of the solar-system body for which the position in + * desired, or -1 if the 'name' should be used instead to identify the + * object. + * @param jd_tdb_high [day] The high-order part of Barycentric Dynamical Time (TDB) based + * Julian date for which to find the position and velocity. Typically + * this may be the integer part of the Julian date for high-precision + * calculations, or else the entire Julian date for reduced precision. + * @param jd_tdb_low [day] The low-order part of Barycentric Dynamical Time (TDB) based + * Julian date for which to find the position and velocity. Typically + * this may be the fractional part of the Julian date for high-precision + * calculations, or else 0.0 if the date is defined entirely by the + * high-order component for reduced precision. + * @param[out] origin Set to NOVAS_BARYCENTER or NOVAS_HELIOCENTER to indicate relative to + * which the ephemeris positions/velocities are reported. + * @param[out] pos [AU] position 3-vector to populate with rectangular equatorial + * coordinates in AU. It may be NULL if position is not required. + * @param[out] vel [AU/day] velocity 3-vector to populate in rectangular equatorial + * coordinates in AU/day. It may be NULL if velocities are not required. + * @return 0 if successful, -1 if any of the pointer arguments are NULL, or some + * non-zero value if the was an error s.t. the position and velocity + * vector should not be used. + * + * @sa novas_cspice_use_ids() + * @sa set_ephem_provider() + * @sa ephemeris() + * @sa NOVAS_EPHEM_OBJECT + * + *@author Attila Kovacs + * @since 1.2 + */ +static int novas_cspice(const char *name, long id, double jd_tdb_high, double jd_tdb_low, enum novas_origin *origin, double *pos, double *vel) { + static const char *fn = "novas_cspice"; + + SpiceDouble pv[6]; + SpiceDouble lt; + SpiceInt target, center; + double tdb2000; + int i, err; + + if(id == -1) { + // Lookup by name... + SpiceChar spiceName[SPICE_WORD_SIZE] = {}; + SpiceInt spiceCode = 0; + SpiceBoolean spiceFound = 0; + + if(!name) + return novas_error(-1, EINVAL, fn, "id=-1 and name is NULL"); + + if(!name[0]) + return novas_error(-1, EINVAL, fn, "id=-1 and name is empty"); + + // Use name to get NAIF ID. + strncpy(spiceName, name, sizeof(spiceName) - 1); + + reset_c(); + bodn2c_c(spiceName, &spiceCode, &spiceFound); + err = return_c(); + reset_c(); + + if(!spiceFound) + return novas_error(1, EINVAL, fn, "CSPICE could not find a NAIF ID for '%s'", name); + + if(err) + return novas_error(1, EINVAL, fn, "CSPICE name lookup error for '%s'", name); + + id = spiceCode; + } + + target = id; + + // Always return positions and velocities w.r.t. the SSB + if(origin) + *origin = NOVAS_SSB; + + center = NAIF_SSB; + + tdb2000 = (jd_tdb_high + jd_tdb_low - NOVAS_JD_J2000) * 86400.0; + + prop_error(fn, mutex_lock(), 0); + + // See https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/frames.html#Reference%20Frames + // "J2000" and "ICRF" are treated the same, with "J2000" being the compatibility label. + reset_c(); + spkez_c(target, tdb2000, "J2000", "NONE", center, pv, <); + err = return_c(); + reset_c(); + + mutex_unlock(); + + if(err) + return novas_error(3, EAGAIN, fn, "spkez_c() failure (name='%s', NAIF=%ld)", name ? name : "", id); + + for(i = 3; --i >= 0;) { + if(pos) pos[i] = pv[i] * NORM_POS; + if(vel) vel[i] = pv[3 + i] * NORM_VEL; + } + + return 0; +} + +/** + * Sets a ephemeris provider for NOVAS_EPHEM_OBJECT types using the NAIF CSPICE library. + * + * CSPICE is configured to suppress error messages and to not exit on errors, since we will check + * errors and handle them ourselves. You can adjust the behavior after this call with the + * CSPICE errprt_c() and erract_c() functions, respectively. + * + * @return 0 + * + * @sa novas_use_cspice_planets() + * @sa novas_use_cspice() + * @sa set_ephem_provider() + * + * @author Attila Kovacs + * @since 1.2 + */ +int novas_use_cspice_ephem() { + errprt_c("SET", 100, "NONE"); // Suppress CSPICE error messages + erract_c("SET", 0, "RETURN"); // Do not exit in case of CSPICE errors. + set_ephem_provider(novas_cspice); + return 0; +} + +/** + * Sets CSPICE as the ephemeris provider for the major planets (and Sun, Moon, and SSB) using the + * NAIF CSPICE library. + * + * CSPICE is configured to suppress error messages and to not exit on errors, since we will check + * errors and handle them ourselves. You can adjust the behavior after this call with the + * CSPICE errprt_c() and erract_c() functions, respectively. + * + * @return 0 + * + * @sa novas_use_cspice_ephem() + * @sa novas_use_cspice() + * @sa set_planet_provider() + * @sa set_planet_provider_hp() + * + * @author Attila Kovacs + * @since 1.2 + */ +int novas_use_cspice_planets() { + errprt_c("SET", 100, "NONE"); // Suppress CSPICE error messages + erract_c("SET", 0, "RETURN"); // Do not exit in case of CSPICE errors. + set_planet_provider_hp(planet_cspice_hp); + set_planet_provider(planet_cspice); + return 0; +} + +/** + * Sets CSPICE as the default ephemeris provider for all types of Solar-system objects (both NOVAS_PLANET and + * NOVAS_EPHEM_OBJECT types). + * + * CSPICE is configured to suppress error messages and to not exit on errors, since we will check + * errors and handle them ourselves. You can adjust the behavior after this call with the + * CSPICE errprt_c() and erract_c() functions, respectively. + * + * @return 0 + * + * @sa novas_use_cspice_planets() + * @sa novas_use_cspice_ephem() + * + * @author Attila Kovacs + * @since 1.2 + */ +int novas_use_cspice() { + novas_use_cspice_planets(); + novas_use_cspice_ephem(); + return 0; +} + + + diff --git a/test/Makefile b/test/Makefile index 65856dbb..cbb9c7e6 100644 --- a/test/Makefile +++ b/test/Makefile @@ -24,6 +24,14 @@ ifeq ($(CALCEPH_SUPPORT), 1) COVFILES += solsys-calceph.c.gcov endif +ifeq ($(CSPICE_SUPPORT), 1) + CPPFLAGS += -DUSE_CSPICE=1 + LDFLAGS += -lcspice + OBJECTS += solsys-cspice.o + TESTS += test-cspice + COVFILES += solsys-cspice.c.gcov +endif + .PHONY: run run: clean-cov clean-data $(TESTS) data ../cio_ra.bin bad-data ./test-compat @@ -37,6 +45,9 @@ run: clean-cov clean-data $(TESTS) data ../cio_ra.bin bad-data ifeq ($(CALCEPH_SUPPORT), 1) ./test-calceph ephem endif +ifeq ($(CSPICE_SUPPORT), 1) + ./test-cspice ephem +endif cov: cov.info genhtml $< -o $@ @@ -93,6 +104,8 @@ test-%.o: src/test-%.c | Makefile test-calceph: LDFLAGS += -lcalceph +test-cspice: LDFLAGS += -lcspice + .PHONY: clean-test clean-test: rm -f test-* diff --git a/test/src/test-calceph.c b/test/src/test-calceph.c index 96217700..a32791fe 100644 --- a/test/src/test-calceph.c +++ b/test/src/test-calceph.c @@ -76,7 +76,8 @@ static int test_calceph() { sprintf(filename, "%s/" PLANET_EPH, prefix); eph = calceph_open(filename); - if(novas_use_calceph(eph)) return 1; + + if(check("calceph:use", 0, novas_use_calceph(eph))) return 1; if(!is_ok("calceph:earth", ephemeris(jd2, &earth, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos, vel))) return 1; earth_sun_calc(jd, NOVAS_EARTH, NOVAS_BARYCENTER, pos0, vel0); @@ -109,11 +110,11 @@ static int test_calceph_planet() { sprintf(filename, "%s/" MARS_EPH, prefix); eph = calceph_open(filename); - if(novas_use_calceph(eph)) return 1; + if(check("calceph_planet:use", 0, novas_use_calceph(eph))) return 1; sprintf(filename, "%s/" PLANET_EPH, prefix); eph = calceph_open(filename); - if(novas_use_calceph_planets(eph)) return 1; + if(check("calceph_planet:use_planets", 0, novas_use_calceph_planets(eph))) return 1; if(!is_ok("calceph_planet:ssb", ephemeris(jd2, &ssb, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos, vel))) return 1; if(!is_ok("calceph_planet:ssb:pos", check_equal_pos(pos, pos0, 1e-5))) return 1; @@ -179,6 +180,7 @@ static int test_errors() { double pos[3], vel[3]; double jd = NOVAS_JD_J2000; double jd2[2] = { jd, 0.0 }; + int n = 0; object earth, phobos; novas_planet_provider_hp pl = get_planet_provider_hp(); @@ -187,28 +189,21 @@ static int test_errors() { make_planet(NOVAS_EARTH, &earth); make_ephem_object("Phobos", 401, &phobos); - if(check("errors:tdb", -1, pl(NULL, NOVAS_MARS, NOVAS_BARYCENTER, pos, vel))) return 1; - - if(check("errors:planet:number:-1", 1, pl(jd2, -1, NOVAS_BARYCENTER, pos, vel))) return 1; - - if(check("errors:planet:number:hi", 1, pl(jd2, NOVAS_PLANETS, NOVAS_BARYCENTER, pos, vel))) return 1; - - if(check("errors:planet:origin", 2, pl(jd2, NOVAS_MARS, -1, pos, vel))) return 1; + if(check("errors:tdb", -1, pl(NULL, NOVAS_MARS, NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:planet:number:-1", 1, pl(jd2, -1, NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:planet:number:hi", 1, pl(jd2, NOVAS_PLANETS, NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:planet:origin", 2, pl(jd2, NOVAS_MARS, -1, pos, vel))) n++; calceph_seterrorhandler(3, dummy_error_handler); - - if(check("errors:body:name:NULL", -1, eph(NULL, -1, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) return 1; - - if(check("errors:body:name:empty", -1, eph("", -1, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) return 1; - - if(check("errors:body:name:nomatch", 1, eph("blah", -1, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) return 1; + if(check("errors:body:name:NULL", -1, eph(NULL, -1, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:body:name:empty", -1, eph("", -1, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:body:name:nomatch", 1, eph("blah", -1, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) n++; jd2[0] = -999999.0; - if(check("errors:planet:time", 3, pl(jd2, NOVAS_MARS, NOVAS_BARYCENTER, pos, vel))) return 1; + if(check("errors:planet:time", 3, pl(jd2, NOVAS_MARS, NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:body:time", 3, eph("phobos", 401, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) n++; - if(check("errors:body:time", 3, eph("phobos", 401, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) return 1; - - return 0; + return n; } diff --git a/test/src/test-cspice.c b/test/src/test-cspice.c new file mode 100644 index 00000000..4dfe2937 --- /dev/null +++ b/test/src/test-cspice.c @@ -0,0 +1,208 @@ +/** + * @file + * + * @date Created on Feb 18, 2024 + * @author Attila Kovacs + */ + +#define _XOPEN_SOURCE 500 /// strdup() + + +#include +#include +#include +#include +#include + +#define __NOVAS_INTERNAL_API__ ///< Use definitions meant for internal use by SuperNOVAS only +#define USE_CSPICE 1 ///< Use cspice routines +#include "novas.h" +#include "solarsystem.h" + +#include "cspice/SpiceUsr.h" +#include "cspice/SpiceZpr.h" // for reset_c + +#define PLANET_EPH "de440s-j2000.bsp" +#define MARS_EPH "mar097-j2000.bsp" + +static const char *prefix; + +static int usage() { + fprintf(stderr, " Syntax: test-cspice \n\n"); + fprintf(stderr, " Path to de440s.bsp and mar097.bsp containing J2000 data.\n\n"); + exit(1); +} + +static int check_equal_pos(const double *posa, const double *posb, double tol) { + int i; + + tol = fabs(tol); + if(tol < 1e-30) tol = 1e-30; + + for(i = 0; i < 3; i++) { + if(fabs(posa[i] - posb[i]) <= tol) continue; + if(isnan(posa[i]) && isnan(posb[i])) continue; + + fprintf(stderr, " A[%d] = %.9g vs B[%d] = %.9g\n", i, posa[i], i, posb[i]); + return i + 1; + } + + return 0; +} + +static int is_ok(const char *func, int error) { + if(error) fprintf(stderr, "ERROR %d! %s\n", error, func); + return !error; +} + +static int check(const char *func, int exp, int error) { + if(error != exp) { + fprintf(stderr, "ERROR! %s: expected %d, got %d\n", func, exp, error); + return 1; + } + return 0; +} + + +static int test_cspice() { + double pos[3], vel[3], pos0[3], vel0[3]; + double jd = NOVAS_JD_J2000; + double jd2[2] = { jd, 0.0 }; + + char filename[1024]; + object earth, mars; + + make_planet(NOVAS_EARTH, &earth); + make_planet(NOVAS_MARS, &mars); + + if(!is_ok("use_cspice", novas_use_cspice() != 0)) return 1; + + if(!is_ok("cspice:earth", ephemeris(jd2, &earth, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos, vel))) return 1; + earth_sun_calc(jd, NOVAS_EARTH, NOVAS_BARYCENTER, pos0, vel0); + + if(!is_ok("cspice:earth:pos", check_equal_pos(pos, pos0, 1e-5))) return 1; + if(!is_ok("cspice:earth:vel", check_equal_pos(vel, vel0, 1e-5))) return 1; + + if(!is_ok("cspice_planet:mars", ephemeris(jd2, &mars, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos0, vel0))) return 1; + + + + return 0; +} + +static int test_cspice_planet() { + double pos[3], vel[3], pos0[3] = {}, vel0[3] = {}; + double jd = NOVAS_JD_J2000; + double jd2[2] = { jd, 0.0 }; + + char filename[1024]; + object ssb, sun, earth, moon, mars, phobos; + + make_planet(NOVAS_SSB, &ssb); + make_planet(NOVAS_SUN, &sun); + make_planet(NOVAS_EARTH, &earth); + make_planet(NOVAS_MOON, &moon); + make_planet(NOVAS_MARS, &mars); + make_ephem_object("Phobos", 401, &phobos); + + if(!is_ok("use_cspice", novas_use_cspice() != 0)) return 1; + + if(!is_ok("cspice_planet:ssb", ephemeris(jd2, &ssb, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos, vel))) return 1; + if(!is_ok("cspice_planet:ssb:pos", check_equal_pos(pos, pos0, 1e-5))) return 1; + + if(!is_ok("cspice_planet:sun_vs_sun", ephemeris(jd2, &sun, NOVAS_HELIOCENTER, NOVAS_REDUCED_ACCURACY, pos, vel))) return 1; + if(!is_ok("cspice_planet:sun_vs_sun:pos", check_equal_pos(pos, pos0, 1e-5))) return 1; + + if(!is_ok("cspice_planet:sun", ephemeris(jd2, &sun, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos, vel))) return 1; + earth_sun_calc(jd, NOVAS_SUN, NOVAS_BARYCENTER, pos0, vel0); + if(!is_ok("cspice_planet:sun:pos", check_equal_pos(pos, pos0, 1e-5))) return 1; + + if(!is_ok("cspice_planet:earth", ephemeris(jd2, &earth, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos, vel))) return 1; + earth_sun_calc(jd, NOVAS_EARTH, NOVAS_BARYCENTER, pos0, vel0); + + if(!is_ok("cspice_planet:earth:pos", check_equal_pos(pos, pos0, 1e-5))) return 1; + if(!is_ok("cspice_planet:earth:vel", check_equal_pos(vel, vel0, 1e-5))) return 1; + + if(!is_ok("cspice_planet:moon", ephemeris(jd2, &moon, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos, vel))) return 1; + earth_sun_calc(jd, NOVAS_MOON, NOVAS_BARYCENTER, pos0, vel0); + + if(!is_ok("cspice_planet:moon:pos", check_equal_pos(pos, pos0, 1e-2))) return 1; + if(!is_ok("cspice_planet:moon:vel", check_equal_pos(vel, vel0, 1e-3))) return 1; + + if(!is_ok("cspice_planet:mars", ephemeris(jd2, &mars, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos0, vel0))) return 1; + if(!is_ok("cspice_planet:phobos", ephemeris(jd2, &phobos, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos, vel))) return -1; + + if(!is_ok("cspice_planet:mars-phobos:pos", check_equal_pos(pos, pos0, 1e-4))) return 1; + + phobos.number = -1; + if(!is_ok("cspice_planet:phobos:byname", ephemeris(jd2, &phobos, NOVAS_BARYCENTER, NOVAS_REDUCED_ACCURACY, pos0, vel0))) return -1; + if(!is_ok("cspice_planet:phobos:match", check_equal_pos(pos, pos0, 1e-6))) return 1; + + + return 0; +} + + +static int test_errors() { + double pos[3], vel[3]; + double jd = NOVAS_JD_J2000; + double jd2[2] = { jd, 0.0 }; + int n = 0; + + object earth, phobos; + novas_planet_provider_hp pl = get_planet_provider_hp(); + novas_ephem_provider eph = get_ephem_provider(); + + make_planet(NOVAS_EARTH, &earth); + make_ephem_object("Phobos", 401, &phobos); + + if(check("errors:tdb", -1, pl(NULL, NOVAS_MARS, NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:planet:number:-1", 1, pl(jd2, -1, NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:planet:number:hi", 1, pl(jd2, NOVAS_PLANETS, NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:planet:origin", 2, pl(jd2, NOVAS_MARS, -1, pos, vel))) n++; + if(check("errors:body:name:NULL", -1, eph(NULL, -1, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:body:name:empty", -1, eph("", -1, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:body:name:nomatch", 1, eph("blah", -1, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) n++; + + jd2[0] = -999999.0; + + if(check("errors:planet:time", 3, pl(jd2, NOVAS_MARS, NOVAS_BARYCENTER, pos, vel))) n++; + if(check("errors:body:time", 3, eph("phobos", 401, jd2[0], jd2[1], NOVAS_BARYCENTER, pos, vel))) n++; + + return n; +} + + + +static void load_eph(const char *name) { + char filename[1024]; + + sprintf(filename, "%s/%s", prefix, name); + reset_c(); + furnsh_c(filename); + if(return_c()) { + fprintf(stderr, "ERROR! furnsh_c() failed for %s\n", name); + exit(1); + } +} + +int main(int argc, char *argv[]) { + int n = 0; + + if(argc < 2) usage(); + + prefix = strdup(argv[1]); + + load_eph(PLANET_EPH); + load_eph(MARS_EPH); + + enable_earth_sun_hp(1); + + if(test_cspice()) n++; + if(test_cspice_planet()) n++; + + novas_debug(NOVAS_DEBUG_OFF); + if(test_errors()) n++; + + return n; +}