From c89e6241faae54056b657a336a9f9f8a1dd7d5c7 Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Thu, 14 Nov 2024 11:37:14 +0100 Subject: [PATCH] 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 36212d11..6344ee2a 100644 --- a/include/novas.h +++ b/include/novas.h @@ -633,7 +633,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 to reference plane (e.g. 90° - Decpole) + double phi; ///< [rad] argument of ascending node (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 @@ -642,9 +696,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; ///< The 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 @@ -661,7 +713,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 e419f10c..58c6719d 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->phi, 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 d61b5a92..e4f709dd 100644 --- a/test/src/test-super.c +++ b/test/src/test-super.c @@ -2209,7 +2209,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 = {}; @@ -2221,10 +2221,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;