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..dfd16367 100644
--- a/src/novas.c
+++ b/src/novas.c
@@ -5863,13 +5863,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 +5899,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..def52153 100644
--- a/test/src/test-errors.c
+++ b/test/src/test-errors.c
@@ -1491,6 +1491,65 @@ 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 +1679,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..dd6853e4 100644
--- a/test/src/test-super.c
+++ b/test/src/test-super.c
@@ -2254,7 +2254,7 @@ 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 +2271,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 +2291,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 +2306,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;
}