Skip to content

Commit

Permalink
Change orbital angle units to deg (NOVAS convention), and add precess…
Browse files Browse the repository at this point in the history 
…ion terms
  • Loading branch information
@attipaci
attipaci committed Nov 15, 2024
1 parent 97a1b93 commit e8aaec3
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 32 deletions.
14 changes: 8 additions & 6 deletions include/novas.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -708,12 +708,14 @@ 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, 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,
///< or 2&pi;/T, where T is orbital period in days.
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 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.
double node_period; ///< [day] Precession period of the ascending node, if known.
} novas_orbital_elements;

/**
Expand Down
34 changes: 23 additions & 11 deletions src/novas.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5799,8 +5799,8 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig
* 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 theta [deg] polar angle of original pole in the new system
* @param phi [deg] 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
*
Expand All @@ -5812,6 +5812,9 @@ static int change_pole(const double *in, double theta, double phi, double *out)
y = in[1];
z = in[2];

theta *= DEGREE;
phi *= DEGREE;

double ca = cos(phi);
double sa = sin(phi);
double cb = cos(theta);
Expand Down Expand Up @@ -5903,7 +5906,7 @@ static int orbit2gcrs(double jd_tdb, const novas_orbital_system *sys, enum novas
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 M, E, nu, r;
double dt, M, E, nu, r, omega, Omega;
double cO, sO, ci, si, co, so;
double xx, yx, zx, xy, yy, zy;
int i = novas_inv_max_iter;
Expand All @@ -5914,7 +5917,8 @@ 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 = remainder(orbit->M0 + orbit->n * (jd_tdb - orbit->jd_tdb), TWOPI);
dt = (jd_tdb - orbit->jd_tdb);
E = M = remainder(orbit->M0 + orbit->n * dt, 360.0) * DEGREE;

// Iteratively determine E, using Newton-Raphson method...
while(--i >= 0) {
Expand All @@ -5933,13 +5937,21 @@ int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orbit, enum
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));

omega = orbit->omega * DEGREE;
if(orbit->apsis_period > 0.0)
omega += TWOPI * remainder(dt / orbit->apsis_period, 1.0);

Omega = orbit->Omega * DEGREE;
if(orbit->node_period > 0.0)
Omega += TWOPI * remainder(dt / orbit->node_period, 1.0);

// pos = Rz(-Omega) . Rx(-i) . Rz(-omega) . orb
cO = cos(orbit->Omega);
sO = sin(orbit->Omega);
ci = cos(orbit->i);
si = sin(orbit->i);
co = cos(orbit->omega);
so = sin(orbit->omega);
cO = cos(Omega);
sO = sin(Omega);
ci = cos(orbit->i * DEGREE);
si = sin(orbit->i * DEGREE);
co = cos(omega);
so = sin(omega);

// Rotation matrix
// See https://en.wikipedia.org/wiki/Euler_angles
Expand All @@ -5965,7 +5977,7 @@ int novas_orbit_posvel(double jd_tdb, const novas_orbital_elements *orbit, enum
}

if(vel) {
double v = orbit->n * orbit->a * orbit->a / r; // [AU/day]
double v = orbit->n * DEGREE * orbit->a * orbit->a / r; // [AU/day]
double x = -v * sin(E);
double y = v * sqrt(1.0 - orbit->e * orbit->e) * cos(E);

Expand Down
4 changes: 2 additions & 2 deletions src/super.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1395,8 +1395,8 @@ int novas_set_orbital_pole(double ra, double dec, novas_orbital_system *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);
sys->obl = remainder(90.0 - dec, 360.0);
sys->Omega = remainder(15.0 * ra + 90.0, 360.0);

return 0;
}
25 changes: 12 additions & 13 deletions test/src/test-super.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2231,18 +2231,18 @@ static int test_orbit_place() {
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;
orbit.i = 10.5879;
orbit.omega = 73.28579;
orbit.Omega = 80.25414;
orbit.M0 = 145.84905;
orbit.n = 0.21418047;

make_observer_at_geocenter(&obs);
make_orbital_object("Ceres", -1, &orbit, &ceres);

if(!is_ok("orbit_place", place(tjd, &ceres, &obs, ut12tt, NOVAS_TOD, NOVAS_REDUCED_ACCURACY, &pos))) return 1;

if(!is_equal("orbit_place:ra", pos.ra, RA0, 1e-5)) n++;
if(!is_equal("orbit_place:ra", pos.ra, RA0, 1e-5 / cos(DEC0 * DEGREE))) 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++;
Expand All @@ -2262,7 +2262,6 @@ static int test_orbit_posvel_callisto() {
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;
Expand All @@ -2281,15 +2280,15 @@ static int test_orbit_posvel_callisto() {
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.omega = 43.8;
orbit.M0 = 87.4;
orbit.i = 0.3;
orbit.Omega = 309.1;
orbit.n = TWOPI / 16.690440;
orbit.apsis_period = 277.921 * 365.25;
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;

Expand Down

0 comments on commit e8aaec3

Please sign in to comment.