More complex orbital systems, such as planet equators, local Laplace …

…planes...

include/novas.h

@@ -633,7 +633,61 @@ typedef struct {
 #define CAT_ENTRY_INIT { {'\0'}, {'\0'}, 0L, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }
 
 /**
- * Orbital elements for NOVAS_ORBITAL_OBJECT type.
+ * Specification of an orbital system, in which orbital elements are defined. Systems can be defined around
+ * all major planets (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 of 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 equator's orientation is defined in
+ * GCRS equatorial RA/Dec). The obliquity is then 90° - Dec (in radians), and `phi` is RA (in radians).
+ *
+ * @author Attila Kovacs
+ * @since 1.2
+ *
+ * @sa novas_orbital_elements
+ */
+typedef struct {
+  enum novas_planet center;          ///< major body at the center of the orbit (default is NOVAS_SUN).
+  enum novas_reference_plane plane;  ///< Reference plane NOVAS_ECLIPTIC_PLANE (default) or NOVAS_EQUATORIAL_PLANE
+  enum novas_equator_type type;      ///< The type of equator in which orientation is referenced (true, mean, or GCRS [default]).
+  double obl;                        ///< [rad] obliquity to reference plane (e.g. 90&deg; - Dec<sub>pole</sub>)
+  double phi;                        ///< [rad] argument of ascending node (e.g. RA<sub>pole</sub>)
+} 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 }
+
+/**
+ * Orbital elements for `NOVAS_ORBITAL_OBJECT` type. Orbital elements can be used to provide approximate
+ * positions for various Solar-system bodies. JPL publishes orbital parameters 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 descrive 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
@@ -642,9 +696,7 @@ typedef struct {
  * @sa enum NOVAS_ORBITAL_OBJECT
  */
 typedef struct {
-  enum novas_planet center;         ///< major body at the center of the orbit (default is NOVAS_SUN).
-  enum novas_reference_plane plane; ///< Reference plane NOVAS_ECLIPTIC_PLANE (0) or NOVAS_EQUATORIAL_PLANE (1)
-  enum novas_equator_type sys;      ///< The type of equator to which coordinates are referenced (true, mean, or GCRS).
+  novas_orbital_system system;      ///< The orbital reference system assumed for the paramtetrization
   double jd_tdb;                    ///< [day] Barycentri Dynamical Time (TDB) based Julian date of parameters.
   double a;                         ///< [AU] semi-major axis
   double e;                         ///< eccentricity
@@ -653,15 +705,15 @@ typedef struct {
   double i;                         ///< [rad] inclination
   double M0;                        ///< [rad] mean anomaly at tjd
   double n;                         ///< [rad/day] mean daily motion, i.e. (_GM_/_a_<sup>3</sup>)<sup>1/2</sup> for the central body
-} novas_orbital_elements;
+} novas_orbital_ents;
 
 /**
  * Initializer for novas_orbital_elements for heliocentric orbits
  *
  * @author Attila Kovacs
  * @since 1.2
  */
-#define NOVAS_ORBIT_INIT { NOVAS_SUN, NOVAS_ECLIPTIC_PLANE, NOVAS_GCRS_EQUATOR, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }
+#define NOVAS_ORBIT_INIT { NOVAS_ORBITAL_SYSTEM_INIT, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }
 
 /**
  * Celestial object of interest.

src/novas.c

@@ -5774,7 +5774,7 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig
       double pos0[3] = {0}, vel0[3] = {0};
       int i;
 
-      prop_error(fn, make_planet(body->orbit.center, &center), 0);
+      prop_error(fn, make_planet(body->orbit.system.center, &center), 0);
       prop_error(fn, ephemeris(jd_tdb, &center, origin, accuracy, pos0, vel0), 0);
       prop_error(fn, novas_orbit_posvel(jd_tdb[0] + jd_tdb[1], &body->orbit, accuracy, pos, vel), 0);
 
@@ -5795,6 +5795,55 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig
   return 0;
 }
 
+
+/**
+ * Change xzy vectors to the new polar orientation. &theta, &phi define the orientation of the input pole in the output system.
+ *
+ *
+ * @param in        input 3-vector in the original system (pole = z)
+ * @param theta     [rad] polar angle of original pole in the new system
+ * @param phi       [rad] azimuthal angle of original pole in the new system
+ * @param[out] out  output 3-vector in the new (rotated) system. It may be the same vector as the input.
+ * @return          0 if successful, or else -1 (rrno set to EINVAL) if either vector parameters are NULL.
+ *
+ */
+static int change_pole(const double *in, double theta, double phi, double *out) {
+  static const char *fn = "novas_change_pole";
+  double x, y, z;
+
+  if(!in)
+    return novas_error(-1, EINVAL, fn, "input 3-vector is NULL");
+
+  if(!out)
+    return novas_error(-1, EINVAL, fn, "output 3-vector is NULL");
+
+  x = in[0];
+  y = in[1];
+  z = in[2];
+
+  // looking in Rz (phi) Rx (theta)
+  double ca = sin(phi);
+  double sa = cos(phi);
+  double cb = cos(theta);
+  double sb = sin(theta);
+
+  out[0] = ca * x - sa * cb * y + sa * sb * z;
+  out[1] = sb * x + ca * cb * y - ca * sb * z;
+  out[3] = sb * y + cb * z;
+
+  return 0;
+}
+
+/**
+ * Converts equatorial coordinates of a given type to GCRS equatorial coordinates
+ *
+ * @param jd_tdb    [day] Barycentric Dynamical Time (TDB) based Julian Date
+ * @param[in] in    input 3-vector
+ * @param sys       the type of equator assumed for the input (mean, true, or GCRS).
+ * @param[out] out  output 3-vector. It may be the same as the input.
+ * @return          0 if successful, or else -1 (errno set to EINVAL) if the 'sys'
+ *                  argument is invalid.
+ */
 static int equ2gcrs(double jd_tdb, const double *in, enum novas_equator_type sys, double *out) {
   switch(sys) {
     case NOVAS_GCRS_EQUATOR:
@@ -5810,6 +5859,29 @@ static int equ2gcrs(double jd_tdb, const double *in, enum novas_equator_type sys
   }
 }
 
+
+/**
+ * Convert coordinates in an orbital system to GCRS equatorial coordinates
+ *
+ * @param jd_tdb        [day] Barycentric Dynamic Time (TDB) based Julian Date
+ * @param sys           Orbital system specification
+ * @param accuracy      NOVAS_FULL_ACCURACY or NOVAS_REDUCED_ACCURACY
+ * @param[in, out] vec  Coordinates
+ * @return              0 if successful, or else an error from ecl2equ_vec().
+ *
+ */
+static int orbit2gcrs(double jd_tdb, const novas_orbital_system *sys, enum novas_accuracy accuracy, double *vec) {
+  if(sys->obl)
+    change_pole(vec, sys->obl, sys->phi, vec);
+
+  if(sys->plane == NOVAS_ECLIPTIC_PLANE)
+    prop_error("orbit2gcrs", ecl2equ_vec(jd_tdb, sys->type, accuracy, vec, vec), 0);
+
+  equ2gcrs(jd_tdb, vec, sys->type, vec);
+
+  return 0;
+}
+
 /**
  * Calculates a rectangular equatorial position and velocity vector for the given orbital elements for the
  * specified time of observation.
@@ -5900,10 +5972,7 @@ int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orbit, enum
     pos[1] = yx * x + yy * y;
     pos[2] = zx * x + zy * y;
 
-    if(orbit->plane == NOVAS_ECLIPTIC_PLANE)
-      prop_error(fn, ecl2equ_vec(jd_tdb, orbit->sys, accuracy, pos, pos), 0);
-
-    equ2gcrs(jd_tdb, pos, orbit->sys, pos);
+    prop_error(fn, orbit2gcrs(jd_tdb, &orbit->system, accuracy, pos), 0);
   }
 
   if(vel) {
@@ -5916,10 +5985,7 @@ int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orbit, enum
     vel[1] = yx * x + yy * y;
     vel[2] = zx * x + zy * y;
 
-    if(orbit->plane == NOVAS_ECLIPTIC_PLANE)
-      prop_error(fn, ecl2equ_vec(jd_tdb, orbit->sys, accuracy, vel, vel), 0);
-
-    equ2gcrs(jd_tdb, vel, orbit->sys, vel);
+    prop_error(fn, orbit2gcrs(jd_tdb, &orbit->system, accuracy, vel), 0);
   }
 
   return 0;

test/src/test-super.c

@@ -2209,7 +2209,7 @@ static int test_planet_for_name() {
 
 static int test_orbit_posvel() {
   object ceres = {};
-  novas_orbital_elements orbit;
+  novas_orbital_elements orbit = NOVAS_ORBIT_INIT;
   observer obs = {};
   sky_pos pos = {};
 
@@ -2221,10 +2221,6 @@ static int test_orbit_posvel() {
   double r = 3.323;               // AU
   int n = 0;
 
-  orbit.center = NOVAS_SUN;
-  orbit.plane = NOVAS_ECLIPTIC_PLANE;
-  orbit.sys = NOVAS_GCRS_EQUATOR;
-
   orbit.jd_tdb = 2460600.5;
   orbit.a = 2.7666197;
   orbit.e = 0.079184;