diff --git a/src/orbits.cpp b/src/orbits.cpp index ef8ed5d..9ee1fee 100644 --- a/src/orbits.cpp +++ b/src/orbits.cpp @@ -52,6 +52,12 @@ ellipseInitAE(double a, double e) return ellipseInitAB(a, b); } +double +ellipseGetEccentricity(double a, double p) +{ + return sqrt(1 - p / a); +} + bool ellipseValidate(const EllipseParameters& ep) { @@ -110,12 +116,36 @@ orbitGetAngularMomentum(double p, double mu) return sqrt(mu * p); } +double +orbitGetSemiMajorAxis(double epsilon, double mu) +{ + return -1 * mu / (2 * epsilon); +} + +double +orbitGetSemiLatusRectum(double h, double mu) +{ + return pow(h, 2) / mu; +} + +double +orbitGetAngularMomentumFromStateVectors(double r, double v, double gamma) +{ + return r * v * cos(gamma); +} + double orbitGetSpecificEnergy(double a, double mu) { return -1 * mu / (2 * a); } +double +orbitGetSpecificEnergyFromStateVectors(double r, double v, double mu) +{ + return pow(v, 2) / 2 - mu / r; +} + double orbitGetVelocity(double epsilon, double mu, double r) { diff --git a/src/orbits.h b/src/orbits.h index f7dd86c..c9ed76c 100644 --- a/src/orbits.h +++ b/src/orbits.h @@ -82,6 +82,9 @@ ellipseInitAB(double a, double b); EllipseParameters ellipseInitAE(double a, double e); +double +ellipseGetEccentricity(double a, double p); + bool ellipseValidate(const EllipseParameters& ep); @@ -157,9 +160,21 @@ orbitGetRadialDistance(double e, double p, double true_anom); double orbitGetAngularMomentum(double p, double mu); +double +orbitGetSemiMajorAxis(double epsilon, double mu); + +double +orbitGetSemiLatusRectum(double h, double mu); + +double +orbitGetAngularMomentumFromStateVectors(double r, double v, double gamma); + double orbitGetSpecificEnergy(double a, double mu); +double +orbitGetSpecificEnergyFromStateVectors(double r, double v, double mu); + glm::vec2 polarToRect(double true_anom, double r); diff --git a/tests/orbit_test.cpp b/tests/orbit_test.cpp index 9935f2c..7f16133 100644 --- a/tests/orbit_test.cpp +++ b/tests/orbit_test.cpp @@ -18,9 +18,28 @@ const double EARTH_GRAVITATIONAL_PARAMETER = 398601.68; // km^3 / s^2 const double EARTH_RADIUS = 6378; // km -TEST_CASE("orbit construction", "[orbits]") +TEST_CASE("orbit determination, example 2.1", "[orbits]") { + double r = 2124 + EARTH_RADIUS; // km + double v = 7.58; // km/s + double gamma = 20 * M_PI / 180; // flight path angle + double mu = EARTH_GRAVITATIONAL_PARAMETER; + + double epsilon = orbitGetSpecificEnergyFromStateVectors(r, v, mu); + REQUIRE_THAT(epsilon, WithinAbs(-18.1549, 2e-4)); + + double h = orbitGetAngularMomentumFromStateVectors(r, v, gamma); + REQUIRE_THAT(h, WithinAbs(60558.64, 1e-2)); + + // NOTE: we're off by an order of magnitude here compared to the inputs... + double a = orbitGetSemiMajorAxis(epsilon, mu); + REQUIRE_THAT(a, WithinAbs(10977.76, 5e-1)); + + double p = orbitGetSemiLatusRectum(h, mu); + REQUIRE_THAT(p, WithinAbs(9200.57, 5e-2)); + double e = ellipseGetEccentricity(a, p); + REQUIRE_THAT(e, WithinAbs(0.4024, 1e-4)); } TEST_CASE("orbit propagation, example 4.6", "[orbits]")