Skip to content

Commit

Permalink
More tests and fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@attipaci
attipaci committed Nov 15, 2024
1 parent fe55c58 commit 027fbd6
Show file tree
Hide file tree
Showing 4 changed files with 112 additions and 15 deletions.
19 changes: 12 additions & 7 deletions include/novas.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
*/
Expand Down Expand Up @@ -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&deg; - &delta;<sub>pole</sub>)
///< (e.g. 90&deg; - &delta;<sub>pole</sub>)
double Omega; ///< [rad] relative argument of ascending node of the orbital reference plane
///< (e.g. &alpha;<sub>pole</sub> + 90&deg;)
///< (e.g. &alpha;<sub>pole</sub> + 90&deg;)
} novas_orbital_system;


Expand Down Expand Up @@ -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_<sup>3</sup>)<sup>1/2</sup> 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.
Expand Down
18 changes: 13 additions & 5 deletions src/novas.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand All @@ -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()
Expand Down
62 changes: 62 additions & 0 deletions test/src/test-errors.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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");
Expand Down
28 changes: 25 additions & 3 deletions test/src/test-super.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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
Expand All @@ -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;
Expand All @@ -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;
}

Expand Down

0 comments on commit 027fbd6

Please sign in to comment.