Browse Source

make use of system_2body time_step property

main
cinnaboot 4 years ago
parent
commit
11fc5fa5e6
  1. 2
      src/main.cpp
  2. 14
      src/orbits.cpp
  3. 3
      src/orbits.h
  4. 8
      tests/orbit_test.cpp

2
src/main.cpp

@ -163,7 +163,7 @@ loadScene(RenderState* rs)
void void
updateSatelliteModel(const system_2body& sys, satellite& sat) updateSatelliteModel(const system_2body& sys, satellite& sat)
{ {
sat.theta = getPropagatedTrueAnomaly(sys, sat.theta, sys.time_step); sat.theta = getPropagatedTrueAnomaly(sys, sat.theta);
sat.gamma = orbitGetFlightPathAngle(sys.ep.e, sat.theta); sat.gamma = orbitGetFlightPathAngle(sys.ep.e, sat.theta);
sat.r = orbitGetRadialDistance(sys.ep.e, sys.ep.p, sat.theta); sat.r = orbitGetRadialDistance(sys.ep.e, sys.ep.p, sat.theta);
sat.v = orbitGetVelocity(sys.epsilon, sys.body.mu, sat.r); sat.v = orbitGetVelocity(sys.epsilon, sys.body.mu, sat.r);

14
src/orbits.cpp

@ -208,14 +208,13 @@ getNextTrialValue(double err, double ecc, double test_anom, double mean_anom)
double double
getPropagatedEccAnomaly(system_2body sys, getPropagatedEccAnomaly(system_2body sys,
double initial_anom, double initial_anom)
unsigned int time_step)
{ {
double e = sys.ep.e; double e = sys.ep.e;
double E1 = getEccAnomFromTrueAnom(e, initial_anom); double E1 = getEccAnomFromTrueAnom(e, initial_anom);
double M1 = getMeanAnomFromEccAnom(E1, e); double M1 = getMeanAnomFromEccAnom(E1, e);
double n = getMeanMotion(sys.body.mu, sys.elements.a); double n = getMeanMotion(sys.body.mu, sys.elements.a);
double M2 = getPropagatedMeanAnom(M1, n, time_step); double M2 = getPropagatedMeanAnom(M1, n, sys.time_step);
double E2_1 = getInitialTrialValue(M2, e); double E2_1 = getInitialTrialValue(M2, e);
// test if guess is a solution to kepler's equation // test if guess is a solution to kepler's equation
@ -236,10 +235,13 @@ getPropagatedEccAnomaly(system_2body sys,
double double
getPropagatedTrueAnomaly(system_2body sys, getPropagatedTrueAnomaly(system_2body sys,
double initial_anom, double initial_anom)
unsigned int time_step)
{ {
double ecc_anom = getPropagatedEccAnomaly(sys, initial_anom, time_step); // NOTE: 'pause' simulation when time_step is set close to 0
if (sys.time_step < 1e-8)
return initial_anom;
double ecc_anom = getPropagatedEccAnomaly(sys, initial_anom);
return getTrueAnomFromEccAnom(sys.ep.e, ecc_anom); return getTrueAnomFromEccAnom(sys.ep.e, ecc_anom);
} }

3
src/orbits.h

@ -145,8 +145,7 @@ orbitUpdate(orbital_elements& o, double a, double e);
*/ */
double double
getPropagatedTrueAnomaly(system_2body sys, getPropagatedTrueAnomaly(system_2body sys,
double initial_anom, double initial_anom);
unsigned int time_step);
// NOTE: aka) the trajectory equation (eq. 2.45) // NOTE: aka) the trajectory equation (eq. 2.45)
double double

8
tests/orbit_test.cpp

@ -24,8 +24,8 @@ TEST_CASE("orbit propagation", "[orbits]")
double mu = 398601.68; // NOTE: gravitational parameter double mu = 398601.68; // NOTE: gravitational parameter
double r = 6378; // NOTE: body radius in km double r = 6378; // NOTE: body radius in km
double initial_anom = 260 * M_PI / 180; // NOTE: radians double initial_anom = 260 * M_PI / 180; // NOTE: radians
unsigned int time_step = 60 * 50; // NOTE: seconds
system_2body sys = systemInit(gravBodyInit(mu, r), orbitInit(a, e)); system_2body sys = systemInit(gravBodyInit(mu, r), orbitInit(a, e));
sys.time_step = 60 * 50; // NOTE: seconds
double E1 = getEccAnomFromTrueAnom(sys.ep.e, initial_anom); double E1 = getEccAnomFromTrueAnom(sys.ep.e, initial_anom);
REQUIRE_THAT(E1, WithinAbs(-0.8615, 1e-4)); REQUIRE_THAT(E1, WithinAbs(-0.8615, 1e-4));
@ -36,17 +36,17 @@ TEST_CASE("orbit propagation", "[orbits]")
double n = getMeanMotion(sys.body.mu, sys.ep.a); double n = getMeanMotion(sys.body.mu, sys.ep.a);
REQUIRE_THAT(n, WithinAbs(0.00014582, 1e-8)); REQUIRE_THAT(n, WithinAbs(0.00014582, 1e-8));
double M2 = getPropagatedMeanAnom(M1, n, time_step); double M2 = getPropagatedMeanAnom(M1, n, sys.time_step);
REQUIRE_THAT(M2, WithinAbs(0.1383, 1e-5)); REQUIRE_THAT(M2, WithinAbs(0.1383, 1e-5));
// TODO: could also test for the other trial values listed in table 4.1 // TODO: could also test for the other trial values listed in table 4.1
double E2_1 = getInitialTrialValue(M2, sys.ep.e); double E2_1 = getInitialTrialValue(M2, sys.ep.e);
REQUIRE_THAT(E2_1, WithinAbs(0.315452, 1e-5)); REQUIRE_THAT(E2_1, WithinAbs(0.315452, 1e-5));
double ecc_anom = getPropagatedEccAnomaly(sys, initial_anom, time_step); double ecc_anom = getPropagatedEccAnomaly(sys, initial_anom);
REQUIRE_THAT(ecc_anom, WithinAbs(0.481518, 1e-5)); REQUIRE_THAT(ecc_anom, WithinAbs(0.481518, 1e-5));
double true_anom = getPropagatedTrueAnomaly(sys, initial_anom, time_step); double true_anom = getPropagatedTrueAnomaly(sys, initial_anom);
REQUIRE_THAT(true_anom, WithinAbs(1.1339, 1e-4)); REQUIRE_THAT(true_anom, WithinAbs(1.1339, 1e-4));
double r2 = orbitGetRadialDistance(sys.ep.e, sys.ep.p, true_anom); double r2 = orbitGetRadialDistance(sys.ep.e, sys.ep.p, true_anom);

Loading…
Cancel
Save