Skip to content

Commit

Permalink
Fixes, tweaks, and satellte orbit test
Browse files Browse the repository at this point in the history
  • Loading branch information
@attipaci
attipaci committed Nov 14, 2024
1 parent 76b6dbe commit 97a1b93
Show file tree
Hide file tree
Showing 4 changed files with 139 additions and 48 deletions.
32 changes: 18 additions & 14 deletions include/novas.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,6 @@
// for compiling your code, which may have relied on these, you can add '-DCOMPAT=1' to the
// compiler options
//
//
#if COMPAT
# include <stdio.h>
# include <ctype.h>
Expand Down Expand Up @@ -664,8 +663,10 @@ 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]).
double obl; ///< [rad] relative obliquity of orbital reference plane (e.g. 90&deg; - Dec<sub>pole</sub>)
double Omega; ///< [rad] relative argument of ascending node of the orbital reference plane (e.g. RA<sub>pole</sub>)
double obl; ///< [rad] relative obliquity of orbital reference plane
///< (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;)
} novas_orbital_system;


Expand All @@ -674,14 +675,14 @@ typedef struct {
#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
* Keplerian orbital elements for `NOVAS_ORBITAL_OBJECT` type. Orbital elements can be used to provide
* approximate positions for various Solar-system bodies. JPL publishes orbital elements (and their evolution)
* 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 describe 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
Expand All @@ -707,11 +708,12 @@ typedef struct {
double jd_tdb; ///< [day] Barycentri Dynamical Time (TDB) based Julian date of the parameters.
double a; ///< [AU] semi-major axis
double e; ///< eccentricity
double omega; ///< [rad] argument of periapsis / perihelion
double Omega; ///< [rad] argument of ascending node
double omega; ///< [rad] argument of periapsis / perihelion, at the reference time
double Omega; ///< [rad] argument of ascending node, at the reference time
double i; ///< [rad] inclination
double M0; ///< [rad] mean anomaly at tjd
double n; ///< [rad/day] mean daily motion, i.e. (_GM_/_a_<sup>3</sup>)<sup>1/2</sup> for the central body
double n; ///< [rad/day] mean daily motion, i.e. (_GM_/_a_<sup>3</sup>)<sup>1/2</sup> for the central body,
///< or 2&pi;/T, where T is orbital period in days.
} novas_orbital_elements;

/**
Expand Down Expand Up @@ -1406,6 +1408,8 @@ 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 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);
Expand Down
29 changes: 9 additions & 20 deletions src/novas.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5795,41 +5795,31 @@ 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.
* @return 0
*
*/
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 ca = cos(phi);
double sa = sin(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;
out[1] = sa * x + ca * cb * y - ca * sb * z;
out[2] = sb * y + cb * z;

return 0;
}
Expand Down Expand Up @@ -5859,7 +5849,6 @@ static int equ2gcrs(double jd_tdb, const double *in, enum novas_equator_type sys
}
}


/**
* Convert coordinates in an orbital system to GCRS equatorial coordinates
*
Expand Down Expand Up @@ -5925,13 +5914,13 @@ int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orbit, enum
if(pos == vel)
return novas_error(-1, EINVAL, fn, "output pos = vel (@ %p)", pos);

E = M = orbit->M0 + orbit->n * (jd_tdb - orbit->jd_tdb);
E = M = remainder(orbit->M0 + orbit->n * (jd_tdb - orbit->jd_tdb), TWOPI);

// Iteratively determine E, using Newton-Raphson method...
while(--i >= 0) {
double es = orbit->e * sin(E);
double ec = orbit->e * cos(E);
double dE = (E - es - M) / (1.0 - ec);
double esE = orbit->e * sin(E);
double ecE = orbit->e * cos(E);
double dE = (E - esE - M) / (1.0 - ecE);

E -= dE;
if(fabs(dE) < EPREC)
Expand Down
37 changes: 37 additions & 0 deletions src/super.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1363,3 +1363,40 @@ enum novas_planet novas_planet_for_name(const char *name) {
return novas_error(-1, EINVAL, fn, "No match for name: '%s'", name);
}

/**
* Sets the orientation of an orbital system using the RA and DEC coordinates of the pole
* of the Laplace (or else equatorial) plane relative to which the orbital elements are
* defined. Orbital parameters of planetary satellites normally include the R.A. and
* declination of the pole of the local Laplace plane in which the Keplerian orbital elements
* are referenced.
*
* The system will become referenced to the equatorial plane, the relative obliquity is set
* to (90&deg; - `dec`), while the argument of the ascending node ('Omega') is set to
* (90&deg; + `ra`).
*
* NOTES:
* <ol>
* <li>You should not expect much precision from the long-range orbital approximations for
* planetary satellites. For applications that require precision at any level, you should rely
* on appropriate ephemerides, or else on up-to-date short-term orbital elements.</li>
* </ol>
*
* @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
* @return 0 if successful, or else -1 (errno will be set to EINVAL) if the output `sys`
* pointer is NULL.
*
* @author Attila Kovacs
* @since 1.2
*/
int novas_set_orbital_pole(double ra, double dec, novas_orbital_system *sys) {
if (!sys)
return novas_error(-1, EINVAL, "novas_set_orbital_pole", "input system is NULL");

sys->plane = NOVAS_EQUATORIAL_PLANE;
sys->obl = remainder((90.0 - dec) * DEGREE, TWOPI);
sys->Omega = remainder((ra + 6.0) * HOURANGLE, TWOPI);

return 0;
}
89 changes: 75 additions & 14 deletions test/src/test-super.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2213,20 +2213,21 @@ static int test_planet_for_name() {
}


static int test_orbit_posvel() {
static int test_orbit_place() {
object ceres = {};
novas_orbital_elements orbit = NOVAS_ORBIT_INIT;
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
// Nov 14 0 UTC, geocentric from JPL Horizons.
double tjd = 2460628.50079861; // 0 UT as TT.
double RA0 = 19.684415;
double DEC0 = -28.62084;
double rv0 = 21.4255198; // km/s
double r = 3.32557776285144; // AU
int n = 0;

// Orbital Parameters for JD 2460600.5 from MPC
orbit.jd_tdb = 2460600.5;
orbit.a = 2.7666197;
orbit.e = 0.079184;
Expand All @@ -2239,13 +2240,71 @@ static int test_orbit_posvel() {
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;
if(!is_ok("orbit_place", 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++;
if(!is_equal("orbit_place:ra", pos.ra, RA0, 1e-5)) n++;
if(!is_equal("orbit_place:dec", pos.dec, DEC0, 1e-4)) n++;
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;
}


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;
int i;

// 2000-01-01 12 UT, geocentric from JPL Horizons.

double dist = 4.62117513332102;
double lt = 0.00577551831217194 * dist; // day
double tjd = 2451545.00079861 - lt; // 0 UT as TT, corrected or light time
double dyr;

double RA0 = 23.86983 * DEGREE;
double DEC0 = 8.59590 * DEGREE;

double dRA = (23.98606 * DEGREE - RA0) / cos(DEC0);
double dDEC = (8.64868 * DEGREE - DEC0);
int n = 0;

// Planet pos;
radec2vector(RA0 / HOURANGLE, DEC0 / DEGREE, dist, pos0);

// Callisto's parameters from JPL Horizons
// 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);

orbit.jd_tdb = NOVAS_JD_J2000;
dyr = (tjd - orbit.jd_tdb) / 365.25; //(yr) time since J2000.0

orbit.a = 1882700.0 * 1e3 / AU;
orbit.e = 0.007;
orbit.omega = remainder(43.8 * DEGREE + TWOPI * dyr / 277.921, TWOPI);
orbit.M0 = 87.4 * DEGREE;
orbit.i = 0.3 * DEGREE;
orbit.Omega = remainder(309.1 * DEGREE + TWOPI * dyr / 577.264, TWOPI);
orbit.n = TWOPI / 16.690440;

if(!is_ok("orbit_posvel_callisto", novas_orbit_posvel(tjd, &orbit, NOVAS_FULL_ACCURACY, pos, vel))) return 1;

for(i = 3 ;--i >= 0; ) pos[i] += pos0[i];
vector2radec(pos, &ra, &dec);

ra *= HOURANGLE;
dec *= DEGREE;

dra = (ra - RA0) * cos(DEC0);
ddec = (dec - DEC0);

if(!is_equal("orbit_posvel_callisto:dist", hypot(dra, ddec) / ARCSEC, hypot(dRA, dDEC) / ARCSEC, 15.0)) n++;
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++;

return n;
}
Expand Down Expand Up @@ -2314,7 +2373,9 @@ int main(int argc, char *argv[]) {
if(test_naif_to_novas_planet()) n++;

if(test_planet_for_name()) n++;
if(test_orbit_posvel()) n++;

if(test_orbit_place()) n++;
if(test_orbit_posvel_callisto()) n++;

n += test_dates();

Expand Down

0 comments on commit 97a1b93

Please sign in to comment.