From d8c58b79cafca929677fc11f72f70631299e9d3b Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Wed, 13 Nov 2024 16:10:49 +0100 Subject: [PATCH 1/7] Add support for orbital element based calculations --- include/novas.h | 42 +++++++++++++-- include/solarsystem.h | 2 + src/frames.c | 4 +- src/novas.c | 118 ++++++++++++++++++++++++++++++++++++++++-- src/super.c | 1 + 5 files changed, 158 insertions(+), 9 deletions(-) diff --git a/include/novas.h b/include/novas.h index 49533893..bf5e485d 100644 --- a/include/novas.h +++ b/include/novas.h @@ -246,11 +246,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 @@ -622,6 +628,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. * @@ -633,7 +668,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 f73a8cea..b7472398 100644 --- a/include/solarsystem.h +++ b/include/solarsystem.h @@ -341,6 +341,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 efea8e62..bf96c119 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); } + From 0065d43d24859f2cf65a0d6c7d58479ce53ca14f Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Wed, 13 Nov 2024 21:24:46 +0100 Subject: [PATCH 2/7] Many fixes and some basic testing --- include/novas.h | 22 +++++++-- include/solarsystem.h | 2 - src/novas.c | 104 +++++++++++++++++++++++++++++------------- src/super.c | 35 +++++++++++++- test/src/test-super.c | 42 +++++++++++++++++ 5 files changed, 167 insertions(+), 38 deletions(-) diff --git a/include/novas.h b/include/novas.h index bf5e485d..6a4a609d 100644 --- a/include/novas.h +++ b/include/novas.h @@ -583,6 +583,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) @@ -638,7 +648,9 @@ 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 @@ -646,7 +658,7 @@ typedef struct { 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_3)1/2 for the central body } novas_orbital_elements; /** @@ -655,7 +667,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. @@ -1341,6 +1353,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 =====================> diff --git a/include/solarsystem.h b/include/solarsystem.h index b7472398..f73a8cea 100644 --- a/include/solarsystem.h +++ b/include/solarsystem.h @@ -341,8 +341,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 diff --git a/src/novas.c b/src/novas.c index bf96c119..db45eebc 100644 --- a/src/novas.c +++ b/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,11 +5776,13 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig 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); + 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; @@ -5787,71 +5790,95 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig 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: *
    + *
  1. https://ssd.jpl.nasa.gov/planets/approx_pos.html
    2. *
  2. https://en.wikipedia.org/wiki/Orbital_elements
    3. - *
  3. https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf
    4. *
  4. https://orbitalofficial.com/
    5. + *
  5. https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf
    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 + * @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 @@ -5859,29 +5886,40 @@ int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orb, double 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() * */ diff --git a/src/super.c b/src/super.c index 153482ea..472e9cc3 100644 --- a/src/super.c +++ b/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 diff --git a/test/src/test-super.c b/test/src/test-super.c index 50615fe4..1378846a 100644 --- a/test/src/test-super.c +++ b/test/src/test-super.c @@ -2213,6 +2213,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; @@ -2277,6 +2318,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(); From 1e503f523e3f64b00471ea2a56aeb007090378b0 Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Thu, 14 Nov 2024 11:37:14 +0100 Subject: [PATCH 3/7] More complex orbital systems, such as planet equators, local Laplace planes... --- include/novas.h | 62 +++++++++++++++++++++++++++++--- src/novas.c | 84 ++++++++++++++++++++++++++++++++++++++----- test/src/test-super.c | 6 +--- 3 files changed, 133 insertions(+), 19 deletions(-) diff --git a/include/novas.h b/include/novas.h index 6a4a609d..4c1eb611 100644 --- a/include/novas.h +++ b/include/novas.h @@ -639,7 +639,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 of orbital reference plane (e.g. 90° - Decpole) + double Omega; ///< [rad] argument of ascending node of the orbital reference plane (e.g. RApole) +} 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: + *
    + *
  1. Up-to-date orbital elements for asteroids, comets, etc from the Minor Planet Center (MPC): + * https://minorplanetcenter.net/data
    2. + *
  2. Mean elements for planetary satellites from JPL Horizons: https://ssd.jpl.nasa.gov/sats/elem/
    3. + *
  3. Low accuracy mean elements for planets from JPL Horizons: + * https://ssd.jpl.nasa.gov/planets/approx_pos.html
    4. + *
* * @author Attila Kovacs * @since 1.2 @@ -648,9 +702,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; ///< 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 @@ -667,7 +719,7 @@ typedef struct { * @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. diff --git a/src/novas.c b/src/novas.c index db45eebc..f017470a 100644 --- a/src/novas.c +++ b/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, ¢er), 0); + prop_error(fn, make_planet(body->orbit.system.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, 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->Omega, 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; diff --git a/test/src/test-super.c b/test/src/test-super.c index 1378846a..50b66106 100644 --- a/test/src/test-super.c +++ b/test/src/test-super.c @@ -2215,7 +2215,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 = {}; @@ -2227,10 +2227,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; From 76b6dbee557504f41fee1340bd17ae861b70a342 Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Thu, 14 Nov 2024 13:37:42 +0100 Subject: [PATCH 4/7] doc updates --- CHANGELOG.md | 16 +++++++++++++--- README.md | 26 ++++++++++++++++++++++++++ include/novas.h | 23 ++++++++++++----------- 3 files changed, 51 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3b0c835c..0f2c0982 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -56,15 +56,25 @@ Changes expected for the next feature release, expected around 1 February 2025. You can also set other standard directory variables, as prescribed in the [GNU standard](https://www.gnu.org/prep/standards/html_node/Directory-Variables.html) to further customize the install locations. - - - SuperNOVAS headers now include each other as system-headers, not local headers. This is unlikely to affect anything - really but it is more proper for an installation of the library, and works with our own `Makefile` too. + + - #95: Added support for using orbital elements. `object.type` can now be set to `NOVAS_ORBITAL_OBJECT`, whose orbit + can be defined by the set of `novas_orbital_elements`, within a `novas_orbital_system`. You can initialize an + `object` with a set of orbital elements using `make_orbital_object()`. For such objects, `ephemeris()` will now + call on the new `novas_orbit_posvel()` to calculate positions. While orbital elements do not always yield precise + positions, they can for shorter periods, provided that the orbital elements are up-to-date. For example, the Minor + Planer Center (MPC) publishes orbital elements for all known asteroids and comets regularly. For newly discovered + objects, this may be the only and/or most accurate information available anywhere. - #97: Added `NOVAS_EMB` (Earth-Moon Barycenter) and `NOVAS_PLUTO_BARYCENTER` to `enum novas_planets` to distinguish from the planet center in calculations. + - SuperNOVAS headers now include each other as system-headers, not local headers. This is unlikely to affect anything + really but it is more proper for an installation of the library, and works with our own `Makefile` too. + ### Changed + - #96: Changed `object` structure to include `novas_orbital_elements` for `NOVAS_ORBITAL_OBJECT` types. + - #97: Updated `NOVAS_PLANETS`, `NOVAS_PLANET_NAMES_INIT` and `NOVAS_RMASS_INIT` macros to include the added planet constants. diff --git a/README.md b/README.md index e728c219..7d6dbeca 100644 --- a/README.md +++ b/README.md @@ -596,6 +596,24 @@ more generic ephemeris handling via a user-provided `novas_ephem_provider`. E.g. make_ephem_object("Ceres", 2000001, &ceres); ``` +As of version 1.2 you can also define solar system sources with orbital elements (such as the most up-to-date ones +provided by the [Minor Planet Center](https://minorplanetcenter.net/data) for asteroids, comets, etc.): + +```c + object NEA; // e.g. a Near-Earth Asteroid + + // Fill in the orbital parameters (pay attention to units!) + novas_orbital_elements orbit = NOVAS_ORBIT_INIT; + orbit.a = ...; + ... + + // Create an object for that orbit + make_orbital_object("NEAxxx", -1, &orbit, object); +``` + +Note, that even with orbital elements, you will, in general, require a planet calculator, to provide precise +positions for the Sun or planet, around which the orbit is defined. + Other than that, it's the same spiel as before, e.g.: ```c @@ -894,6 +912,14 @@ before that level of accuracy is reached. - Added `novas_planet_for_name()` function to return the NOVAS planet ID for a given (case insensitive) name. + - Added support for using orbital elements. `object.type` can now be set to `NOVAS_ORBITAL_OBJECT`, whose orbit + can be defined by the set of `novas_orbital_elements`, within a `novas_orbital_system`. You can initialize an + `object` with a set of orbital elements using `make_orbital_object()`. For such objects, `ephemeris()` will now + call on the new `novas_orbit_posvel()` to calculate positions. While orbital elements do not always yield precise + positions, they can for shorter periods, provided that the orbital elements are up-to-date. For example, the Minor + Planer Center (MPC) publishes orbital elements for all known asteroids and comets regularly. For newly discovered + objects, this may be the only and/or most accurate information available anywhere. + - Added `NOVAS_EMB` (Earth-Moon Barycenter) and `NOVAS_PLUTO_BARYCENTER` to `enum novas_planets` to distinguish from the corresponding planet centers in calculations. diff --git a/include/novas.h b/include/novas.h index 4c1eb611..101c6472 100644 --- a/include/novas.h +++ b/include/novas.h @@ -640,8 +640,8 @@ typedef struct { /** * 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. + * all major planets and barycenters (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` @@ -649,10 +649,11 @@ typedef struct { * * 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. + * North poles of planets from JPL Horizons, and use them as a proxy for the Laplace planes for 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). + * `NOVAS_EQUATORIAL_PLANE` and the `type` to `NOVAS_GCRS_EQUATOR` (since the plane is defined by the North pole + * orientation in GCRS equatorial RA/Dec). The obliquity is then 90° - Decpole (in radians), and `phi` + * is RApole (in radians). * * @author Attila Kovacs * @since 1.2 @@ -660,11 +661,11 @@ typedef struct { * @sa novas_orbital_elements */ typedef struct { - enum novas_planet center; ///< major body at the center of the orbit (default is NOVAS_SUN). + enum novas_planet center; ///< major planet or barycenter 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 of orbital reference plane (e.g. 90° - Decpole) - double Omega; ///< [rad] argument of ascending node of the orbital reference plane (e.g. RApole) + double obl; ///< [rad] relative obliquity of orbital reference plane (e.g. 90° - Decpole) + double Omega; ///< [rad] relative argument of ascending node of the orbital reference plane (e.g. RApole) } novas_orbital_system; @@ -702,8 +703,8 @@ typedef struct { * @sa enum NOVAS_ORBITAL_OBJECT */ typedef struct { - novas_orbital_system system; ///< orbital reference system assumed for the paramtetrization - double jd_tdb; ///< [day] Barycentri Dynamical Time (TDB) based Julian date of parameters. + novas_orbital_system system; ///< orbital reference system assumed for the parametrization + double jd_tdb; ///< [day] Barycentri Dynamical Time (TDB) based Julian date of the parameters. double a; ///< [AU] semi-major axis double e; ///< eccentricity double omega; ///< [rad] argument of periapsis / perihelion @@ -714,7 +715,7 @@ typedef struct { } novas_orbital_elements; /** - * Initializer for novas_orbital_elements for heliocentric orbits + * Initializer for novas_orbital_elements for heliocentric orbits using GCRS ecliptic pqrametrization. * * @author Attila Kovacs * @since 1.2 From 97a1b9390f830d6e620824a7393318a65de4723a Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Thu, 14 Nov 2024 23:42:18 +0100 Subject: [PATCH 5/7] Fixes, tweaks, and satellte orbit test --- include/novas.h | 32 +++++++++------- src/novas.c | 29 +++++--------- src/super.c | 37 ++++++++++++++++++ test/src/test-super.c | 89 ++++++++++++++++++++++++++++++++++++------- 4 files changed, 139 insertions(+), 48 deletions(-) diff --git a/include/novas.h b/include/novas.h index 101c6472..315f83f2 100644 --- a/include/novas.h +++ b/include/novas.h @@ -48,7 +48,6 @@ // for compiling your code, which may have relied on these, you can add '-DCOMPAT=1' to the // compiler options // -// #if COMPAT # include # include @@ -664,8 +663,10 @@ typedef struct { enum novas_planet center; ///< major planet or barycenter 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] relative obliquity of orbital reference plane (e.g. 90° - Decpole) - double Omega; ///< [rad] relative argument of ascending node of the orbital reference plane (e.g. RApole) + double obl; ///< [rad] relative obliquity of orbital reference plane + ///< (e.g. 90° - δpole) + double Omega; ///< [rad] relative argument of ascending node of the orbital reference plane + ///< (e.g. αpole + 90°) } novas_orbital_system; @@ -674,14 +675,14 @@ typedef struct { #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 + * Keplerian orbital elements for `NOVAS_ORBITAL_OBJECT` type. Orbital elements can be used to provide + * approximate positions for various Solar-system bodies. JPL publishes orbital elements (and their evolution) + * 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 describe 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 @@ -707,11 +708,12 @@ typedef struct { double jd_tdb; ///< [day] Barycentri Dynamical Time (TDB) based Julian date of the 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 omega; ///< [rad] argument of periapsis / perihelion, at the reference time + double Omega; ///< [rad] argument of ascending node, at the reference time double i; ///< [rad] inclination double M0; ///< [rad] mean anomaly at tjd - double n; ///< [rad/day] mean daily motion, i.e. (_GM_/_a_3)1/2 for the central body + double n; ///< [rad/day] mean daily motion, i.e. (_GM_/_a_3)1/2 for the central body, + ///< or 2π/T, where T is orbital period in days. } novas_orbital_elements; /** @@ -1406,6 +1408,8 @@ double novas_z_inv(double z); enum novas_planet novas_planet_for_name(const char *name); +int novas_set_orbital_pole(double ra, double dec, novas_orbital_system *sys); + 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); diff --git a/src/novas.c b/src/novas.c index f017470a..4194bba9 100644 --- a/src/novas.c +++ b/src/novas.c @@ -5795,41 +5795,31 @@ 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. + * @return 0 * */ 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 ca = cos(phi); + double sa = sin(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; + out[1] = sa * x + ca * cb * y - ca * sb * z; + out[2] = sb * y + cb * z; return 0; } @@ -5859,7 +5849,6 @@ static int equ2gcrs(double jd_tdb, const double *in, enum novas_equator_type sys } } - /** * Convert coordinates in an orbital system to GCRS equatorial coordinates * @@ -5925,13 +5914,13 @@ int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orbit, enum if(pos == vel) return novas_error(-1, EINVAL, fn, "output pos = vel (@ %p)", pos); - E = M = orbit->M0 + orbit->n * (jd_tdb - orbit->jd_tdb); + E = M = remainder(orbit->M0 + orbit->n * (jd_tdb - orbit->jd_tdb), TWOPI); // Iteratively determine E, using Newton-Raphson method... while(--i >= 0) { - double es = orbit->e * sin(E); - double ec = orbit->e * cos(E); - double dE = (E - es - M) / (1.0 - ec); + double esE = orbit->e * sin(E); + double ecE = orbit->e * cos(E); + double dE = (E - esE - M) / (1.0 - ecE); E -= dE; if(fabs(dE) < EPREC) diff --git a/src/super.c b/src/super.c index 472e9cc3..0f849640 100644 --- a/src/super.c +++ b/src/super.c @@ -1363,3 +1363,40 @@ enum novas_planet novas_planet_for_name(const char *name) { return novas_error(-1, EINVAL, fn, "No match for name: '%s'", name); } +/** + * Sets the orientation of an orbital system using the RA and DEC coordinates of the pole + * of the Laplace (or else equatorial) plane relative to which the orbital elements are + * defined. Orbital parameters of planetary satellites normally include the R.A. and + * declination of the pole of the local Laplace plane in which the Keplerian orbital elements + * are referenced. + * + * The system will become referenced to the equatorial plane, the relative obliquity is set + * to (90° - `dec`), while the argument of the ascending node ('Omega') is set to + * (90° + `ra`). + * + * NOTES: + *
    + *
  1. You should not expect much precision from the long-range orbital approximations for + * planetary satellites. For applications that require precision at any level, you should rely + * on appropriate ephemerides, or else on up-to-date short-term orbital elements.
    2. + *
+ * + * @param ra [h] the R.A. of the pole of the oribtal reference plane. + * @param dec [deg] the declination of the pole of the oribtal reference plane. + * @param[out] sys The orbital system + * @return 0 if successful, or else -1 (errno will be set to EINVAL) if the output `sys` + * pointer is NULL. + * + * @author Attila Kovacs + * @since 1.2 + */ +int novas_set_orbital_pole(double ra, double dec, novas_orbital_system *sys) { + if (!sys) + return novas_error(-1, EINVAL, "novas_set_orbital_pole", "input system is NULL"); + + sys->plane = NOVAS_EQUATORIAL_PLANE; + sys->obl = remainder((90.0 - dec) * DEGREE, TWOPI); + sys->Omega = remainder((ra + 6.0) * HOURANGLE, TWOPI); + + return 0; +} diff --git a/test/src/test-super.c b/test/src/test-super.c index 50b66106..e7198219 100644 --- a/test/src/test-super.c +++ b/test/src/test-super.c @@ -2213,20 +2213,21 @@ static int test_planet_for_name() { } -static int test_orbit_posvel() { +static int test_orbit_place() { object ceres = {}; novas_orbital_elements orbit = NOVAS_ORBIT_INIT; 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 + // Nov 14 0 UTC, geocentric from JPL Horizons. + double tjd = 2460628.50079861; // 0 UT as TT. + double RA0 = 19.684415; + double DEC0 = -28.62084; + double rv0 = 21.4255198; // km/s + double r = 3.32557776285144; // AU int n = 0; + // Orbital Parameters for JD 2460600.5 from MPC orbit.jd_tdb = 2460600.5; orbit.a = 2.7666197; orbit.e = 0.079184; @@ -2239,13 +2240,71 @@ static int test_orbit_posvel() { 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; + if(!is_ok("orbit_place", 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++; + if(!is_equal("orbit_place:ra", pos.ra, RA0, 1e-5)) n++; + if(!is_equal("orbit_place:dec", pos.dec, DEC0, 1e-4)) n++; + if(!is_equal("orbit_place:dist", pos.dis, r, 1e-4)) n++; + if(!is_equal("orbit_place:vrad", pos.rv, rv0, 1e-2)) n++; + + return n; +} + + +static int test_orbit_posvel_callisto() { + novas_orbital_elements orbit = NOVAS_ORBIT_INIT; + novas_orbital_system *sys = &orbit.system; + double pos0[3] = {}, pos[3] = {}, vel[3] = {}, ra, dec, dra, ddec; + int i; + + // 2000-01-01 12 UT, geocentric from JPL Horizons. + + double dist = 4.62117513332102; + double lt = 0.00577551831217194 * dist; // day + double tjd = 2451545.00079861 - lt; // 0 UT as TT, corrected or light time + double dyr; + + double RA0 = 23.86983 * DEGREE; + double DEC0 = 8.59590 * DEGREE; + + double dRA = (23.98606 * DEGREE - RA0) / cos(DEC0); + double dDEC = (8.64868 * DEGREE - DEC0); + int n = 0; + + // Planet pos; + radec2vector(RA0 / HOURANGLE, DEC0 / DEGREE, dist, pos0); + + // Callisto's parameters from JPL Horizons + // https://ssd.jpl.nasa.gov/sats/elem/sep.html + // 1882700. 0.007 43.8 87.4 0.3 309.1 16.690440 277.921 577.264 268.7 64.8 + sys->center = NOVAS_JUPITER; + novas_set_orbital_pole(268.7 / 15.0, 64.8, sys); + + orbit.jd_tdb = NOVAS_JD_J2000; + dyr = (tjd - orbit.jd_tdb) / 365.25; //(yr) time since J2000.0 + + orbit.a = 1882700.0 * 1e3 / AU; + orbit.e = 0.007; + orbit.omega = remainder(43.8 * DEGREE + TWOPI * dyr / 277.921, TWOPI); + orbit.M0 = 87.4 * DEGREE; + orbit.i = 0.3 * DEGREE; + orbit.Omega = remainder(309.1 * DEGREE + TWOPI * dyr / 577.264, TWOPI); + orbit.n = TWOPI / 16.690440; + + if(!is_ok("orbit_posvel_callisto", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, pos, vel))) return 1; + + for(i = 3 ;--i >= 0; ) pos[i] += pos0[i]; + vector2radec(pos, &ra, &dec); + + ra *= HOURANGLE; + dec *= DEGREE; + + dra = (ra - RA0) * cos(DEC0); + ddec = (dec - DEC0); + + if(!is_equal("orbit_posvel_callisto:dist", hypot(dra, ddec) / ARCSEC, hypot(dRA, dDEC) / ARCSEC, 15.0)) n++; + if(!is_equal("orbit_posvel_callisto:ra", dra / ARCSEC, dRA / ARCSEC, 15.0)) n++; + if(!is_equal("orbit_posvel_callisto:dec", ddec / ARCSEC, dDEC / ARCSEC, 15.0)) n++; return n; } @@ -2314,7 +2373,9 @@ int main(int argc, char *argv[]) { if(test_naif_to_novas_planet()) n++; if(test_planet_for_name()) n++; - if(test_orbit_posvel()) n++; + + if(test_orbit_place()) n++; + if(test_orbit_posvel_callisto()) n++; n += test_dates(); From fe55c5819bbf7db4e7fa80b571f81150ac722a7f Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Fri, 15 Nov 2024 08:49:00 +0100 Subject: [PATCH 6/7] Change orbital angle units to deg (NOVAS convention), and add precession terms --- include/novas.h | 16 +++++++++------- src/novas.c | 34 +++++++++++++++++++++++----------- src/super.c | 4 ++-- test/src/test-super.c | 25 ++++++++++++------------- 4 files changed, 46 insertions(+), 33 deletions(-) diff --git a/include/novas.h b/include/novas.h index 315f83f2..1f38001f 100644 --- a/include/novas.h +++ b/include/novas.h @@ -708,12 +708,14 @@ typedef struct { double jd_tdb; ///< [day] Barycentri Dynamical Time (TDB) based Julian date of the parameters. double a; ///< [AU] semi-major axis double e; ///< eccentricity - double omega; ///< [rad] argument of periapsis / perihelion, at the reference time - double Omega; ///< [rad] argument of ascending node, at the reference time - double i; ///< [rad] inclination - double M0; ///< [rad] mean anomaly at tjd - double n; ///< [rad/day] mean daily motion, i.e. (_GM_/_a_3)1/2 for the central body, - ///< or 2π/T, where T is orbital period in days. + double omega; ///< [deg] argument of periapsis / perihelion, at the reference time + double Omega; ///< [deg] argument of ascending node on the reference plane, at the reference time + double i; ///< [deg] inclination of orbit to reference plane + double M0; ///< [deg] mean anomaly at tjd + double n; ///< [deg/day] mean daily motion, i.e. (_GM_/_a_3)1/2 for the central body, + ///< or 360/T, where T is orbital period in days. + double apsis_period; ///< [day] Precession period of the apsis, if known. + double node_period; ///< [day] Precession period of the ascending node, if known. } novas_orbital_elements; /** @@ -722,7 +724,7 @@ typedef struct { * @author Attila Kovacs * @since 1.2 */ -#define NOVAS_ORBIT_INIT { NOVAS_ORBITAL_SYSTEM_INIT, 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, 0.0, 0.0 } /** * Celestial object of interest. diff --git a/src/novas.c b/src/novas.c index 4194bba9..32be3101 100644 --- a/src/novas.c +++ b/src/novas.c @@ -5799,8 +5799,8 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig * 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 theta [deg] polar angle of original pole in the new system + * @param phi [deg] 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 * @@ -5812,6 +5812,9 @@ static int change_pole(const double *in, double theta, double phi, double *out) y = in[1]; z = in[2]; + theta *= DEGREE; + phi *= DEGREE; + double ca = cos(phi); double sa = sin(phi); double cb = cos(theta); @@ -5903,7 +5906,7 @@ static int orbit2gcrs(double jd_tdb, const novas_orbital_system *sys, enum novas 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 M, E, nu, r; + double dt, M, E, nu, r, omega, Omega; double cO, sO, ci, si, co, so; double xx, yx, zx, xy, yy, zy; int i = novas_inv_max_iter; @@ -5914,7 +5917,8 @@ int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orbit, enum if(pos == vel) return novas_error(-1, EINVAL, fn, "output pos = vel (@ %p)", pos); - E = M = remainder(orbit->M0 + orbit->n * (jd_tdb - orbit->jd_tdb), TWOPI); + dt = (jd_tdb - orbit->jd_tdb); + E = M = remainder(orbit->M0 + orbit->n * dt, 360.0) * DEGREE; // Iteratively determine E, using Newton-Raphson method... while(--i >= 0) { @@ -5933,13 +5937,21 @@ int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orbit, enum 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)); + omega = orbit->omega * DEGREE; + if(orbit->apsis_period > 0.0) + omega += TWOPI * remainder(dt / orbit->apsis_period, 1.0); + + Omega = orbit->Omega * DEGREE; + if(orbit->node_period > 0.0) + Omega += TWOPI * remainder(dt / orbit->node_period, 1.0); + // pos = Rz(-Omega) . Rx(-i) . Rz(-omega) . orb - cO = cos(orbit->Omega); - sO = sin(orbit->Omega); - ci = cos(orbit->i); - si = sin(orbit->i); - co = cos(orbit->omega); - so = sin(orbit->omega); + cO = cos(Omega); + sO = sin(Omega); + ci = cos(orbit->i * DEGREE); + si = sin(orbit->i * DEGREE); + co = cos(omega); + so = sin(omega); // Rotation matrix // See https://en.wikipedia.org/wiki/Euler_angles @@ -5965,7 +5977,7 @@ int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orbit, enum } if(vel) { - double v = orbit->n * orbit->a * orbit->a / r; // [AU/day] + double v = orbit->n * DEGREE * orbit->a * orbit->a / r; // [AU/day] double x = -v * sin(E); double y = v * sqrt(1.0 - orbit->e * orbit->e) * cos(E); diff --git a/src/super.c b/src/super.c index 0f849640..0bdd1258 100644 --- a/src/super.c +++ b/src/super.c @@ -1395,8 +1395,8 @@ int novas_set_orbital_pole(double ra, double dec, novas_orbital_system *sys) { return novas_error(-1, EINVAL, "novas_set_orbital_pole", "input system is NULL"); sys->plane = NOVAS_EQUATORIAL_PLANE; - sys->obl = remainder((90.0 - dec) * DEGREE, TWOPI); - sys->Omega = remainder((ra + 6.0) * HOURANGLE, TWOPI); + sys->obl = remainder(90.0 - dec, 360.0); + sys->Omega = remainder(15.0 * ra + 90.0, 360.0); return 0; } diff --git a/test/src/test-super.c b/test/src/test-super.c index e7198219..8e871ca4 100644 --- a/test/src/test-super.c +++ b/test/src/test-super.c @@ -2231,18 +2231,18 @@ static int test_orbit_place() { 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; + orbit.i = 10.5879; + orbit.omega = 73.28579; + orbit.Omega = 80.25414; + orbit.M0 = 145.84905; + orbit.n = 0.21418047; make_observer_at_geocenter(&obs); make_orbital_object("Ceres", -1, &orbit, &ceres); if(!is_ok("orbit_place", place(tjd, &ceres, &obs, ut12tt, NOVAS_TOD, NOVAS_REDUCED_ACCURACY, &pos))) return 1; - if(!is_equal("orbit_place:ra", pos.ra, RA0, 1e-5)) n++; + if(!is_equal("orbit_place:ra", pos.ra, RA0, 1e-5 / cos(DEC0 * DEGREE))) n++; if(!is_equal("orbit_place:dec", pos.dec, DEC0, 1e-4)) n++; if(!is_equal("orbit_place:dist", pos.dis, r, 1e-4)) n++; if(!is_equal("orbit_place:vrad", pos.rv, rv0, 1e-2)) n++; @@ -2262,7 +2262,6 @@ static int test_orbit_posvel_callisto() { double dist = 4.62117513332102; double lt = 0.00577551831217194 * dist; // day double tjd = 2451545.00079861 - lt; // 0 UT as TT, corrected or light time - double dyr; double RA0 = 23.86983 * DEGREE; double DEC0 = 8.59590 * DEGREE; @@ -2281,15 +2280,15 @@ static int test_orbit_posvel_callisto() { novas_set_orbital_pole(268.7 / 15.0, 64.8, sys); orbit.jd_tdb = NOVAS_JD_J2000; - dyr = (tjd - orbit.jd_tdb) / 365.25; //(yr) time since J2000.0 - orbit.a = 1882700.0 * 1e3 / AU; orbit.e = 0.007; - orbit.omega = remainder(43.8 * DEGREE + TWOPI * dyr / 277.921, TWOPI); - orbit.M0 = 87.4 * DEGREE; - orbit.i = 0.3 * DEGREE; - orbit.Omega = remainder(309.1 * DEGREE + TWOPI * dyr / 577.264, TWOPI); + orbit.omega = 43.8; + orbit.M0 = 87.4; + orbit.i = 0.3; + orbit.Omega = 309.1; orbit.n = TWOPI / 16.690440; + orbit.apsis_period = 277.921 * 365.25; + orbit.node_period = 577.264 * 365.25; if(!is_ok("orbit_posvel_callisto", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, pos, vel))) return 1; From 15911b9d80abda82665cc6d7525ecb16fdb125cb Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Fri, 15 Nov 2024 13:11:06 +0100 Subject: [PATCH 7/7] More tests and fixes --- CHANGELOG.md | 13 ++++----- README.md | 13 ++++----- include/novas.h | 19 ++++++++----- src/novas.c | 34 +++++++++++++++-------- test/src/test-errors.c | 61 ++++++++++++++++++++++++++++++++++++++++++ test/src/test-super.c | 30 ++++++++++++++++++--- 6 files changed, 137 insertions(+), 33 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0f2c0982..87fc9857 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -58,12 +58,13 @@ Changes expected for the next feature release, expected around 1 February 2025. install locations. - #95: Added support for using orbital elements. `object.type` can now be set to `NOVAS_ORBITAL_OBJECT`, whose orbit - can be defined by the set of `novas_orbital_elements`, within a `novas_orbital_system`. You can initialize an - `object` with a set of orbital elements using `make_orbital_object()`. For such objects, `ephemeris()` will now - call on the new `novas_orbit_posvel()` to calculate positions. While orbital elements do not always yield precise - positions, they can for shorter periods, provided that the orbital elements are up-to-date. For example, the Minor - Planer Center (MPC) publishes orbital elements for all known asteroids and comets regularly. For newly discovered - objects, this may be the only and/or most accurate information available anywhere. + can be defined by the set of `novas_orbital_elements`, relative to a `novas_orbital_system`. You can initialize an + `object` with a set of orbital elements using `make_orbital_object()`, and for planetary satellite orbits you might + use `novas_set_orbital_pole()`. For orbital objects, `ephemeris()` will call on the new `novas_orbit_posvel()` to + calculate positions. While orbital elements do not always yield precise positions, they can for shorter periods, + provided that the orbital elements are up-to-date. For example, the Minor Planer Center (MPC) publishes accurate + orbital elements for all known asteroids and comets regularly. For newly discovered objects, this may be the only + and/or most accurate information available anywhere. - #97: Added `NOVAS_EMB` (Earth-Moon Barycenter) and `NOVAS_PLUTO_BARYCENTER` to `enum novas_planets` to distinguish from the planet center in calculations. diff --git a/README.md b/README.md index 7d6dbeca..296f4578 100644 --- a/README.md +++ b/README.md @@ -913,12 +913,13 @@ before that level of accuracy is reached. - Added `novas_planet_for_name()` function to return the NOVAS planet ID for a given (case insensitive) name. - Added support for using orbital elements. `object.type` can now be set to `NOVAS_ORBITAL_OBJECT`, whose orbit - can be defined by the set of `novas_orbital_elements`, within a `novas_orbital_system`. You can initialize an - `object` with a set of orbital elements using `make_orbital_object()`. For such objects, `ephemeris()` will now - call on the new `novas_orbit_posvel()` to calculate positions. While orbital elements do not always yield precise - positions, they can for shorter periods, provided that the orbital elements are up-to-date. For example, the Minor - Planer Center (MPC) publishes orbital elements for all known asteroids and comets regularly. For newly discovered - objects, this may be the only and/or most accurate information available anywhere. + can be defined by the set of `novas_orbital_elements`, relative to a `novas_orbital_system`. You can initialize an + `object` with a set of orbital elements using `make_orbital_object()`, and for planetary satellite orbits you might + use `novas_set_orbital_pole()`. For orbital objects, `ephemeris()` will call on the new `novas_orbit_posvel()` to + calculate positions. While orbital elements do not always yield precise positions, they can for shorter periods, + provided that the orbital elements are up-to-date. For example, the Minor Planer Center (MPC) publishes accurate + orbital elements for all known asteroids and comets regularly. For newly discovered objects, this may be the only + and/or most accurate information available anywhere. - Added `NOVAS_EMB` (Earth-Moon Barycenter) and `NOVAS_PLUTO_BARYCENTER` to `enum novas_planets` to distinguish from the corresponding planet centers in calculations. diff --git a/include/novas.h b/include/novas.h index 1f38001f..d9571a40 100644 --- a/include/novas.h +++ b/include/novas.h @@ -376,6 +376,10 @@ enum novas_equator_type { NOVAS_GCRS_EQUATOR ///< Geocentric Celestial Reference System (GCRS). }; +/// The number of equator types we define. +/// @since 1.2 +#define NOVAS_EQUATOR_TYPES (NOVAS_GCRS_EQUATOR + 1) + /** * Constants that determine the type of dynamical system type for gcrs2equ() */ @@ -660,13 +664,14 @@ typedef struct { * @sa novas_orbital_elements */ typedef struct { - enum novas_planet center; ///< major planet or barycenter 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]). + enum novas_planet center; ///< major planet or barycenter at the center of the orbit. + enum novas_reference_plane plane; ///< reference plane NOVAS_ECLIPTIC_PLANE or NOVAS_EQUATORIAL_PLANE + enum novas_equator_type type; ///< the type of equator in which orientation is referenced (NOVAS_TRUE_EQUATOR, + ///< NOVAS_MEAN_EQUATOR, or NOVAS_GCRS_EQUATOR). double obl; ///< [rad] relative obliquity of orbital reference plane - ///< (e.g. 90° - δpole) + ///< (e.g. 90° - δpole) double Omega; ///< [rad] relative argument of ascending node of the orbital reference plane - ///< (e.g. αpole + 90°) + ///< (e.g. αpole + 90°) } novas_orbital_system; @@ -710,8 +715,8 @@ typedef struct { double e; ///< eccentricity double omega; ///< [deg] argument of periapsis / perihelion, at the reference time double Omega; ///< [deg] argument of ascending node on the reference plane, at the reference time - double i; ///< [deg] inclination of orbit to reference plane - double M0; ///< [deg] mean anomaly at tjd + double i; ///< [deg] inclination of orbit to the reference plane + double M0; ///< [deg] mean anomaly at the reference time double n; ///< [deg/day] mean daily motion, i.e. (_GM_/_a_3)1/2 for the central body, ///< or 360/T, where T is orbital period in days. double apsis_period; ///< [day] Precession period of the apsis, if known. diff --git a/src/novas.c b/src/novas.c index 32be3101..d57007d4 100644 --- a/src/novas.c +++ b/src/novas.c @@ -864,8 +864,14 @@ 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 && ss_body->type != NOVAS_ORBITAL_OBJECT) - return novas_error(-1, EINVAL, fn, "object is not solar-system type: type=%d", ss_body->type); + switch(ss_body->type) { + case NOVAS_PLANET: + case NOVAS_EPHEM_OBJECT: + case NOVAS_ORBITAL_OBJECT: + break; + default: + 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); @@ -5779,10 +5785,8 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig 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]; + pos[i] += pos0[i]; + vel[i] += vel0[i]; } break; @@ -5863,13 +5867,20 @@ static int equ2gcrs(double jd_tdb, const double *in, enum novas_equator_type sys * */ static int orbit2gcrs(double jd_tdb, const novas_orbital_system *sys, enum novas_accuracy accuracy, double *vec) { + static const char *fn = "orbit2gcrs"; + if(sys->obl) change_pole(vec, sys->obl, sys->Omega, 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); + if(sys->plane == NOVAS_ECLIPTIC_PLANE) { + if(ecl2equ_vec(jd_tdb, sys->type, accuracy, vec, vec) != 0) + return novas_trace(fn, -1, 0); + } + else if(sys->plane != NOVAS_EQUATORIAL_PLANE) + return novas_error(-1, EINVAL, fn, "invalid orbital system reference plane type: %d", sys->type); + + prop_error(fn, equ2gcrs(jd_tdb, vec, sys->type, vec), 0); return 0; } @@ -5892,8 +5903,9 @@ static int orbit2gcrs(double jd_tdb, const novas_orbital_system *sys, enum novas * @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). + * or if the position and velocity output vectors are the same or the orbital + * system is ill defined (errno set to EINVAL), or if the calculation did not converge (errno set to + * ECANCELED), or * * @sa ephemeris() * @sa novas_geom_posvel() diff --git a/test/src/test-errors.c b/test/src/test-errors.c index 86a4378d..81bf5a65 100644 --- a/test/src/test-errors.c +++ b/test/src/test-errors.c @@ -1491,6 +1491,64 @@ static int test_planet_for_name() { return n; } +static int test_make_orbital_object() { + int n = 0; + novas_orbital_elements orbit = {}; + object body = {}; + + if(check("make_orbital_object:orbit", -1, make_orbital_object("blah", -1, NULL, &body))) n++; + if(check("make_orbital_object:body", -1, make_orbital_object("blah", -1, &orbit, NULL))) n++; + if(check("make_orbital_object:orbit+body", -1, make_orbital_object("blah", -1, NULL, NULL))) n++; + + return n; +} + +static int test_orbit_posvel() { + extern int novas_inv_max_iter; + + int n = 0; + double pos[3] = {}, vel[3] = {}; + int saved = novas_inv_max_iter; + novas_orbital_elements orbit = NOVAS_ORBIT_INIT; + + orbit.a = 1.0; + + if(check("set_orbital_pole:orbit", -1, novas_orbit_posvel(0.0, NULL, NOVAS_REDUCED_ACCURACY, pos, vel))) n++; + if(check("set_orbital_pole:pos=vel", -1, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, pos, pos))) n++; + if(check("set_orbital_pole:pos=vel:NULL", -1, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, NULL, NULL))) n++; + + if(check("set_orbital_pole:orbit:converge", 0, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, pos, vel))) n++; + + novas_inv_max_iter = 0; + if(check("set_orbital_pole:orbit:converge", -1, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, pos, vel))) n++; + else if(check("set_orbital_pole:orbit:converge:errno", ECANCELED, errno)) n++; + novas_inv_max_iter = saved; + + orbit.system.type = -1; + if(check("set_orbital_pole:orbit:type:-1", -1, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, pos, vel))) n++; + + orbit.system.type = NOVAS_EQUATOR_TYPES; + if(check("set_orbital_pole:orbit:type:hi", -1, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, pos, vel))) n++; + + orbit.system.plane = NOVAS_EQUATORIAL_PLANE; + orbit.system.type = NOVAS_EQUATOR_TYPES; + if(check("set_orbital_pole:orbit:type:-1:eq", -1, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, pos, vel))) n++; + + orbit.system.type = NOVAS_GCRS_EQUATOR; + orbit.system.plane = -1; + if(check("set_orbital_pole:orbit:plane:-1", -1, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, pos, vel))) n++; + + return n; +} + +static int test_set_orbital_pole() { + int n = 0; + + if(check("set_orbital_pole:orbit", -1, novas_set_orbital_pole(0.0, 0.0, NULL))) n++; + + return n; +} + int main() { int n = 0; @@ -1620,6 +1678,9 @@ int main() { if(test_naif_to_novas_planet()) n++; if(test_planet_for_name()) n++; + if(test_make_orbital_object()) n++; + if(test_set_orbital_pole()) n++; + if(test_orbit_posvel()) n++; if(n) fprintf(stderr, " -- FAILED %d tests\n", n); else fprintf(stderr, " -- OK\n"); diff --git a/test/src/test-super.c b/test/src/test-super.c index 8e871ca4..d427fcf7 100644 --- a/test/src/test-super.c +++ b/test/src/test-super.c @@ -2251,10 +2251,12 @@ static int test_orbit_place() { } + + static int test_orbit_posvel_callisto() { novas_orbital_elements orbit = NOVAS_ORBIT_INIT; novas_orbital_system *sys = &orbit.system; - double pos0[3] = {}, pos[3] = {}, vel[3] = {}, ra, dec, dra, ddec; + double pos0[3] = {}, pos[3] = {}, vel[3] = {}, pos1[3] = {}, vel1[3] = {}, ra, dec, dra, ddec; int i; // 2000-01-01 12 UT, geocentric from JPL Horizons. @@ -2271,7 +2273,7 @@ static int test_orbit_posvel_callisto() { int n = 0; // Planet pos; - radec2vector(RA0 / HOURANGLE, DEC0 / DEGREE, dist, pos0); + radec2vector(RA0 / HOURANGLE, DEC0 / DEGREE, dist, pos1); // Callisto's parameters from JPL Horizons // https://ssd.jpl.nasa.gov/sats/elem/sep.html @@ -2291,8 +2293,9 @@ static int test_orbit_posvel_callisto() { orbit.node_period = 577.264 * 365.25; if(!is_ok("orbit_posvel_callisto", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, pos, vel))) return 1; + memcpy(pos0, pos, sizeof(pos)); - for(i = 3 ;--i >= 0; ) pos[i] += pos0[i]; + for(i = 3 ;--i >= 0; ) pos[i] += pos1[i]; vector2radec(pos, &ra, &dec); ra *= HOURANGLE; @@ -2305,6 +2308,27 @@ static int test_orbit_posvel_callisto() { if(!is_equal("orbit_posvel_callisto:ra", dra / ARCSEC, dRA / ARCSEC, 15.0)) n++; if(!is_equal("orbit_posvel_callisto:dec", ddec / ARCSEC, dDEC / ARCSEC, 15.0)) n++; + + + if(!is_ok("orbit_posvel_callisto:vel:null", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, pos1, NULL))) n++; + if(!is_ok("orbit_posvel_callisto:vel:null:check", check_equal_pos(pos1, pos0, 1e-8))) n++; + + if(!is_ok("orbit_posvel_callisto:pos:null", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, NULL, vel1))) n++; + if(!is_ok("orbit_posvel_callisto:pos:null:check", check_equal_pos(vel1, vel, 1e-8))) n++; + + + sys->type = NOVAS_MEAN_EQUATOR; + if(!is_ok("orbit_posvel_callisto:mod", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, pos1, NULL))) n++; + precession(tjd, pos0, NOVAS_JD_J2000, pos); + j2000_to_gcrs(pos, pos); + if(!is_ok("orbit_posvel_callisto:mod:check", check_equal_pos(pos1, pos, 1e-8))) n++; + + sys->type = NOVAS_TRUE_EQUATOR; + if(!is_ok("orbit_posvel_callisto:mod", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, pos1, NULL))) n++; + tod_to_j2000(tjd, NOVAS_FULL_ACCURACY, pos0, pos); + j2000_to_gcrs(pos, pos); + if(!is_ok("orbit_posvel_callisto:mod:check", check_equal_pos(pos1, pos, 1e-8))) n++; + return n; }