From 292eaf193350a0336d71e9ea0b28b5c13fa3e517 Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Wed, 13 Nov 2024 21:24:46 +0100 Subject: [PATCH] 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 566d1668..36212d11 100644 --- a/include/novas.h +++ b/include/novas.h @@ -577,6 +577,16 @@ enum novas_nutation_direction { NUTATE_MEAN_TO_TRUE }; +/** + * The plane in which values, such as orbital parameyters are referenced. + * + * @since 1.2 + */ +enum novas_reference_plane { + NOVAS_ECLIPTIC_PLANE = 0, ///< the plane of the ecliptic + NOVAS_EQUATORIAL_PLANE ///< The plane of the equator +}; + /** * Fundamental Delaunay arguments of the Sun and Moon, from Simon section 3.4(b.3), * precession = 5028.8200 arcsec/cy) @@ -632,7 +642,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 @@ -640,7 +652,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; /** @@ -649,7 +661,7 @@ typedef struct { * @author Attila Kovacs * @since 1.2 */ -#define NOVAS_ORBIT_INIT { NOVAS_SUN, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } +#define NOVAS_ORBIT_INIT { NOVAS_SUN, NOVAS_ECLIPTIC_PLANE, NOVAS_GCRS_EQUATOR, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } /** * Celestial object of interest. @@ -1335,6 +1347,10 @@ double novas_z_inv(double z); enum novas_planet novas_planet_for_name(const char *name); +int make_orbital_object(const char *name, long num, const novas_orbital_elements *orbit, object *body); + +int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orb, enum novas_accuracy accuracy, double *pos, double *vel); + // <================= END of SuperNOVAS API =====================> diff --git a/include/solarsystem.h b/include/solarsystem.h index 7153ca3a..bae376f6 100644 --- a/include/solarsystem.h +++ b/include/solarsystem.h @@ -340,8 +340,6 @@ long novas_to_naif_planet(enum novas_planet id); long novas_to_dexxx_planet(enum novas_planet id); -int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orb, double *pos, double *vel); - /// \cond PRIVATE diff --git a/src/novas.c b/src/novas.c index 948c9674..e419f10c 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 1d2deb98..d61b5a92 100644 --- a/test/src/test-super.c +++ b/test/src/test-super.c @@ -2207,6 +2207,47 @@ static int test_planet_for_name() { } +static int test_orbit_posvel() { + object ceres = {}; + novas_orbital_elements orbit; + observer obs = {}; + sky_pos pos = {}; + + // This is for SMA, so topocentric with Earth rot... + double tjd = NOVAS_JD_MJD0 + 60627.796269; + double RA0 = 19.679722; // 19:40:47.000 + double DEC0 = -28.633014; // -28:37:58.849 + double rv0 = 21.138; // km/s + double r = 3.323; // AU + int n = 0; + + orbit.center = NOVAS_SUN; + orbit.plane = NOVAS_ECLIPTIC_PLANE; + orbit.sys = NOVAS_GCRS_EQUATOR; + + orbit.jd_tdb = 2460600.5; + orbit.a = 2.7666197; + orbit.e = 0.079184; + orbit.i = 10.5879 * DEGREE; + orbit.omega = 73.28579 * DEGREE; + orbit.Omega = 80.25414 * DEGREE; + orbit.M0 = 145.84905 * DEGREE; + orbit.n = 0.21418047 * DEGREE; + + make_observer_at_geocenter(&obs); + make_orbital_object("Ceres", -1, &orbit, &ceres); + + if(!is_ok("orbit_posvel", place(tjd, &ceres, &obs, ut12tt, NOVAS_TOD, NOVAS_REDUCED_ACCURACY, &pos))) return 1; + + // TODO refine with HORIZONS data... + if(!is_equal("orbit_posvel:ra", pos.ra, RA0, 5e-5)) n++; + if(!is_equal("orbit_posvel:dec", pos.dec, DEC0, 1e-3)) n++; + if(!is_equal("orbit_posvel:dist", pos.dis, r, 2e-3)) n++; + if(!is_equal("orbit_posvel:vrad", pos.rv, rv0, 0.5)) n++; + + return n; +} + int main(int argc, char *argv[]) { int n = 0; @@ -2271,6 +2312,7 @@ int main(int argc, char *argv[]) { if(test_naif_to_novas_planet()) n++; if(test_planet_for_name()) n++; + if(test_orbit_posvel()) n++; n += test_dates();