From e85a28bdbd842d8fd9d58b7788df920c4fb55234 Mon Sep 17 00:00:00 2001 From: Attila Kovacs Date: Sun, 17 Nov 2024 11:05:44 +0100 Subject: [PATCH] Tweak novas_orbital_system and rename to novas_set_obsys_pole(). --- CHANGELOG.md | 12 ++-- README.md | 6 +- include/novas.h | 16 +++-- src/novas.c | 135 +++++++++++++++++++++++++++++++---------- src/super.c | 6 +- test/src/test-errors.c | 81 ++++++++++++++++++++----- test/src/test-super.c | 105 +++++++++++++++++++++++++++++--- 7 files changed, 292 insertions(+), 69 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b791513c..1ed985ec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,9 +7,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [1.2.0-rc2] - 2024-11-18 +## [Unreleased] -Release candidate for the next feature release, expected around 1 February 2025. +Changes exzpected for the next feature release, expected around 1 February 2025. ### Added @@ -57,10 +57,10 @@ Release candidate for the next feature release, expected around 1 February 2025. [GNU standard](https://www.gnu.org/prep/standards/html_node/Directory-Variables.html) to further customize the install locations. - - #95: Added support for using orbital elements. `object.type` can now be set to `NOVAS_ORBITAL_OBJECT`, whose orbit + - #95, #98: 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`, 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 + use `novas_set_orbsys_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 @@ -69,6 +69,9 @@ Release candidate for the next feature release, expected around 1 February 2025. - #97: Added `NOVAS_EMB` (Earth-Moon Barycenter) and `NOVAS_PLUTO_BARYCENTER` to `enum novas_planets` to distinguish from the planet center in calculations. + - #98: Added `gcrs_to_tod()` / `tod_to_gcrs()` and `gcrs_to_mod()` / `mod_to_gcrs()` vector conversion functions for + convenience. + - 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. @@ -80,6 +83,7 @@ Release candidate for the next feature release, expected around 1 February 2025. constants. + ## [1.1.1] - 2024-10-28 Bug fix release. Nothing too scary, mainly just a collection of smaller fixes and improvements. diff --git a/README.md b/README.md index 296f4578..2ed7103c 100644 --- a/README.md +++ b/README.md @@ -915,7 +915,7 @@ before that level of accuracy is reached. - 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`, 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 + use `novas_set_orbsys_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 @@ -923,7 +923,9 @@ before that level of accuracy is reached. - Added `NOVAS_EMB` (Earth-Moon Barycenter) and `NOVAS_PLUTO_BARYCENTER` to `enum novas_planets` to distinguish from the corresponding planet centers in calculations. - + + - Added `gcrs_to_tod()` / `tod_to_gcrs()` and `gcrs_to_mod()` / `mod_to_gcrs()` vector conversion functions for + convenience. ### Refinements to the NOVAS C API diff --git a/include/novas.h b/include/novas.h index 5753b80a..8b238f1c 100644 --- a/include/novas.h +++ b/include/novas.h @@ -715,15 +715,14 @@ typedef struct { * @author Attila Kovacs * @since 1.2 * - * @sa novas_set_orbital_pole() + * @sa novas_set_obsys_pole() * @sa novas_orbital_elements * @sa NOVAS_ORBITAL_SYSTEM_INIT */ typedef struct { 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). + enum novas_reference_system type; ///< the coordinate reference system used for the reference plane and orbitals. 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 @@ -737,7 +736,7 @@ typedef struct { * @author Attila Kovacs * @since 1.2 */ -#define NOVAS_ORBITAL_SYSTEM_INIT { NOVAS_SUN, NOVAS_ECLIPTIC_PLANE, NOVAS_GCRS_EQUATOR, 0.0, 0.0 } +#define NOVAS_ORBITAL_SYSTEM_INIT { NOVAS_SUN, NOVAS_ECLIPTIC_PLANE, NOVAS_GCRS, 0.0, 0.0 } /** * Keplerian orbital elements for `NOVAS_ORBITAL_OBJECT` type. Orbital elements can be used to provide @@ -1486,12 +1485,19 @@ 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 novas_set_orbsys_pole(enum novas_reference_system type, 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); +int gcrs_to_tod(double jd_tdb, enum novas_accuracy accuracy, const double *in, double *out); + +int tod_to_gcrs(double jd_tdb, enum novas_accuracy accuracy, const double *in, double *out); + +int gcrs_to_mod(double jd_tdb, const double *in, double *out); + +int mod_to_gcrs(double jd_tdb, const double *in, double *out); // <================= END of SuperNOVAS API =====================> diff --git a/src/novas.c b/src/novas.c index a7fe2f87..008d876f 100644 --- a/src/novas.c +++ b/src/novas.c @@ -470,6 +470,58 @@ int gcrs_to_j2000(const double *in, double *out) { return 0; } +/** + * Transforms a rectangular equatorial (x, y, z) vector from the Geocentric Celestial + * Reference System (GCRS) to the Mean of Date (MOD) reference frame at the given epoch + * + * @param jd_tdb [day] Barycentric Dynamical Time (TT) based Julian date that defines the + * output epoch. Typically it does not require much precision, and Julian + * dates in other time measures will be unlikely to affect the result + * @param in GCRS Input (x, y, z) position or velocity vector + * @param[out] out Output position or velocity 3-vector in the Mean wquinox of Date coordinate + * frame. It can be the same vector as the input. + * @return 0 if successful, or -1 if either of the vector arguments is NULL. + * + * @sa mod_to_gcrs() + * @sa gcrs_to_tod() + * + * @since 1.2 + * @author Attila Kovacs + */ +int gcrs_to_mod(double jd_tdb, const double *in, double *out) { + static const char *fn = "gcrs_to_tod [internal]"; + prop_error(fn, frame_tie(in, ICRS_TO_J2000, out), 0); + prop_error(fn, precession(NOVAS_JD_J2000, out, jd_tdb, out), 0); + return 0; +} + +/** + * Transforms a rectangular equatorial (x, y, z) vector from Mean of Date (MOD) reference + * frame at the given epoch to the Geocentric Celestial Reference System(GCRS) + * + * @param jd_tdb [day] Barycentric Dynamical Time (TDB) based Julian date that defines the + * input epoch. Typically it does not require much precision, and Julian dates + * in other time measures will be unlikely to affect the result + * @param in Input (x, y, z) position or velocity 3-vector in the Mean equinox of Date + * coordinate frame. + * @param[out] out Output GCRS position or velocity vector. It can be the same vector as the + * input. + * @return 0 if successful, or -1 if either of the vector arguments is NULL. + * + * @sa gcrs_to_mod() + * @sa tod_to_gcrs() + * + * @since 1.2 + * @author Attila Kovacs + */ +int mod_to_gcrs(double jd_tdb, const double *in, double *out) { + static const char *fn = "tod_to_gcrs [internal]"; + prop_error(fn, precession(jd_tdb, in, NOVAS_JD_J2000, out), 0); + prop_error(fn, frame_tie(out, J2000_TO_ICRS, out), 0); + return 0; +} + + /** * Transforms a rectangular equatorial (x, y, z) vector from the Geocentric Celestial * Reference System (GCRS) to the True of Date (TOD) reference frame at the given epoch @@ -487,10 +539,10 @@ int gcrs_to_j2000(const double *in, double *out) { * @sa tod_to_gcrs() * @sa j2000_to_tod() * - * @since 1.0 + * @since 1.2 * @author Attila Kovacs */ -static int gcrs_to_tod(double jd_tdb, enum novas_accuracy accuracy, const double *in, double *out) { +int gcrs_to_tod(double jd_tdb, enum novas_accuracy accuracy, const double *in, double *out) { static const char *fn = "gcrs_to_tod [internal]"; prop_error(fn, frame_tie(in, ICRS_TO_J2000, out), 0); prop_error(fn, j2000_to_tod(jd_tdb, accuracy, out, out), 0); @@ -516,10 +568,10 @@ static int gcrs_to_tod(double jd_tdb, enum novas_accuracy accuracy, const double * @sa tod_to_j2000() * @sa tod_to_itrs() * - * @since 1.0 + * @since 1.2 * @author Attila Kovacs */ -static int tod_to_gcrs(double jd_tdb, enum novas_accuracy accuracy, const double *in, double *out) { +int tod_to_gcrs(double jd_tdb, enum novas_accuracy accuracy, const double *in, double *out) { static const char *fn = "tod_to_gcrs [internal]"; prop_error(fn, tod_to_j2000(jd_tdb, accuracy, in, out), 0); prop_error(fn, frame_tie(out, J2000_TO_ICRS, out), 0); @@ -1678,7 +1730,7 @@ short place(double jd_tt, const object *source, const observer *location, double prop_error(fn, grav_planets(pos, pob, &planets, pos), 70); // --------------------------------------------------------------------- - // Apply light and aberration. + // Apply aberration correction. // --------------------------------------------------------------------- aberration(pos, vob, t_light, pos); } @@ -1695,8 +1747,7 @@ short place(double jd_tt, const object *source, const observer *location, double case NOVAS_MOD: { // Transform to equator and equinox of date. - gcrs_to_j2000(pos, pos); - precession(NOVAS_JD_J2000, pos, jd_tdb, pos); + gcrs_to_mod(jd_tdb, pos, pos); break; } @@ -2237,9 +2288,7 @@ short gcrs2equ(double jd_tt, enum novas_dynamical_type sys, enum novas_accuracy break; case NOVAS_DYNAMICAL_MOD: { - double pos3[3]; - frame_tie(pos1, ICRS_TO_J2000, pos3); - precession(JD_J2000, pos3, jd_tdb, pos2); + gcrs_to_mod(jd_tdb, pos1, pos2); break; } @@ -2349,7 +2398,7 @@ short sidereal_time(double jd_ut1_high, double jd_ut1_low, double ut1_to_tt, enu case EROT_ERA: { // Use 'CIO-TIO-theta' method. See Circular 179, Section 6.5.4. const double ux[3] = { 1.0, 0.0, 0.0 }; - double ra_cio, ha_eq, x[3], y[3], z[3], w1[3], w2[3], eq[3]; + double ra_cio, ha_eq, x[3], y[3], z[3], eq[3]; short ref_sys; // Obtain the basis vectors, in the GCRS, of the celestial intermediate system. @@ -2358,9 +2407,7 @@ short sidereal_time(double jd_ut1_high, double jd_ut1_low, double ut1_to_tt, enu cio_basis(jd_tdb, ra_cio, ref_sys, accuracy, x, y, z); // Compute the direction of the true equinox in the GCRS. - nutation(jd_tdb, NUTATE_TRUE_TO_MEAN, accuracy, ux, w1); - precession(jd_tdb, w1, JD_J2000, w2); - frame_tie(w2, J2000_TO_ICRS, eq); + tod_to_gcrs(jd_tdb, accuracy, ux, eq); // Compute the hour angle of the equinox wrt the TIO meridian // (near Greenwich, but passes through the CIP and TIO). @@ -2996,8 +3043,7 @@ int polar_dxdy_to_dpsideps(double jd_tt, double dx, double dy, double *dpsi, dou double dp[3] = { dx * MAS, dy * MAS, dz * MAS }; // Precess pole offset vector to mean equator and equinox of date. - frame_tie(dp, ICRS_TO_J2000, dp); - precession(JD_J2000, dp, jd_tt, dp); + gcrs_to_mod(jd_tt, dp, dp); // Compute delta-delta-psi and delta-delta-epsilon in arcseconds. if(dpsi) { @@ -5835,27 +5881,29 @@ static int change_pole(const double *in, double theta, double phi, double *out) * 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. + * @param[in, out] vec vector to change to GCRS. * @return 0 if successful, or else -1 (errno set to EINVAL) if the 'sys' * argument is invalid. * * @author Attila Kovacs * @since 1.2 */ -static int equ2gcrs(double jd_tdb, const double *in, enum novas_equator_type sys, double *out) { +static int equ2gcrs(double jd_tdb, enum novas_reference_system sys, double *vec) { switch(sys) { - case NOVAS_GCRS_EQUATOR: - memcpy(out, in, XYZ_VECTOR_SIZE); + case NOVAS_GCRS: + case NOVAS_ICRS: 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); + case NOVAS_CIRS: + return cirs_to_gcrs(jd_tdb, NOVAS_REDUCED_ACCURACY, vec, vec); + case NOVAS_J2000: + return j2000_to_gcrs(vec, vec); + case NOVAS_TOD: + return tod_to_gcrs(jd_tdb, NOVAS_REDUCED_ACCURACY, vec, vec); + case NOVAS_MOD: + return mod_to_gcrs(jd_tdb, vec, vec); default: - return novas_error(-1, EINVAL, "equ2gcrs", "invalid equator type: %d", sys); + return novas_error(-1, EINVAL, "equ2gcrs", "invalid reference system: %d", sys); } } @@ -5874,18 +5922,43 @@ 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(accuracy != NOVAS_FULL_ACCURACY && accuracy != NOVAS_REDUCED_ACCURACY) + return novas_error(-1, EINVAL, fn, "invalid accuracy: %d", accuracy); + if(sys->obl) change_pole(vec, sys->obl, sys->Omega, vec); - if(sys->plane == NOVAS_ECLIPTIC_PLANE) { - if(ecl2equ_vec(jd_tdb, sys->type, accuracy, vec, vec) != 0) - return novas_trace(fn, -1, 0); + enum novas_equator_type eq; + double jd = jd_tdb; + + switch(sys->type) { + case NOVAS_GCRS: + case NOVAS_ICRS: + eq = NOVAS_GCRS_EQUATOR; + jd = NOVAS_JD_J2000; + break; + case NOVAS_J2000: + eq = NOVAS_TRUE_EQUATOR; + jd = NOVAS_JD_J2000; + break; + case NOVAS_TOD: + case NOVAS_CIRS: + eq = NOVAS_TRUE_EQUATOR; + break; + case NOVAS_MOD: + eq = NOVAS_MEAN_EQUATOR; + break; + default: + return novas_error(-1, EINVAL, fn, "invalid reference system: %d", sys->type); + } + + ecl2equ_vec(jd, eq, accuracy, vec, vec); } 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); + prop_error(fn, equ2gcrs(jd_tdb, sys->type, vec), 0); return 0; } diff --git a/src/super.c b/src/super.c index 3376c533..f2a22a35 100644 --- a/src/super.c +++ b/src/super.c @@ -1381,6 +1381,7 @@ enum novas_planet novas_planet_for_name(const char *name) { * on appropriate ephemerides, or else on up-to-date short-term orbital elements. * * + * @param type Coordinate reference system in which `ra` and `dec` are defined (e.g. NOVAS_GCRS). * @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 @@ -1392,11 +1393,12 @@ enum novas_planet novas_planet_for_name(const char *name) { * * @sa make_orbital_object() */ -int novas_set_orbital_pole(double ra, double dec, novas_orbital_system *sys) { +int novas_set_orbsys_pole(enum novas_reference_system type, double ra, double dec, novas_orbital_system *sys) { if (!sys) - return novas_error(-1, EINVAL, "novas_set_orbital_pole", "input system is NULL"); + return novas_error(-1, EINVAL, "novas_set_obsys_pole", "input system is NULL"); sys->plane = NOVAS_EQUATORIAL_PLANE; + sys->type = type; sys->obl = remainder(90.0 - dec, 360.0); sys->Omega = remainder(15.0 * ra + 90.0, 360.0); diff --git a/test/src/test-errors.c b/test/src/test-errors.c index dca1b562..98815fca 100644 --- a/test/src/test-errors.c +++ b/test/src/test-errors.c @@ -307,6 +307,48 @@ static int test_tod_to_j2000() { return n; } +static int test_gcrs_to_tod() { + double p[3]; + int n = 0; + + if(check("gcrs_to_tod:in", -1, gcrs_to_tod(0.0, NOVAS_FULL_ACCURACY, NULL, p))) n++; + if(check("gcrs_to_tod:out", -1, gcrs_to_tod(0.0, NOVAS_FULL_ACCURACY, p, NULL))) n++; + if(check("gcrs_to_tod:accuracy", -1, gcrs_to_tod(0.0, -1, p, p))) n++; + + return n; +} + +static int test_tod_to_gcrs() { + double p[3] = {1.0}; + int n = 0; + + if(check("tod_to_gcrs:in", -1, tod_to_gcrs(0.0, NOVAS_FULL_ACCURACY, NULL, p))) n++; + if(check("tod_to_gcrs:out", -1, tod_to_gcrs(0.0, NOVAS_FULL_ACCURACY, p, NULL))) n++; + if(check("tod_to_gcrs:accuracy", -1, tod_to_gcrs(0.0, -1, p, p))) n++; + + return n; +} + +static int test_gcrs_to_mod() { + double p[3]; + int n = 0; + + if(check("gcrs_to_mod:in", -1, gcrs_to_mod(0.0, NULL, p))) n++; + if(check("gcrs_to_mod:out", -1, gcrs_to_mod(0.0, p, NULL))) n++; + + return n; +} + +static int test_mod_to_gcrs() { + double p[3] = {1.0}; + int n = 0; + + if(check("mod_to_gcrs:in", -1, mod_to_gcrs(0.0, NULL, p))) n++; + if(check("mod_to_gcrs:out", -1, mod_to_gcrs(0.0, p, NULL))) n++; + + return n; +} + static int test_gcrs_to_cirs() { double p[3] = {1.0}; int n = 0; @@ -1511,38 +1553,40 @@ static int test_orbit_posvel() { 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_obsys_pole:orbit", -1, novas_orbit_posvel(0.0, NULL, NOVAS_REDUCED_ACCURACY, pos, vel))) n++; + if(check("set_obsys_pole:pos=vel", -1, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, pos, pos))) n++; + if(check("set_obsys_pole:pos=vel:NULL", -1, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, NULL, NULL))) n++; + if(check("set_obsys_pole:accuracy:-1", -1, novas_orbit_posvel(0.0, &orbit, -1, pos, vel))) n++; + if(check("set_obsys_pole:accuracy:2", -1, novas_orbit_posvel(0.0, &orbit, 2, pos, vel))) n++; - if(check("set_orbital_pole:orbit:converge", 0, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, pos, vel))) n++; + if(check("set_obsys_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++; + if(check("set_obsys_pole:orbit:converge", -1, novas_orbit_posvel(0.0, &orbit, NOVAS_REDUCED_ACCURACY, pos, vel))) n++; + else if(check("set_obsys_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++; + if(check("set_obsys_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.type = NOVAS_REFERENCE_SYSTEMS; + if(check("set_obsys_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_REFERENCE_SYSTEMS; + if(check("set_obsys_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.type = NOVAS_GCRS; 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++; + if(check("set_obsys_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() { +static int test_set_obsys_pole() { int n = 0; - if(check("set_orbital_pole:orbit", -1, novas_set_orbital_pole(0.0, 0.0, NULL))) n++; + if(check("set_obsys_pole:orbit", -1, novas_set_orbsys_pole(NOVAS_GCRS, 0.0, 0.0, NULL))) n++; return n; } @@ -1677,9 +1721,14 @@ int main() { if(test_planet_for_name()) n++; if(test_make_orbital_object()) n++; - if(test_set_orbital_pole()) n++; + if(test_set_obsys_pole()) n++; if(test_orbit_posvel()) n++; + if(test_gcrs_to_tod()) n++; + if(test_tod_to_gcrs()) n++; + if(test_gcrs_to_mod()) n++; + if(test_mod_to_gcrs()) 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 d427fcf7..48d3609a 100644 --- a/test/src/test-super.c +++ b/test/src/test-super.c @@ -76,7 +76,7 @@ static int check_equal_pos(const double *posa, const double *posb, double tol) { if(fabs(posa[i] - posb[i]) <= tol) continue; if(isnan(posa[i]) && isnan(posb[i])) continue; - fprintf(stderr, " A[%d] = %.9g vs B[%d] = %.9g\n", i, posa[i], i, posb[i]); + fprintf(stderr, " A[%d] = %.9g vs B[%d] = %.9g (delta=%.1g)\n", i, posa[i], i, posb[i], posa[i] - posb[i]); return i + 1; } @@ -744,7 +744,43 @@ static int test_transform_inv() { return 0; } +static int test_gcrs_to_tod() { + double pos1[3] = {}, pos2[3] = {}, d; + int n = 0; + + d = novas_vlen(pos0); + + if(!is_ok("gcrs_to_tod", gcrs_to_tod(tdb, NOVAS_FULL_ACCURACY, pos0, pos1))) n++; + + gcrs_to_j2000(pos0, pos2); + j2000_to_tod(tdb, NOVAS_FULL_ACCURACY, pos2, pos2); + + if(!is_ok("gcrs_to_tod:check", check_equal_pos(pos1, pos2, 1e-9 * d))) n++; + + if(!is_ok("gcrs_to_tod:tod_to_gcrs", tod_to_gcrs(tdb, NOVAS_FULL_ACCURACY, pos1, pos2))) n++; + if(!is_ok("gcrs_to_tod:tod_to_gcrs:check", check_equal_pos(pos2, pos0, 1e-9 * d))) n++; + + return n; +} +static int test_gcrs_to_mod() { + double pos1[3] = {}, pos2[3] = {}, d; + int n = 0; + + d = novas_vlen(pos0); + + if(!is_ok("gcrs_to_mod", gcrs_to_mod(tdb, pos0, pos1))) n++; + + gcrs_to_j2000(pos0, pos2); + precession(NOVAS_JD_J2000, pos2, tdb, pos2); + + if(!is_ok("gcrs_to_mod:check", check_equal_pos(pos1, pos2, 1e-9 * d))) n++; + + if(!is_ok("gcrs_to_mod:mod_to_gcrs", mod_to_gcrs(tdb, pos1, pos2))) n++; + if(!is_ok("gcrs_to_mod:mod_to_gcrs:check", check_equal_pos(pos2, pos0, 1e-9 * d))) n++; + + return n; +} static int test_source() { int k, n = 0; @@ -782,6 +818,9 @@ static int test_source() { if(test_transform_mod_tod()) n++; if(test_transform_inv()) n++; + if(test_gcrs_to_tod()) n++; + if(test_gcrs_to_mod()) n++; + for(k = 0; k < NOVAS_REFERENCE_SYSTEMS; k++) if(test_app_hor(k)) n++; for(k = 0; k < NOVAS_REFERENCE_SYSTEMS; k++) if(test_app_geom(k)) n++; @@ -2218,6 +2257,7 @@ static int test_orbit_place() { novas_orbital_elements orbit = NOVAS_ORBIT_INIT; observer obs = {}; sky_pos pos = {}; + double p0[3] = {}, p1[3] = {}; // Nov 14 0 UTC, geocentric from JPL Horizons. double tjd = 2460628.50079861; // 0 UT as TT. @@ -2247,10 +2287,50 @@ static int test_orbit_place() { 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; -} + if(!is_ok("orbit_place", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, p0, NULL))) return 1; + equ2ecl_vec(tjd, NOVAS_GCRS_EQUATOR, NOVAS_FULL_ACCURACY, p0, p0); + + orbit.system.type = NOVAS_ICRS; + if(!is_ok("orbit_place:icrs", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, p1, NULL))) n++; + else { + equ2ecl_vec(tjd, NOVAS_GCRS_EQUATOR, NOVAS_FULL_ACCURACY, p1, p1); + if(!is_ok("orbit_place:icrs:check", check_equal_pos(p1, p0, 1e-9))) n++; + } + orbit.system.type = NOVAS_CIRS; + if(!is_ok("orbit_place:cirs", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, p1, NULL))) n++; + else { + gcrs_to_cirs(tjd, NOVAS_REDUCED_ACCURACY, p1, p1); + equ2ecl_vec(tjd, NOVAS_TRUE_EQUATOR, NOVAS_FULL_ACCURACY, p1, p1); + if(!is_ok("orbit_place:cirs:check", check_equal_pos(p1, p0, 1e-9))) n++; + } + orbit.system.type = NOVAS_J2000; + if(!is_ok("orbit_place:j2000", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, p1, NULL))) n++; + else { + gcrs_to_j2000(p1, p1); + equ2ecl_vec(NOVAS_JD_J2000, NOVAS_TRUE_EQUATOR, NOVAS_FULL_ACCURACY, p1, p1); + if(!is_ok("orbit_place:j2000:check", check_equal_pos(p1, p0, 1e-9))) n++; + } + + orbit.system.type = NOVAS_MOD; + if(!is_ok("orbit_place:mod", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, p1, NULL))) n++; + else { + gcrs_to_mod(tjd, p1, p1); + equ2ecl_vec(tjd, NOVAS_MEAN_EQUATOR, NOVAS_FULL_ACCURACY, p1, p1); + if(!is_ok("orbit_place:mod:check", check_equal_pos(p1, p0, 1e-9))) n++; + } + + orbit.system.type = NOVAS_TOD; + if(!is_ok("orbit_place:tod", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, p1, NULL))) n++; + else { + gcrs_to_tod(tjd, NOVAS_FULL_ACCURACY, p1, p1); + equ2ecl_vec(tjd, NOVAS_TRUE_EQUATOR, NOVAS_FULL_ACCURACY, p1, p1); + if(!is_ok("orbit_place:tod:check", check_equal_pos(p1, p0, 1e-9))) n++; + } + + return n; +} static int test_orbit_posvel_callisto() { @@ -2279,7 +2359,7 @@ static int test_orbit_posvel_callisto() { // 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); + novas_set_orbsys_pole(NOVAS_GCRS, 268.7 / 15.0, 64.8, sys); orbit.jd_tdb = NOVAS_JD_J2000; orbit.a = 1882700.0 * 1e3 / AU; @@ -2308,27 +2388,34 @@ 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; + sys->type = NOVAS_MOD; 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; + sys->type = NOVAS_TOD; 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++; + sys->type = NOVAS_CIRS; + if(!is_ok("orbit_posvel_callisto:cirs", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, pos1, NULL))) n++; + cirs_to_gcrs(tjd, NOVAS_FULL_ACCURACY, pos0, pos); + if(!is_ok("orbit_posvel_callisto:cirs:check", check_equal_pos(pos1, pos, 1e-8))) n++; + + sys->type = NOVAS_J2000; + if(!is_ok("orbit_posvel_callisto:j2000", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, pos1, NULL))) n++; + j2000_to_gcrs(pos0, pos); + if(!is_ok("orbit_posvel_callisto:j2000:check", check_equal_pos(pos1, pos, 1e-8))) n++; + return n; }