Add support for orbital element based calculations

include/novas.h

@@ -244,11 +244,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
@@ -616,6 +622,35 @@ 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.
+ *
+ * @author Attila Kovacs
+ * @since 1.2
+ *
+ * @sa object
+ * @sa enum NOVAS_ORBITAL_OBJECT
+ */
+typedef struct {
+  enum novas_planet center;         ///< major body at orbital center (default is NOVAS_SUN).
+  double jd_tdb;                    ///< [day] Barycentri Dynamical Time (TDB) based Julian date of parameters.
+  double a;                         ///< [AU] semi-major axis
+  double e;                         ///< eccentricity
+  double omega;                     ///< [rad] argument of periapsis / perihelion
+  double Omega;                     ///< [rad] argument of ascending node
+  double i;                         ///< [rad] inclination
+  double M0;                        ///< [rad] mean anomaly at tjd
+  double n;                         ///< [rad/day] mean daily motion
+} novas_orbital_elements;
+
+/**
+ * Initializer for novas_orbital_elements for heliocentric orbits
+ *
+ * @author Attila Kovacs
+ * @since 1.2
+ */
+#define NOVAS_ORBIT_INIT = { NOVAS_SUN, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }
+
 /**
  * Celestial object of interest.
  *
@@ -627,7 +662,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;
 
 /**

include/solarsystem.h

@@ -340,6 +340,8 @@ long novas_to_naif_planet(enum novas_planet id);
 
 long novas_to_dexxx_planet(enum novas_planet id);
 
+int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orb, double *pos, double *vel);
+
 /// \cond PRIVATE
 
 

src/novas.c

@@ -40,6 +40,7 @@
 #endif
 
 #define C2                  (C * C)   ///< [m<sup>2</sup>/s<sup>2</sup>] Speed of light squared
+#define EPREC               1e-12     ///< Required precision for eccentric anomaly in orbital calculation
 
 // <---------- GLOBAL VARIABLES -------------->
 
@@ -5769,6 +5770,23 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig
       break;
     }
 
+    case NOVAS_ORBITAL_OBJECT: {
+      object center;
+      double pos0[3] = {0}, vel0[3] = {0};
+      int i;
+
+      prop_error(fn, make_planet(body->orbit.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, pos, vel), 0);
+
+      for(i = 3; --i >= 0; ) {
+        if(pos) pos[i] += pos0[i];
+        if(vel) vel[i] += vel0[i];
+      }
+
+      break;
+    }
+
     default:
       return novas_error(2, EINVAL, fn, "invalid Solar-system body type: %d", body->type);
 
@@ -5777,6 +5795,99 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig
   return 0;
 }
 
+/**
+ * Calculates a rectangular equatorial position and velocity vector for the given orbital elements for the
+ * specified time of observation.
+ *
+ * REFERENCES:
+ * <ol>
+ * <li>https://en.wikipedia.org/wiki/Orbital_elements</li>
+ * <li>https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf</li>
+ * <li>https://orbitalofficial.com/</li>
+ * </ol>
+ *
+ * @param jd_tdb    [day] Barycentric Dynamic Time (TDB) based Julian date
+ * @param orb       Orbital parameters
+ * @param[out] pos  [AU] Output position vector, or NULL if not required
+ * @param[out] vel  [AU/day] Output velocity vector, or NULL if not required
+ * @return          0 if successful, or else -1 if the orbital parameters are NULL
+ *                  or if the position and velocity output vectors are the same (errno set to EINVAL),
+ *                  or if the calculation did not converge (errno set to ECANCELED).
+ */
+int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orb, double *pos, double *vel) {
+  static const char *fn = "novas_orbit_posvel";
+
+  double dt, M, E, nu, r;
+  double cO, sO, ci, si, co, so;
+  double xx, yx, zx, xy, yy, zy;
+  int i = novas_inv_max_iter;
+
+  if(!orb)
+    return novas_error(-1, EINVAL, fn, "input orbital elements is NULL");
+
+  if(pos == vel)
+    return novas_error(-1, EINVAL, fn, "output pos = vel (@ %p)", pos);
+
+  dt = jd_tdb - orb->jd_tdb;
+  E = M = orb->M0 + orb->n * dt;
+
+  while(--i >= 0) {
+    double s = sin(E);
+    double c = cos(E);
+    double dE = (E - s - M) / (1.0 - c);
+    E -= dE;
+
+    if(fabs(dE) < EPREC) break;
+  }
+
+  if(i < 0)
+    return novas_error(-1, ECANCELED, fn, "Eccentric anomaly convergence failure");
+
+  nu = 2.0 * atan2(sqrt(1.0 + orb->e) * sin(0.5 * E), sqrt(1.0 - orb->e) * cos(0.5 * E));
+  r = orb->a * (1.0 - orb->e * E);
+
+  // pos = Rz(-Omega) . Rx(-i) . Rz(-omega) . orb
+  cO = cos(orb->Omega);
+  sO = sin(orb->Omega);
+  ci = cos(orb->i);
+  si = sin(orb->i);
+  co = cos(orb->omega);
+  so = sin(orb->omega);
+
+  // Rotation matrix
+  // See https://en.wikipedia.org/wiki/Euler_angles
+  // (note the Wikipedia has opposite sign convention for angles...)
+  xx = cO * co - sO * ci * so;
+  yx = sO * co + cO * ci * so;
+  zx = si * so;
+  xy = -cO * so - sO * ci * co;
+  yy = -sO * so + cO * ci * co;
+  zy = si * co;
+
+  if(pos) {
+    double x = cos(nu);
+    double y = sin(nu);
+
+    // Perform rotation
+    pos[0] = xx * x + xy * y;
+    pos[1] = yx * x + yy * y;
+    pos[2] = zx * x + zy * y;
+  }
+
+  if(vel) {
+    double vt = orb->n / r;
+    double x = -vt * sin(E);
+    double y = vt * sqrt(1.0 - orb->e * orb->e) * cos(E);
+
+    // Perform rotation
+    vel[0] = xx * x + xy * y;
+    vel[1] = yx * x + yy * y;
+    vel[2] = zx * x + zy * y;
+  }
+
+  return 0;
+}
+
 /**
  * Convert Hipparcos catalog data at epoch J1991.25 to epoch J2000.0, for use within NOVAS.
  * To be used only for Hipparcos or Tycho stars with linear space motion.  Both input and

src/super.c

@@ -1329,3 +1329,4 @@ enum novas_planet novas_planet_for_name(const char *name) {
 
   return novas_error(-1, EINVAL, fn, "No match for name: '%s'", name);
 }
+