diff --git a/config.mk b/config.mk index 8324f941..5d8e21ba 100644 --- a/config.mk +++ b/config.mk @@ -159,7 +159,7 @@ ifeq ($(DEFAULT_SOLSYS), 3) CPPFLAGS += -DDEFAULT_SOLSYS=3 endif -SOURCES = $(SRC)/novas.c $(SRC)/nutation.c $(SRC)/super.c $(SRC)/timescale.c $(SRC)/frames.c $(SRC)/refract.c $(SRC)/naif.c +SOURCES = $(SRC)/novas.c $(SRC)/nutation.c $(SRC)/super.c $(SRC)/timescale.c $(SRC)/frames.c $(SRC)/refract.c $(SRC)/naif.c $(SRC)/orbit.c ifeq ($(BUILTIN_SOLSYS1), 1) SOURCES += $(SRC)/solsys1.c $(SRC)/eph_manager.c diff --git a/include/novas.h b/include/novas.h index 4df4082b..7d3fffa0 100644 --- a/include/novas.h +++ b/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; /** diff --git a/include/solarsystem.h b/include/solarsystem.h index bae376f6..7153ca3a 100644 --- a/include/solarsystem.h +++ b/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 diff --git a/src/novas.c b/src/novas.c index 0f90cd2d..339b0872 100644 --- a/src/novas.c +++ b/src/novas.c @@ -5769,6 +5769,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, ¢er), 0); + prop_error(fn, ephemeris(jd_tdb, ¢er, 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); diff --git a/src/orbit.c b/src/orbit.c new file mode 100644 index 00000000..008f4702 --- /dev/null +++ b/src/orbit.c @@ -0,0 +1,115 @@ +/** + * @file + * + * @date Created on Nov 13, 2024 + * @author Attila Kovacs + */ + +#include +#include +#include +#include + +#define __NOVAS_INTERNAL_API__ ///< Use definitions meant for internal use by SuperNOVAS only +#include "novas.h" + +/// \cond PRIVATE + +/// Required precision for the eccentric anomaly +#define EPREC 1e-12 + +/// \endcond + +/** + * Calculates a rectangular equatorial position and velocity vector for the given orbital elements for the + * specified time of observation. + * + * REFERENCES: + *
    + *
  1. https://en.wikipedia.org/wiki/Orbital_elements
  2. + *
  3. https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf
  4. + *
  5. https://orbitalofficial.com/
  6. + *
+ * + * @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"; + extern int novas_inv_max_iter; + + 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; +} diff --git a/src/super.c b/src/super.c index 9e9424d9..153482ea 100644 --- a/src/super.c +++ b/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); } +