Orbital elements based calculations (#96)

* Add support for orbital element based calculations

* Many fixes and some basic testing

* More complex orbital systems, such as planet equators, local Laplace planes...

* doc updates

* Fixes, tweaks, and satellte orbit test

* Change orbital angle units to deg (NOVAS convention), and add precession terms

* More tests and fixes

CHANGELOG.md

@@ -56,15 +56,26 @@ Changes expected for the next feature release, expected around 1 February 2025.
    You can also set other standard directory variables, as prescribed in the 
    [GNU standard](https://www.gnu.org/prep/standards/html_node/Directory-Variables.html) to further customize the
    install locations.
-
- - SuperNOVAS headers now include each other as system-headers, not local headers. This is unlikely to affect anything
-   really but it is more proper for an installation of the library, and works with our own `Makefile` too.
+
+ - #95: Added support for using orbital elements. `object.type` can now be set to `NOVAS_ORBITAL_OBJECT`, whose orbit
+   can be defined by the set of `novas_orbital_elements`, relative to a `novas_orbital_system`. You can initialize an 
+   `object` with a set of orbital elements using `make_orbital_object()`, and for planetary satellite orbits you might
+   use `novas_set_orbital_pole()`. For orbital objects, `ephemeris()` will call on the new `novas_orbit_posvel()` to 
+   calculate positions. While orbital elements do not always yield precise positions, they can for shorter periods, 
+   provided that the orbital elements are up-to-date. For example, the Minor Planer Center (MPC) publishes accurate 
+   orbital elements for all known asteroids and comets regularly. For newly discovered objects, this may be the only 
+   and/or most accurate information available anywhere.
 
  - #97: Added `NOVAS_EMB` (Earth-Moon Barycenter) and `NOVAS_PLUTO_BARYCENTER` to `enum novas_planets` to distinguish
    from the planet center in calculations.
 
+ - SuperNOVAS headers now include each other as system-headers, not local headers. This is unlikely to affect anything
+   really but it is more proper for an installation of the library, and works with our own `Makefile` too.
+
 ### Changed
 
+ - #96: Changed `object` structure to include `novas_orbital_elements` for `NOVAS_ORBITAL_OBJECT` types.
+
  - #97: Updated `NOVAS_PLANETS`, `NOVAS_PLANET_NAMES_INIT` and `NOVAS_RMASS_INIT` macros to include the added planet 
    constants.
 

README.md

@@ -596,6 +596,24 @@ more generic ephemeris handling via a user-provided `novas_ephem_provider`. E.g.
  make_ephem_object("Ceres", 2000001, &ceres);
 ```
 
+As of version 1.2 you can also define solar system sources with orbital elements (such as the most up-to-date ones 
+provided by the [Minor Planet Center](https://minorplanetcenter.net/data) for asteroids, comets, etc.):
+
+```c
+  object NEA;		// e.g. a Near-Earth Asteroid
+
+  // Fill in the orbital parameters (pay attention to units!)
+  novas_orbital_elements orbit = NOVAS_ORBIT_INIT;
+  orbit.a = ...;
+  ...
+
+  // Create an object for that orbit
+  make_orbital_object("NEAxxx", -1, &orbit, object);
+```
+
+Note, that even with orbital elements, you will, in general, require a planet calculator, to provide precise
+positions for the Sun or planet, around which the orbit is defined.
+
 Other than that, it's the same spiel as before, e.g.:
 
 ```c
@@ -894,6 +912,15 @@ 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.
 
+ - Added support for using orbital elements. `object.type` can now be set to `NOVAS_ORBITAL_OBJECT`, whose orbit
+   can be defined by the set of `novas_orbital_elements`, relative to a `novas_orbital_system`. You can initialize an 
+   `object` with a set of orbital elements using `make_orbital_object()`, and for planetary satellite orbits you might
+   use `novas_set_orbital_pole()`. For orbital objects, `ephemeris()` will call on the new `novas_orbit_posvel()` to 
+   calculate positions. While orbital elements do not always yield precise positions, they can for shorter periods, 
+   provided that the orbital elements are up-to-date. For example, the Minor Planer Center (MPC) publishes accurate 
+   orbital elements for all known asteroids and comets regularly. For newly discovered objects, this may be the only 
+   and/or most accurate information available anywhere.
+
  - Added `NOVAS_EMB` (Earth-Moon Barycenter) and `NOVAS_PLUTO_BARYCENTER` to `enum novas_planets` to distinguish
    from the corresponding planet centers in calculations.
 

include/novas.h

@@ -48,7 +48,6 @@
 // for compiling your code, which may have relied on these, you can add '-DCOMPAT=1' to the
 // compiler options
 //
-//
 #if COMPAT
 #  include <stdio.h>
 #  include <ctype.h>
@@ -246,11 +245,17 @@ enum novas_object_type {
 
   /// Any non-solar system object that may be handled via 'catalog' coordinates, such as a star
   /// or a quasar.
-  NOVAS_CATALOG_OBJECT
+  /// @sa cat_entry
+  NOVAS_CATALOG_OBJECT,
+
+  /// Any Solar-system body, whose position is determined by a set of orbital elements
+  /// @since 1.2
+  /// @sa novas_orbital_elements
+  NOVAS_ORBITAL_OBJECT
 };
 
 /// The number of object types distinguished by NOVAS.
-#define NOVAS_OBJECT_TYPES      (NOVAS_CATALOG_OBJECT + 1)
+#define NOVAS_OBJECT_TYPES      (NOVAS_ORBITAL_OBJECT + 1)
 
 /**
  * Enumeration for the 'major planet' numbers in NOVAS to use as the solar-system body number whenever
@@ -371,6 +376,10 @@ enum novas_equator_type {
   NOVAS_GCRS_EQUATOR      ///< Geocentric Celestial Reference System (GCRS).
 };
 
+/// The number of equator types we define.
+/// @since 1.2
+#define NOVAS_EQUATOR_TYPES (NOVAS_GCRS_EQUATOR + 1)
+
 /**
  * Constants that determine the type of dynamical system type for gcrs2equ()
  */
@@ -577,6 +586,16 @@ enum novas_nutation_direction {
   NUTATE_MEAN_TO_TRUE
 };
 
+/**
+ * The plane in which values, such as orbital parameyters are referenced.
+ *
+ * @since 1.2
+ */
+enum novas_reference_plane {
+  NOVAS_ECLIPTIC_PLANE = 0,     ///< the plane of the ecliptic
+  NOVAS_EQUATORIAL_PLANE        ///< The plane of the equator
+};
+
 /**
  * Fundamental Delaunay arguments of the Sun and Moon, from Simon section 3.4(b.3),
  * precession = 5028.8200 arcsec/cy)
@@ -622,6 +641,96 @@ typedef struct {
  */
 #define CAT_ENTRY_INIT { {'\0'}, {'\0'}, 0L, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }
 
+/**
+ * Specification of an orbital system, in which orbital elements are defined. Systems can be defined around
+ * all major planets and barycenters (and Sun, Moon, SSB..). They may be referenced to the GCRS, mean, or true equator
+ * or ecliptic of date, or to a plane that is tilted relative to that.
+ *
+ * For example, The Minor Planet Center (MPC) publishes up-to-date orbital elements for asteroids and comets,
+ * which are heliocentric and referenced to the GCRS ecliptic. Hence 'center' for these is `NOVAS_SUN`, the `plane`
+ * is `NOVAS_ECLIPTIC_PLANE` and the `type` is `NOVAS_GCRS_EQUATOR`.
+ *
+ * The orbits of planetary satellites may be parametrized in their local Laplace planes, which are typically close
+ * to the host planet's equatorial planes. You can, for example, obtain the RA/Dec orientation of the planetary
+ * North poles of planets from JPL Horizons, and use them as a proxy for the Laplace planes for their satellite orbits.
+ * In this case you would set the `center` to the host planet (e.g. `NOVAS_SATURN`), the reference plane to
+ * `NOVAS_EQUATORIAL_PLANE` and the `type` to `NOVAS_GCRS_EQUATOR` (since the plane is defined by the North pole
+ * orientation in GCRS equatorial RA/Dec). The obliquity is then 90&deg; - Dec<sub>pole</sub> (in radians), and `phi`
+ * is RA<sub>pole</sub> (in radians).
+ *
+ * @author Attila Kovacs
+ * @since 1.2
+ *
+ * @sa novas_orbital_elements
+ */
+typedef struct {
+  enum novas_planet center;          ///< major planet or barycenter at the center of the orbit.
+  enum novas_reference_plane plane;  ///< reference plane NOVAS_ECLIPTIC_PLANE or NOVAS_EQUATORIAL_PLANE
+  enum novas_equator_type type;      ///< the type of equator in which orientation is referenced (NOVAS_TRUE_EQUATOR,
+                                     ///< NOVAS_MEAN_EQUATOR, or NOVAS_GCRS_EQUATOR).
+  double obl;                        ///< [rad] relative obliquity of orbital reference plane
+                                     ///<       (e.g. 90&deg; - &delta;<sub>pole</sub>)
+  double Omega;                      ///< [rad] relative argument of ascending node of the orbital reference plane
+                                     ///<       (e.g. &alpha;<sub>pole</sub> + 90&deg;)
+} novas_orbital_system;
+
+
+/// Default orbital system initializer for heliocentric GCRS ecliptic orbits.
+/// @since 1.2
+#define NOVAS_ORBITAL_SYSTEM_INIT { NOVAS_SUN, NOVAS_ECLIPTIC_PLANE, NOVAS_GCRS_EQUATOR, 0.0, 0.0 }
+
+/**
+ * Keplerian orbital elements for `NOVAS_ORBITAL_OBJECT` type. Orbital elements can be used to provide
+ * approximate positions for various Solar-system bodies. JPL publishes orbital elements (and their evolution)
+ * for the major planets and their satellites. However, these are suitable only for very approximate
+ * calculations, with up to degree scale errors for the gas giants for the time range between 1850 AD and
+ * 2050 AD. Accurate positions and velocities for planets and their satellites should generally require the
+ * use of precise ephemeris data instead, such as obtained from the JPL Horizons system.
+ *
+ * Orbital elements describe motion from a purely Keplerian perspective. However, for short periods, for
+ * which the perturbing bodies can be ignored, this description can be quite accurate provided that an
+ * up-to-date set of values are used. The Minor Planet Center (MPC) thus regularly publishes orbital
+ * elements for all known asteroids and comets. For such objects, orbital elements can offer precise, and
+ * the most up-to-date positions and velocities.
+ *
+ * REFERENCES:
+ * <ol>
+ * <li>Up-to-date orbital elements for asteroids, comets, etc from the Minor Planet Center (MPC):
+ * https://minorplanetcenter.net/data</li>
+ * <li>Mean elements for planetary satellites from JPL Horizons: https://ssd.jpl.nasa.gov/sats/elem/</li>
+ * <li>Low accuracy mean elements for planets from JPL Horizons:
+ * https://ssd.jpl.nasa.gov/planets/approx_pos.html</li>
+ * </ol>
+ *
+ * @author Attila Kovacs
+ * @since 1.2
+ *
+ * @sa object
+ * @sa enum NOVAS_ORBITAL_OBJECT
+ */
+typedef struct {
+  novas_orbital_system system;      ///< orbital reference system assumed for the parametrization
+  double jd_tdb;                    ///< [day] Barycentri Dynamical Time (TDB) based Julian date of the parameters.
+  double a;                         ///< [AU] semi-major axis
+  double e;                         ///< eccentricity
+  double omega;                     ///< [deg] argument of periapsis / perihelion, at the reference time
+  double Omega;                     ///< [deg] argument of ascending node on the reference plane, at the reference time
+  double i;                         ///< [deg] inclination of orbit to the reference plane
+  double M0;                        ///< [deg] mean anomaly at the reference time
+  double n;                         ///< [deg/day] mean daily motion, i.e. (_GM_/_a_<sup>3</sup>)<sup>1/2</sup> for the central body,
+                                    ///< or 360/T, where T is orbital period in days.
+  double apsis_period;              ///< [day] Precession period of the apsis, if known.
+  double node_period;               ///< [day] Precession period of the ascending node, if known.
+} novas_orbital_elements;
+
+/**
+ * Initializer for novas_orbital_elements for heliocentric orbits using GCRS ecliptic pqrametrization.
+ *
+ * @author Attila Kovacs
+ * @since 1.2
+ */
+#define NOVAS_ORBIT_INIT { NOVAS_ORBITAL_SYSTEM_INIT, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }
+
 /**
  * Celestial object of interest.
  *
@@ -633,7 +742,8 @@ typedef struct {
   enum novas_object_type type;    ///< NOVAS object type
   long number;                    ///< enum novas_planet, or minor planet ID (e.g. NAIF), or star catalog ID.
   char name[SIZE_OF_OBJ_NAME];    ///< name of the object (0-terminated)
-  cat_entry star;                 ///< basic astrometric data for any 'catalog' object.
+  cat_entry star;                 ///< basic astrometric data for NOVAS_CATALOG_OBJECT type.
+  novas_orbital_elements orbit;   ///< orbital data for NOVAS_ORBITAL_OBJECT TYPE
 } object;
 
 /**
@@ -1305,6 +1415,12 @@ double novas_z_inv(double z);
 
 enum novas_planet novas_planet_for_name(const char *name);
 
+int novas_set_orbital_pole(double ra, double dec, novas_orbital_system *sys);
+
+int make_orbital_object(const char *name, long num, const novas_orbital_elements *orbit, object *body);
+
+int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orb, enum novas_accuracy accuracy, double *pos, double *vel);
+
 
 // <================= END of SuperNOVAS API =====================>
 

src/frames.c

@@ -315,8 +315,8 @@ static int is_frame_initialized(const novas_frame *frame) {
 int novas_make_frame(enum novas_accuracy accuracy, const observer *obs, const novas_timespec *time, double dx, double dy,
         novas_frame *frame) {
   static const char *fn = "novas_make_frame";
-  static const object earth = { NOVAS_PLANET, NOVAS_EARTH, "Earth", CAT_ENTRY_INIT };
-  static const object sun = { NOVAS_PLANET, NOVAS_SUN, "Sun", CAT_ENTRY_INIT };
+  static const object earth = { NOVAS_PLANET, NOVAS_EARTH, "Earth", CAT_ENTRY_INIT, NOVAS_ORBIT_INIT };
+  static const object sun = { NOVAS_PLANET, NOVAS_SUN, "Sun", CAT_ENTRY_INIT, NOVAS_ORBIT_INIT };
 
   double tdb2[2];
   double mobl, tobl, ee, dpsi, deps;