add test case for orbit determination

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)
 {

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);
 

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]")