Many fixes and some basic testing

include/novas.h

@@ -577,6 +577,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)
@@ -632,15 +642,17 @@ typedef struct {
  * @sa enum NOVAS_ORBITAL_OBJECT
  */
 typedef struct {
-  enum novas_planet center;         ///< major body at orbital center (default is NOVAS_SUN).
+  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).
   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
+  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;
 
 /**
@@ -649,7 +661,7 @@ typedef struct {
  * @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 }
+#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 }
 
 /**
  * Celestial object of interest.
@@ -1335,6 +1347,10 @@ double novas_z_inv(double z);
 
 enum novas_planet novas_planet_for_name(const char *name);
 
+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 =====================>
 

include/solarsystem.h

@@ -340,8 +340,6 @@ 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

@@ -864,7 +864,7 @@ int radec_planet(double jd_tt, const object *ss_body, const observer *obs, doubl
   if(rv)
     *rv = NAN;
 
-  if(ss_body->type != NOVAS_PLANET && ss_body->type != NOVAS_EPHEM_OBJECT)
+  if(ss_body->type != NOVAS_PLANET && ss_body->type != NOVAS_EPHEM_OBJECT && ss_body->type != NOVAS_ORBITAL_OBJECT)
     return novas_error(-1, EINVAL, fn, "object is not solar-system type: type=%d", ss_body->type);
 
   prop_error(fn, place(jd_tt, ss_body, obs, ut1_to_tt, sys, accuracy, &output), 10);
@@ -4377,6 +4377,7 @@ double rad_vel2(const object *source, const double *pos_emit, const double *vel_
       /* fallthrough */
 
     case NOVAS_EPHEM_OBJECT:
+    case NOVAS_ORBITAL_OBJECT:
       // Solar potential at source (bodies strictly outside the Sun's volume)
       if(d_src_sun * AU > NOVAS_SOLAR_RADIUS)
         rel /= 1.0 - GS / (d_src_sun * AU) / C2;
@@ -5775,113 +5776,150 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig
 
       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);
+      prop_error(fn, novas_orbit_posvel(jd_tdb[0] + jd_tdb[1], &body->orbit, accuracy, pos, vel), 0);
 
       for(i = 3; --i >= 0; ) {
-        if(pos) pos[i] += pos0[i];
-        if(vel) vel[i] += vel0[i];
+        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);
-
   }
 
   return 0;
 }
 
+static int equ2gcrs(double jd_tdb, const double *in, enum novas_equator_type sys, double *out) {
+  switch(sys) {
+    case NOVAS_GCRS_EQUATOR:
+      memcpy(out, in, XYZ_VECTOR_SIZE);
+      return 0;
+    case NOVAS_TRUE_EQUATOR:
+      return tod_to_gcrs(jd_tdb, NOVAS_REDUCED_ACCURACY, in, out);
+    case NOVAS_MEAN_EQUATOR:
+      precession(jd_tdb, in, NOVAS_JD_J2000, out);
+      return j2000_to_gcrs(out, out);
+    default:
+      return novas_error(-1, EINVAL, "equ2gcrs", "invalid equator type: %d", sys);
+  }
+}
+
 /**
  * Calculates a rectangular equatorial position and velocity vector for the given orbital elements for the
  * specified time of observation.
  *
  * REFERENCES:
  * <ol>
+ * <li>https://ssd.jpl.nasa.gov/planets/approx_pos.html</li>
  * <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>
+ * <li>https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf</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
+ * @param orbit     Orbital parameters
+ * @param accuracy  NOVAS_FULL_ACCURACY (0) or NOVAS_REDUCED_ACCURACY (1).
+ * @param[out] pos  [AU] Output ICRS equatorial position vector, or NULL if not required
+ * @param[out] vel  [AU/day] Output ICRS equatorial 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).
+ *
+ * @sa ephemeris()
+ * @sa novas_geom_posvel()
+ * @sa place()
+ * @sa make_orbital_object()
+ *
+ * @author Attila Kovacs
+ * @since 1.2
  */
-int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orb, double *pos, double *vel) {
+int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orbit, enum novas_accuracy accuracy, double *pos, double *vel) {
   static const char *fn = "novas_orbit_posvel";
 
-  double dt, M, E, nu, r;
+  double 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)
+  if(!orbit)
     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;
+  E = M = orbit->M0 + orbit->n * (jd_tdb - orbit->jd_tdb);
 
   // Iteratively determine E, using Newton-Raphson method...
   while(--i >= 0) {
-    double s = sin(E);
-    double c = cos(E);
-    double dE = (E - s - M) / (1.0 - c);
+    double es = orbit->e * sin(E);
+    double ec = orbit->e * cos(E);
+    double dE = (E - es - M) / (1.0 - ec);
 
     E -= dE;
-    if(fabs(dE) < EPREC) break;
+    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);
+  nu = 2.0 * atan2(sqrt(1.0 + orbit->e) * sin(0.5 * E), sqrt(1.0 - orbit->e) * cos(0.5 * E));
+  r = orbit->a * (1.0 - orbit->e * cos(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);
+  cO = cos(orbit->Omega);
+  sO = sin(orbit->Omega);
+  ci = cos(orbit->i);
+  si = sin(orbit->i);
+  co = cos(orbit->omega);
+  so = sin(orbit->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);
+    double x = r * cos(nu);
+    double y = r * sin(nu);
 
     // Perform rotation
     pos[0] = xx * x + xy * y;
     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);
   }
 
   if(vel) {
-    double vt = orb->n / r;
-    double x = -vt * sin(E);
-    double y = vt * sqrt(1.0 - orb->e * orb->e) * cos(E);
+    double v = orbit->n * orbit->a * orbit->a / r;    // [AU/day]
+    double x = -v * sin(E);
+    double y = v * sqrt(1.0 - orbit->e * orbit->e) * cos(E);
 
     // Perform rotation
     vel[0] = xx * x + xy * y;
     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);
   }
 
   return 0;
@@ -6521,7 +6559,7 @@ void novas_case_sensitive(int value) {
  * however enable case-sensitive processing by calling novas_case_sensitive() before.
  *
  * @param type          The type of object. NOVAS_PLANET (0), NOVAS_EPHEM_OBJECT (1) or
- *                      NOVAS_CATALOG_OBJECT (2)
+ *                      NOVAS_CATALOG_OBJECT (2), or NOVAS_ORBITAL_OBJECT (3).
  * @param number        The novas ID number (for solar-system bodies only, otherwise ignored)
  * @param name          The name of the object (case insensitive). It should be shorter than
  *                      SIZE_OF_OBJ_NAME or else an error will be returned. The name is
@@ -6540,6 +6578,8 @@ void novas_case_sensitive(int value) {
  * @sa make_redshifted_object()
  * @sa make_planet()
  * @sa make_ephem_object()
+ * @sa make_orbital_object()
+ * @sa novas_geom_posvel()
  * @sa place()
  *
  */

src/super.c

@@ -990,6 +990,7 @@ int grav_undef(double jd_tdb, enum novas_accuracy accuracy, const double *pos_ap
  * @sa make_cat_entry()
  * @sa make_planet()
  * @sa make_ephem_object()
+ * @sa novas_geom_posvel()
  * @sa place()
  *
  * @since 1.1
@@ -1015,13 +1016,14 @@ int make_cat_object(const cat_entry *star, object *source) {
  *                      ephemeris provider used with NOVAS. (If the ephemeris provider is by name and not
  *                      ID number, then the number here is not important).
  * @param[out] body     Pointer to structure to populate.
- * @return              0 if successful, or else -1 if the 'planet' pointer is NULL or the name
+ * @return              0 if successful, or else -1 if the 'body' pointer is NULL or the name
  *                      is too long.
  *
  *
  * @sa set_ephem_provider()
  * @sa make_planet()
  * @sa make_cat_entry()
+ * @sa novas_geom_posvel()
  * @sa place()
  *
  * @since 1.0
@@ -1032,6 +1034,37 @@ int make_ephem_object(const char *name, long num, object *body) {
   return 0;
 }
 
+/**
+ * Sets a celestial object to be a Solar-system orbital body. Typically this would be used to define
+ * minor planets, asteroids, comets, or even planetary satellites.
+ *
+ * @param name          Name of object. It may be NULL if not relevant.
+ * @param num           Solar-system body ID number (e.g. NAIF). It is not required and can be set e.g. to
+ *                      -1 if not relevant to the caller.
+ * @param orbit         The orbital parameters to adopt. The data will be copied, not referenced.
+ * @param[out] body     Pointer to structure to populate.
+ * @return              0 if successful, or else -1 if the 'orbit' or 'body' pointer is NULL or the name
+ *                      is too long.
+ *
+ *
+ * @sa novas_orbit_posvel()
+ * @sa make_planet()
+ * @sa make_ephem_object()
+ * @sa novas_geom_posvel()
+ * @sa place()
+ *
+ * @since 1.2
+ * @author Attila Kovacs
+ */
+int make_orbital_object(const char *name, long num, const novas_orbital_elements *orbit, object *body) {
+  static const char *fn = "make_orbital_object";
+  if(!orbit)
+    return novas_error(-1, EINVAL, fn, "Input orbital elements is NULL");
+  prop_error(fn, (make_object(NOVAS_ORBITAL_OBJECT, num, name, NULL, body) ? -1 : 0), 0);
+  memcpy(&body->orbit, orbit, sizeof(novas_orbital_elements));
+  return 0;
+}
+
 /**
  * Populates a celestial object data structure with the parameters for a redhifted catalog
  * source, such as a distant quasar or galaxy. It is similar to `make_cat_object()` except

test/src/test-super.c

@@ -2207,6 +2207,47 @@ static int test_planet_for_name() {
 }
 
 
+static int test_orbit_posvel() {
+  object ceres = {};
+  novas_orbital_elements orbit;
+  observer obs = {};
+  sky_pos pos = {};
+
+  // This is for SMA, so topocentric with Earth rot...
+  double tjd = NOVAS_JD_MJD0 + 60627.796269;
+  double RA0 = 19.679722;         // 19:40:47.000
+  double DEC0 = -28.633014;       // -28:37:58.849
+  double rv0 = 21.138;            // km/s
+  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;
+  orbit.i = 10.5879 * DEGREE;
+  orbit.omega = 73.28579 * DEGREE;
+  orbit.Omega = 80.25414 * DEGREE;
+  orbit.M0 = 145.84905 * DEGREE;
+  orbit.n = 0.21418047 * DEGREE;
+
+  make_observer_at_geocenter(&obs);
+  make_orbital_object("Ceres", -1, &orbit, &ceres);
+
+  if(!is_ok("orbit_posvel", place(tjd, &ceres, &obs, ut12tt, NOVAS_TOD, NOVAS_REDUCED_ACCURACY, &pos))) return 1;
+
+  // TODO refine with HORIZONS data...
+  if(!is_equal("orbit_posvel:ra", pos.ra, RA0, 5e-5)) n++;
+  if(!is_equal("orbit_posvel:dec", pos.dec, DEC0, 1e-3)) n++;
+  if(!is_equal("orbit_posvel:dist", pos.dis, r, 2e-3)) n++;
+  if(!is_equal("orbit_posvel:vrad", pos.rv, rv0, 0.5)) n++;
+
+  return n;
+}
+
 int main(int argc, char *argv[]) {
   int n = 0;
 
@@ -2271,6 +2312,7 @@ int main(int argc, char *argv[]) {
   if(test_naif_to_novas_planet()) n++;
 
   if(test_planet_for_name()) n++;
+  if(test_orbit_posvel()) n++;
 
   n += test_dates();