diff --git a/include/novas.h b/include/novas.h index 4df4082b..566d1668 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 } + /** * 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/frames.c b/src/frames.c index c6cdd545..ee636073 100644 --- a/src/frames.c +++ b/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; diff --git a/src/novas.c b/src/novas.c index 0f90cd2d..948c9674 100644 --- a/src/novas.c +++ b/src/novas.c @@ -40,6 +40,7 @@ #endif #define C2 (C * C) ///< [m2/s2] Speed of light squared +#define EPREC 1e-12 ///< Required precision for eccentric anomaly in orbital calculation // <---------- GLOBAL VARIABLES --------------> @@ -1394,7 +1395,6 @@ short mean_star(double jd_tt, double tra, double tdec, enum novas_accuracy accur int obs_posvel(double jd_tdb, double ut1_to_tt, enum novas_accuracy accuracy, const observer *obs, const double *geo_pos, const double *geo_vel, double *pos, double *vel) { static const char *fn = "get_obs_posvel"; - static const cat_entry zero_star = CAT_ENTRY_INIT; if(!obs) return novas_error(-1, EINVAL, fn, "NULL observer parameter"); @@ -1415,7 +1415,7 @@ int obs_posvel(double jd_tdb, double ut1_to_tt, enum novas_accuracy accuracy, co if(!geo_pos || !geo_vel) { const double tdb2[2] = { jd_tdb }; - object earth = { NOVAS_PLANET, NOVAS_EARTH, "Earth", zero_star }; + object earth = { NOVAS_PLANET, NOVAS_EARTH, "Earth", CAT_ENTRY_INIT, NOVAS_ORBIT_INIT }; double gpos[3], gvel[3]; prop_error(fn, ephemeris(tdb2, &earth, NOVAS_BARYCENTER, accuracy, gpos, gvel), 0); if(pos) @@ -3412,7 +3412,6 @@ short geo_posvel(double jd_tt, double ut1_to_tt, enum novas_accuracy accuracy, c static THREAD_LOCAL double t_last = 0; static THREAD_LOCAL enum novas_accuracy acc_last = -1; static THREAD_LOCAL double gast; - static const cat_entry zero_star = CAT_ENTRY_INIT; double gmst, eqeq, pos1[3], vel1[3], jd_tdb, jd_ut1; @@ -3487,7 +3486,7 @@ short geo_posvel(double jd_tt, double ut1_to_tt, enum novas_accuracy accuracy, c } case NOVAS_SOLAR_SYSTEM_OBSERVER: { // Observer in Solar orbit - const object earth = { NOVAS_PLANET, NOVAS_EARTH, "Earth", zero_star }; + const object earth = { NOVAS_PLANET, NOVAS_EARTH, "Earth", CAT_ENTRY_INIT, NOVAS_ORBIT_INIT }; const double tdb[2] = { jd_tdb, 0.0 }; int i; @@ -5769,6 +5768,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); @@ -5777,6 +5793,100 @@ 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: + *
    + *
  1. https://en.wikipedia.org/wiki/Orbital_elements
    2. + *
  2. https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf
    3. + *
  3. https://orbitalofficial.com/
    4. + *
+ * + * @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; + + // 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); + + 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 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); } +