9 changed files with 0 additions and 1238 deletions
@ -1,155 +0,0 @@ |
|||||||
#include <catch2/catch_test_macros.hpp> |
|
||||||
#include <catch2/matchers/catch_matchers_floating_point.hpp> |
|
||||||
#include "../src/orbital_mechanics.h" |
|
||||||
#include <cmath> |
|
||||||
|
|
||||||
TEST_CASE("Barker's equation - zero mean anomaly", "[barker][analytical]") { |
|
||||||
double M = 0.0; |
|
||||||
double nu = solve_barker_equation(M); |
|
||||||
double nu_expected = 0.0; |
|
||||||
REQUIRE_THAT(nu, Catch::Matchers::WithinAbs(nu_expected, 1e-15)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Barker's equation - small positive mean anomaly", "[barker][analytical]") { |
|
||||||
double M = 0.1; |
|
||||||
double nu = solve_barker_equation(M); |
|
||||||
REQUIRE(nu > 0.0); |
|
||||||
REQUIRE(nu < M_PI); |
|
||||||
double D = tan(nu / 2.0); |
|
||||||
double M_recovered = D + (D * D * D) / 3.0; |
|
||||||
REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Barker's equation - moderate positive mean anomaly", "[barker][analytical]") { |
|
||||||
double M = 1.0; |
|
||||||
double nu = solve_barker_equation(M); |
|
||||||
REQUIRE(nu > 0.0); |
|
||||||
REQUIRE(nu < M_PI); |
|
||||||
double D = tan(nu / 2.0); |
|
||||||
double M_recovered = D + (D * D * D) / 3.0; |
|
||||||
REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Barker's equation - large positive mean anomaly", "[barker][analytical]") { |
|
||||||
double M = 5.0; |
|
||||||
double nu = solve_barker_equation(M); |
|
||||||
REQUIRE(nu > 0.0); |
|
||||||
REQUIRE(nu < M_PI); |
|
||||||
double D = tan(nu / 2.0); |
|
||||||
double M_recovered = D + (D * D * D) / 3.0; |
|
||||||
REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Barker's equation - very large mean anomaly", "[barker][analytical]") { |
|
||||||
double M = 20.0; |
|
||||||
double nu = solve_barker_equation(M); |
|
||||||
REQUIRE(nu > 0.0); |
|
||||||
REQUIRE(nu < M_PI); |
|
||||||
double D = tan(nu / 2.0); |
|
||||||
double M_recovered = D + (D * D * D) / 3.0; |
|
||||||
REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-13)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Barker's equation - small negative mean anomaly", "[barker][analytical]") { |
|
||||||
double M = -0.1; |
|
||||||
double nu = solve_barker_equation(M); |
|
||||||
REQUIRE(nu < 0.0); |
|
||||||
REQUIRE(nu > -M_PI); |
|
||||||
double D = tan(nu / 2.0); |
|
||||||
double M_recovered = D + (D * D * D) / 3.0; |
|
||||||
REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Barker's equation - moderate negative mean anomaly", "[barker][analytical]") { |
|
||||||
double M = -1.0; |
|
||||||
double nu = solve_barker_equation(M); |
|
||||||
REQUIRE(nu < 0.0); |
|
||||||
REQUIRE(nu > -M_PI); |
|
||||||
double D = tan(nu / 2.0); |
|
||||||
double M_recovered = D + (D * D * D) / 3.0; |
|
||||||
REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Barker's equation - large negative mean anomaly", "[barker][analytical]") { |
|
||||||
double M = -5.0; |
|
||||||
double nu = solve_barker_equation(M); |
|
||||||
REQUIRE(nu < 0.0); |
|
||||||
REQUIRE(nu > -M_PI); |
|
||||||
double D = tan(nu / 2.0); |
|
||||||
double M_recovered = D + (D * D * D) / 3.0; |
|
||||||
REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Barker's equation - round-trip conversion", "[barker][analytical]") { |
|
||||||
std::vector<double> test_values = {-10.0, -5.0, -1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0, 5.0, 10.0}; |
|
||||||
|
|
||||||
for (double M_original : test_values) { |
|
||||||
double nu = solve_barker_equation(M_original); |
|
||||||
double D = tan(nu / 2.0); |
|
||||||
double M_recovered = D + (D * D * D) / 3.0; |
|
||||||
REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M_original, 1e-13)); |
|
||||||
} |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Barker's equation - true anomaly range", "[barker][analytical]") { |
|
||||||
for (double M = -50.0; M <= 50.0; M += 1.0) { |
|
||||||
double nu = solve_barker_equation(M); |
|
||||||
REQUIRE(nu > -M_PI * 0.99); |
|
||||||
REQUIRE(nu < M_PI * 0.99); |
|
||||||
} |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Parabolic orbit propagation using Barker's equation", "[barker][propagation]") { |
|
||||||
const double PARENT_MASS = 1.989e30; |
|
||||||
const double TIME_STEP = 3600.0; |
|
||||||
const int NUM_STEPS = 24; |
|
||||||
|
|
||||||
OrbitalElements initial; |
|
||||||
initial.semi_latus_rectum = 2.992e11; |
|
||||||
initial.eccentricity = 1.0; |
|
||||||
initial.true_anomaly = 0.0; |
|
||||||
initial.inclination = 0.0; |
|
||||||
initial.longitude_of_ascending_node = 0.0; |
|
||||||
initial.argument_of_periapsis = 0.0; |
|
||||||
|
|
||||||
Vec3 pos, vel; |
|
||||||
orbital_elements_to_cartesian(initial, PARENT_MASS, &pos, &vel); |
|
||||||
|
|
||||||
double initial_distance = vec3_magnitude(pos); |
|
||||||
double initial_velocity = vec3_magnitude(vel); |
|
||||||
double escape_velocity = sqrt(2.0 * G * PARENT_MASS / initial_distance); |
|
||||||
|
|
||||||
INFO("Initial distance: " << initial_distance / 1.496e11 << " AU"); |
|
||||||
INFO("Initial velocity: " << initial_velocity / 1000.0 << " km/s"); |
|
||||||
INFO("Escape velocity: " << escape_velocity / 1000.0 << " km/s"); |
|
||||||
|
|
||||||
REQUIRE_THAT(initial_velocity, Catch::Matchers::WithinAbs(escape_velocity, 1.0)); |
|
||||||
|
|
||||||
OrbitalElements current = initial; |
|
||||||
double total_time = 0.0; |
|
||||||
|
|
||||||
for (int step = 0; step < NUM_STEPS; step++) { |
|
||||||
OrbitalElements next = propagate_orbital_elements(current, TIME_STEP, PARENT_MASS); |
|
||||||
current = next; |
|
||||||
total_time += TIME_STEP; |
|
||||||
} |
|
||||||
|
|
||||||
Vec3 pos_final, vel_final; |
|
||||||
orbital_elements_to_cartesian(current, PARENT_MASS, &pos_final, &vel_final); |
|
||||||
|
|
||||||
double final_distance = vec3_magnitude(pos_final); |
|
||||||
double final_velocity = vec3_magnitude(vel_final); |
|
||||||
|
|
||||||
INFO("Final true anomaly: " << current.true_anomaly << " rad"); |
|
||||||
INFO("Final distance: " << final_distance / 1.496e11 << " AU"); |
|
||||||
INFO("Final velocity: " << final_velocity / 1000.0 << " km/s"); |
|
||||||
|
|
||||||
REQUIRE(final_distance > initial_distance); |
|
||||||
|
|
||||||
REQUIRE(final_velocity < initial_velocity); |
|
||||||
|
|
||||||
double final_escape_velocity = sqrt(2.0 * G * PARENT_MASS / final_distance); |
|
||||||
INFO("Final escape velocity: " << final_escape_velocity / 1000.0 << " km/s"); |
|
||||||
|
|
||||||
REQUIRE_THAT(final_velocity, Catch::Matchers::WithinAbs(final_escape_velocity, 1.0)); |
|
||||||
} |
|
||||||
@ -1,508 +0,0 @@ |
|||||||
#include <catch2/catch_test_macros.hpp> |
|
||||||
#include <catch2/matchers/catch_matchers_floating_point.hpp> |
|
||||||
#include <cmath> |
|
||||||
#include "../src/orbital_mechanics.h" |
|
||||||
#include "../src/orbital_objects.h" |
|
||||||
#include "../src/test_utilities.h" |
|
||||||
#include "../src/config_loader.h" |
|
||||||
#include "../src/simulation.h" |
|
||||||
|
|
||||||
using Catch::Matchers::WithinAbs; |
|
||||||
|
|
||||||
TEST_CASE("Cartesian to Elements - Advanced Tests", "[orbital_mechanics]") { |
|
||||||
const double G = 6.67430e-11; |
|
||||||
const double M_sun = 1.989e30; |
|
||||||
const double mu = G * M_sun; |
|
||||||
|
|
||||||
SECTION("Circular orbit conversion preserves exact circular parameters") { |
|
||||||
double r = 1.496e11; |
|
||||||
double v_circular = sqrt(mu / r); |
|
||||||
Vec3 position = {r, 0.0, 0.0}; |
|
||||||
Vec3 velocity = {0.0, v_circular, 0.0}; |
|
||||||
|
|
||||||
OrbitalElements elements = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(elements.eccentricity, WithinAbs(0.0, 1e-10)); |
|
||||||
REQUIRE_THAT(elements.semi_major_axis, WithinAbs(r, 1e3)); |
|
||||||
|
|
||||||
Vec3 converted_position, converted_velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &converted_position, &converted_velocity); |
|
||||||
|
|
||||||
REQUIRE(compare_vec3(position, converted_position, 1e3)); |
|
||||||
REQUIRE(compare_vec3(velocity, converted_velocity, 1e-3)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Near-circular orbit (e=0.001) recovers small eccentricity") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.496e11, |
|
||||||
.eccentricity = 0.001, |
|
||||||
.true_anomaly = 0.5, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.001, 1e-6)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.496e11, 1e3)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Elliptical orbit (e=0.5) preserves orbital shape") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = 0.8, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e6)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Highly elliptical orbit (e=0.95) preserves extreme eccentricity") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.95, |
|
||||||
.true_anomaly = 0.1, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.95, 1e-3)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e6)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Near-parabolic orbit (e=0.999) recovers near-escape trajectory") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.999, |
|
||||||
.true_anomaly = 0.05, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.999, 1e-3)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Parabolic orbit (e=1.0) recovers escape trajectory") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_latus_rectum = 1.0e11, |
|
||||||
.eccentricity = 1.0, |
|
||||||
.true_anomaly = 0.5, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(1.0, 1e-2)); |
|
||||||
REQUIRE_THAT(recovered.semi_latus_rectum, WithinAbs(1.0e11, 1e3)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Hyperbolic orbit (e=2.0) preserves unbound trajectory") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = -1.0e11, |
|
||||||
.eccentricity = 2.0, |
|
||||||
.true_anomaly = 0.5, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(2.0, 1e-3)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(-1.0e11, 1e6)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Highly hyperbolic orbit (e=10.0) preserves extreme unbound trajectory") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = -1.0e10, |
|
||||||
.eccentricity = 10.0, |
|
||||||
.true_anomaly = 0.8, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(10.0, 1e-3)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(-1.0e10, 1e8)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Zero inclination (i=0) preserves equatorial orbit") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.3, |
|
||||||
.true_anomaly = 0.5, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.inclination, WithinAbs(0.0, 1e-6)); |
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.3, 1e-4)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("90-degree inclination (i=pi/2) preserves polar orbit") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.2, |
|
||||||
.true_anomaly = 0.6, |
|
||||||
.inclination = M_PI / 2.0, |
|
||||||
.longitude_of_ascending_node = 0.5, |
|
||||||
.argument_of_periapsis = 0.3 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.inclination, WithinAbs(M_PI / 2.0, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.longitude_of_ascending_node, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.argument_of_periapsis, WithinAbs(0.3, 1e-4)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("180-degree inclination (i=pi) preserves retrograde orbit") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.2, |
|
||||||
.true_anomaly = 0.6, |
|
||||||
.inclination = M_PI, |
|
||||||
.longitude_of_ascending_node = 0.5, |
|
||||||
.argument_of_periapsis = 0.3 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.inclination, WithinAbs(M_PI, 1e-4)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Periapsis (nu=0) recovers true anomaly correctly") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = 0.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(0.0, 1e-6)); |
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Apoapsis (nu=pi) recovers true anomaly correctly") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = M_PI, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(M_PI, 1e-6)); |
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Quadrature point nu=pi/2 (90 deg) preserves orbital elements") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = M_PI / 2.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e6)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(M_PI / 2.0, 1e-6)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Quadrature point nu=-pi/2 (-90 deg) preserves orbital elements") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = -M_PI / 2.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e6)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(3.0 * M_PI / 2.0, 1e-6)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Quadrature point nu=3pi/2 (270 deg) preserves orbital elements") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = 3.0 * M_PI / 2.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e6)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(3.0 * M_PI / 2.0, 1e-6)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Quadrature point nu=-3pi/2 (-270 deg) preserves orbital elements") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = -3.0 * M_PI / 2.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e6)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(M_PI / 2.0, 1e-6)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Quadrature point with high eccentricity (e=0.9) preserves accuracy") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.9, |
|
||||||
.true_anomaly = M_PI / 2.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.9, 1e-3)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e7)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(M_PI / 2.0, 1e-5)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Quadrature point with low eccentricity (e=0.1) preserves accuracy") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.1, |
|
||||||
.true_anomaly = M_PI / 2.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.1, 1e-5)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e4)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(M_PI / 2.0, 1e-6)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Large true anomaly nu=5.0 rad (approx 286 deg) preserves accuracy") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = 5.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e6)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(5.0, 1e-6)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Large negative true anomaly nu=-5.0 rad (approx -286 deg) preserves accuracy") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = -5.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e6)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(1.28318530717958623, 1e-6)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Very large true anomaly nu=10.0 rad (approx 573 deg) preserves accuracy") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = 10.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e5)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(10.0 - 2.0 * M_PI, 1e-5)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Quadrature point with 3D orientation preserves all elements") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = M_PI / 2.0, |
|
||||||
.inclination = M_PI / 3.0, |
|
||||||
.longitude_of_ascending_node = M_PI / 4.0, |
|
||||||
.argument_of_periapsis = M_PI / 6.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e6)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(M_PI / 2.0, 1e-5)); |
|
||||||
REQUIRE_THAT(recovered.inclination, WithinAbs(M_PI / 3.0, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.longitude_of_ascending_node, WithinAbs(M_PI / 4.0, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.argument_of_periapsis, WithinAbs(M_PI / 6.0, 1e-4)); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Multiple quadrature points in sequence maintain accuracy") { |
|
||||||
double true_anomalies[] = {0.0, M_PI/4.0, M_PI/2.0, 3.0*M_PI/4.0, M_PI}; |
|
||||||
|
|
||||||
for (int i = 0; i < 5; i++) { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = 1.0e11, |
|
||||||
.eccentricity = 0.5, |
|
||||||
.true_anomaly = true_anomalies[i], |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, 1e-4)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, 1e6)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(true_anomalies[i], 1e-6)); |
|
||||||
} |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Hyperbolic orbit at quadrature point nu=pi/2") { |
|
||||||
OrbitalElements elements = { |
|
||||||
.semi_major_axis = -1.0e11, |
|
||||||
.eccentricity = 2.0, |
|
||||||
.true_anomaly = M_PI / 2.0, |
|
||||||
.inclination = 0.0, |
|
||||||
.longitude_of_ascending_node = 0.0, |
|
||||||
.argument_of_periapsis = 0.0 |
|
||||||
}; |
|
||||||
|
|
||||||
Vec3 position, velocity; |
|
||||||
orbital_elements_to_cartesian(elements, M_sun, &position, &velocity); |
|
||||||
|
|
||||||
OrbitalElements recovered = cartesian_to_orbital_elements(position, velocity, M_sun); |
|
||||||
|
|
||||||
REQUIRE_THAT(recovered.eccentricity, WithinAbs(2.0, 1e-3)); |
|
||||||
REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(-1.0e11, 1e6)); |
|
||||||
REQUIRE_THAT(recovered.true_anomaly, WithinAbs(M_PI / 2.0, 1e-5)); |
|
||||||
} |
|
||||||
} |
|
||||||
@ -1,34 +0,0 @@ |
|||||||
#include <catch2/catch_test_macros.hpp> |
|
||||||
#include "../src/physics.h" |
|
||||||
#include "../src/simulation.h" |
|
||||||
#include "../src/config_loader.h" |
|
||||||
#include "../src/test_utilities.h" |
|
||||||
#include <cmath> |
|
||||||
|
|
||||||
TEST_CASE("Energy conservation - Earth circular orbit", "[energy][rk4]") { |
|
||||||
const double TIME_STEP = 60.0; |
|
||||||
const double DAYS_TO_SIMULATE = 10.0; |
|
||||||
const double SECONDS_PER_DAY = 86400.0; |
|
||||||
|
|
||||||
SimulationState* sim = create_simulation(10, 0, 0, TIME_STEP); |
|
||||||
|
|
||||||
REQUIRE(load_system_config(sim, "tests/test_energy.toml")); |
|
||||||
|
|
||||||
double initial_energy = calculate_system_total_energy(sim); |
|
||||||
|
|
||||||
double total_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY; |
|
||||||
while (sim->time < total_time) { |
|
||||||
update_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
double final_energy = calculate_system_total_energy(sim); |
|
||||||
double energy_drift_percent = fabs((final_energy - initial_energy) / initial_energy) * 100.0; |
|
||||||
|
|
||||||
INFO("Initial energy: " << initial_energy << " J"); |
|
||||||
INFO("Final energy: " << final_energy << " J"); |
|
||||||
INFO("Energy drift: " << energy_drift_percent << "%"); |
|
||||||
|
|
||||||
REQUIRE(energy_drift_percent < 5.0); |
|
||||||
|
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
@ -1,27 +0,0 @@ |
|||||||
# Test Configuration: Sun + Earth (circular orbit) |
|
||||||
# Earth at 1 AU with circular orbit |
|
||||||
# Expected orbital period: ~365 days |
|
||||||
|
|
||||||
[[bodies]] |
|
||||||
name = "Sun" |
|
||||||
mass = 1.989e30 |
|
||||||
radius = 6.96e8 |
|
||||||
parent_index = -1 |
|
||||||
color = { r = 1.0, g = 1.0, b = 0.0 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 0.0, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
|
|
||||||
[[bodies]] |
|
||||||
name = "Earth" |
|
||||||
mass = 5.972e24 |
|
||||||
radius = 6.371e6 |
|
||||||
parent_index = 0 |
|
||||||
color = { r = 0.0, g = 0.5, b = 1.0 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 1.496e11, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
@ -1,205 +0,0 @@ |
|||||||
#include <catch2/catch_test_macros.hpp> |
|
||||||
#include <catch2/matchers/catch_matchers_floating_point.hpp> |
|
||||||
#include "../src/physics.h" |
|
||||||
#include "../src/simulation.h" |
|
||||||
#include "../src/config_loader.h" |
|
||||||
#include "../src/test_utilities.h" |
|
||||||
#include <cmath> |
|
||||||
|
|
||||||
const double POSITION_TOLERANCE_METERS = 10000.0; |
|
||||||
const double PERIOD_TOLERANCE_SECONDS = 600.0; |
|
||||||
|
|
||||||
TEST_CASE("Molniya orbit - position verification at multiple true anomalies", "[inclined][molniya]") { |
|
||||||
const double TIME_STEP = 60.0; |
|
||||||
const double SEMI_MAJOR_AXIS = 26540000.0; |
|
||||||
const double ECCENTRICITY = 0.74; |
|
||||||
|
|
||||||
SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); |
|
||||||
|
|
||||||
REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml")); |
|
||||||
|
|
||||||
Spacecraft* molniya = &sim->spacecraft[0]; |
|
||||||
CelestialBody* earth = &sim->bodies[0]; |
|
||||||
|
|
||||||
SECTION("Position at perigee (true_anomaly = 0)") { |
|
||||||
double expected_radius = SEMI_MAJOR_AXIS * (1.0 - ECCENTRICITY); |
|
||||||
double actual_radius = vec3_magnitude(vec3_sub(molniya->global_position, earth->global_position)); |
|
||||||
double radius_error = fabs(actual_radius - expected_radius); |
|
||||||
|
|
||||||
INFO("Expected radius at perigee: " << expected_radius << " m"); |
|
||||||
INFO("Actual radius: " << actual_radius << " m"); |
|
||||||
INFO("Error: " << radius_error << " m"); |
|
||||||
|
|
||||||
REQUIRE(radius_error < POSITION_TOLERANCE_METERS); |
|
||||||
|
|
||||||
CHECK(molniya->global_position.z != 0.0); |
|
||||||
INFO("Z-coordinate should be non-zero for inclined orbit (currently deferred)"); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Position at true_anomaly = π/2 (90°)") { |
|
||||||
molniya->orbit.true_anomaly = M_PI / 2.0; |
|
||||||
initialize_orbital_objects(sim); |
|
||||||
|
|
||||||
double expected_radius = SEMI_MAJOR_AXIS * (1.0 - ECCENTRICITY * ECCENTRICITY) / (1.0 + ECCENTRICITY * cos(M_PI / 2.0)); |
|
||||||
double actual_radius = vec3_magnitude(vec3_sub(molniya->global_position, earth->global_position)); |
|
||||||
double radius_error = fabs(actual_radius - expected_radius); |
|
||||||
|
|
||||||
INFO("Expected radius at ν=π/2: " << expected_radius << " m"); |
|
||||||
INFO("Actual radius: " << actual_radius << " m"); |
|
||||||
INFO("Error: " << radius_error << " m"); |
|
||||||
|
|
||||||
REQUIRE(radius_error < POSITION_TOLERANCE_METERS); |
|
||||||
|
|
||||||
CHECK(molniya->global_position.z != 0.0); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Position at apogee (true_anomaly = π)") { |
|
||||||
molniya->orbit.true_anomaly = M_PI; |
|
||||||
initialize_orbital_objects(sim); |
|
||||||
|
|
||||||
double expected_radius = SEMI_MAJOR_AXIS * (1.0 + ECCENTRICITY); |
|
||||||
double actual_radius = vec3_magnitude(vec3_sub(molniya->global_position, earth->global_position)); |
|
||||||
double radius_error = fabs(actual_radius - expected_radius); |
|
||||||
|
|
||||||
INFO("Expected radius at apogee: " << expected_radius << " m"); |
|
||||||
INFO("Actual radius: " << actual_radius << " m"); |
|
||||||
INFO("Error: " << radius_error << " m"); |
|
||||||
|
|
||||||
REQUIRE(radius_error < POSITION_TOLERANCE_METERS); |
|
||||||
|
|
||||||
CHECK(molniya->global_position.z != 0.0); |
|
||||||
INFO("At apogee, satellite should be at northernmost point (max z)"); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Position at true_anomaly = 3π/2 (270°)") { |
|
||||||
molniya->orbit.true_anomaly = 3.0 * M_PI / 2.0; |
|
||||||
initialize_orbital_objects(sim); |
|
||||||
|
|
||||||
double expected_radius = SEMI_MAJOR_AXIS * (1.0 - ECCENTRICITY * ECCENTRICITY) / (1.0 + ECCENTRICITY * cos(3.0 * M_PI / 2.0)); |
|
||||||
double actual_radius = vec3_magnitude(vec3_sub(molniya->global_position, earth->global_position)); |
|
||||||
double radius_error = fabs(actual_radius - expected_radius); |
|
||||||
|
|
||||||
INFO("Expected radius at ν=3π/2: " << expected_radius << " m"); |
|
||||||
INFO("Actual radius: " << actual_radius << " m"); |
|
||||||
INFO("Error: " << radius_error << " m"); |
|
||||||
|
|
||||||
REQUIRE(radius_error < POSITION_TOLERANCE_METERS); |
|
||||||
|
|
||||||
CHECK(molniya->global_position.z != 0.0); |
|
||||||
INFO("At ν=270°, satellite should be at southernmost point (min z)"); |
|
||||||
} |
|
||||||
|
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Molniya orbit - orbital period verification", "[inclined][molniya][period]") { |
|
||||||
const double TIME_STEP = 60.0; |
|
||||||
const double SECONDS_PER_HOUR = 3600.0; |
|
||||||
const double MAX_SIMULATION_HOURS = 15.0; |
|
||||||
// Relaxed tolerance for highly elliptical orbit with 60s timestep
|
|
||||||
const double MOLNIYA_PERIOD_TOLERANCE_SECONDS = 1800.0; // 30 minutes
|
|
||||||
|
|
||||||
SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); |
|
||||||
|
|
||||||
REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml")); |
|
||||||
|
|
||||||
Spacecraft* molniya = &sim->spacecraft[0]; |
|
||||||
CelestialBody* earth = &sim->bodies[0]; |
|
||||||
|
|
||||||
double semi_major_axis = molniya->orbit.semi_major_axis; |
|
||||||
double mu = G * earth->mass; |
|
||||||
double theoretical_period_seconds = 2.0 * M_PI * sqrt(pow(semi_major_axis, 3) / mu); |
|
||||||
double theoretical_period_hours = theoretical_period_seconds / SECONDS_PER_HOUR; |
|
||||||
|
|
||||||
INFO("Semi-major axis: " << semi_major_axis << " m"); |
|
||||||
INFO("Theoretical period from Kepler's 3rd law: " << theoretical_period_hours << " hours"); |
|
||||||
|
|
||||||
OrbitTracker* tracker = create_orbit_tracker_3d(0, 0.01, |
|
||||||
molniya->orbit.inclination, |
|
||||||
molniya->orbit.longitude_of_ascending_node, |
|
||||||
molniya->orbit.argument_of_periapsis); |
|
||||||
|
|
||||||
double max_time = MAX_SIMULATION_HOURS * SECONDS_PER_HOUR; |
|
||||||
while (sim->time < max_time && !tracker->orbit_completed) { |
|
||||||
update_simulation(sim); |
|
||||||
update_orbit_tracker(tracker, (CelestialBody*)molniya, earth, sim->time); |
|
||||||
} |
|
||||||
|
|
||||||
REQUIRE(tracker->orbit_completed); |
|
||||||
|
|
||||||
double measured_period_hours = tracker->time_at_completion / SECONDS_PER_HOUR; |
|
||||||
double period_error_hours = fabs(measured_period_hours - theoretical_period_hours); |
|
||||||
|
|
||||||
INFO("Measured period: " << measured_period_hours << " hours"); |
|
||||||
INFO("Period error: " << period_error_hours << " hours"); |
|
||||||
INFO("Period error: " << (period_error_hours / theoretical_period_hours * 100.0) << "%"); |
|
||||||
|
|
||||||
REQUIRE(period_error_hours * SECONDS_PER_HOUR < MOLNIYA_PERIOD_TOLERANCE_SECONDS); |
|
||||||
|
|
||||||
destroy_orbit_tracker(tracker); |
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Generic inclined orbit - moderate inclination", "[inclined][generic]") { |
|
||||||
const double TIME_STEP = 60.0; |
|
||||||
const double SEMI_MAJOR_AXIS = 10000000.0; |
|
||||||
const double ECCENTRICITY = 0.5; |
|
||||||
const double INCLINATION_DEG = 45.0; |
|
||||||
const double INCLINATION_RAD = INCLINATION_DEG * M_PI / 180.0; |
|
||||||
|
|
||||||
SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); |
|
||||||
|
|
||||||
REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml")); |
|
||||||
|
|
||||||
Spacecraft* craft = &sim->spacecraft[0]; |
|
||||||
CelestialBody* earth = &sim->bodies[0]; |
|
||||||
|
|
||||||
craft->orbit.semi_major_axis = SEMI_MAJOR_AXIS; |
|
||||||
craft->orbit.eccentricity = ECCENTRICITY; |
|
||||||
craft->orbit.true_anomaly = 0.0; |
|
||||||
craft->orbit.inclination = INCLINATION_RAD; |
|
||||||
craft->orbit.longitude_of_ascending_node = 0.0; |
|
||||||
craft->orbit.argument_of_periapsis = M_PI / 2.0; |
|
||||||
|
|
||||||
initialize_orbital_objects(sim); |
|
||||||
|
|
||||||
SECTION("Z-coordinate is non-zero for inclined orbit") { |
|
||||||
double z_position = craft->global_position.z; |
|
||||||
INFO("Z-coordinate: " << z_position << " m"); |
|
||||||
|
|
||||||
REQUIRE(z_position != 0.0); |
|
||||||
} |
|
||||||
|
|
||||||
SECTION("Position magnitude matches orbital radius") { |
|
||||||
double position_vector_mag = vec3_magnitude(craft->global_position); |
|
||||||
double orbital_radius = vec3_magnitude(vec3_sub(craft->global_position, earth->global_position)); |
|
||||||
double magnitude_error = fabs(position_vector_mag - orbital_radius); |
|
||||||
|
|
||||||
INFO("Position vector magnitude: " << position_vector_mag << " m"); |
|
||||||
INFO("Orbital radius: " << orbital_radius << " m"); |
|
||||||
INFO("Error: " << magnitude_error << " m"); |
|
||||||
|
|
||||||
REQUIRE(magnitude_error < POSITION_TOLERANCE_METERS); |
|
||||||
} |
|
||||||
|
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Inclined orbit - inclination parameter is preserved", "[inclined][config]") { |
|
||||||
const double TIME_STEP = 60.0; |
|
||||||
const double EXPECTED_INCLINATION_RAD = 1.107; |
|
||||||
const double EXPECTED_INCLINATION_DEG = EXPECTED_INCLINATION_RAD * 180.0 / M_PI; |
|
||||||
|
|
||||||
SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); |
|
||||||
|
|
||||||
REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml")); |
|
||||||
|
|
||||||
Spacecraft* molniya = &sim->spacecraft[0]; |
|
||||||
|
|
||||||
INFO("Loaded inclination: " << (molniya->orbit.inclination * 180.0 / M_PI) << " degrees"); |
|
||||||
INFO("Expected inclination: " << EXPECTED_INCLINATION_DEG << " degrees"); |
|
||||||
|
|
||||||
REQUIRE_THAT(molniya->orbit.inclination, Catch::Matchers::WithinAbs(EXPECTED_INCLINATION_RAD, 0.01)); |
|
||||||
|
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
@ -1,35 +0,0 @@ |
|||||||
# Test Configuration: Molniya Orbit |
|
||||||
# Earth as root body with highly elliptical, highly inclined satellite orbit |
|
||||||
# Molniya orbit parameters: |
|
||||||
# - Period: ~718 minutes (~12 hours) |
|
||||||
# - Eccentricity: 0.74 |
|
||||||
# - Inclination: 63.4° |
|
||||||
# - Argument of perigee: 270° (apogee at northernmost point) |
|
||||||
# - Perigee altitude: ~600 km |
|
||||||
# - Apogee altitude: ~39,700 km |
|
||||||
# - Semi-major axis: ~26,600 km |
|
||||||
|
|
||||||
[[bodies]] |
|
||||||
name = "Earth" |
|
||||||
mass = 5.972e24 |
|
||||||
radius = 6.371e6 |
|
||||||
parent_index = -1 |
|
||||||
color = { r = 0.0, g = 0.5, b = 1.0 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 0.0, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
|
|
||||||
[[spacecraft]] |
|
||||||
name = "Molniya_Satellite" |
|
||||||
mass = 1000.0 |
|
||||||
parent_index = 0 |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 26540000.0, |
|
||||||
eccentricity = 0.74, |
|
||||||
true_anomaly = 0.0, |
|
||||||
inclination = 1.107, |
|
||||||
longitude_of_ascending_node = 0.0, |
|
||||||
argument_of_periapsis = 4.71 |
|
||||||
} |
|
||||||
@ -1,102 +0,0 @@ |
|||||||
#include <catch2/catch_test_macros.hpp> |
|
||||||
#include "../src/physics.h" |
|
||||||
#include "../src/simulation.h" |
|
||||||
#include "../src/config_loader.h" |
|
||||||
#include "../src/test_utilities.h" |
|
||||||
#include <cmath> |
|
||||||
|
|
||||||
TEST_CASE("Orbital period - Earth (RK4)", "[period][rk4]") { |
|
||||||
const double TIME_STEP = 60.0; |
|
||||||
const double EXPECTED_PERIOD_DAYS = 365.0; |
|
||||||
const double SECONDS_PER_DAY = 86400.0; |
|
||||||
const double MAX_SIMULATION_DAYS = 400.0; |
|
||||||
|
|
||||||
SimulationState* sim = create_simulation(10, 0, 0, TIME_STEP); |
|
||||||
|
|
||||||
REQUIRE(load_system_config(sim, "tests/test_orbital_period.toml")); |
|
||||||
|
|
||||||
OrbitTracker* tracker = create_orbit_tracker(1); |
|
||||||
|
|
||||||
double max_time = MAX_SIMULATION_DAYS * SECONDS_PER_DAY; |
|
||||||
while (sim->time < max_time && !tracker->orbit_completed) { |
|
||||||
update_simulation(sim); |
|
||||||
update_orbit_tracker(tracker, &sim->bodies[1], &sim->bodies[0], sim->time); |
|
||||||
} |
|
||||||
|
|
||||||
REQUIRE(tracker->orbit_completed); |
|
||||||
|
|
||||||
double measured_period_days = tracker->time_at_completion / SECONDS_PER_DAY; |
|
||||||
double period_error_days = fabs(measured_period_days - EXPECTED_PERIOD_DAYS); |
|
||||||
|
|
||||||
INFO("Expected period: " << EXPECTED_PERIOD_DAYS << " days"); |
|
||||||
INFO("Measured period: " << measured_period_days << " days"); |
|
||||||
INFO("Error: " << period_error_days << " days"); |
|
||||||
|
|
||||||
REQUIRE(period_error_days < 5.0); |
|
||||||
|
|
||||||
destroy_orbit_tracker(tracker); |
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Orbital period - Mars (RK4)", "[period][rk4]") { |
|
||||||
const double TIME_STEP = 60.0; |
|
||||||
const double EXPECTED_PERIOD_DAYS = 687.0; |
|
||||||
const double SECONDS_PER_DAY = 86400.0; |
|
||||||
const double MAX_SIMULATION_DAYS = 750.0; |
|
||||||
|
|
||||||
SimulationState* sim = create_simulation(10, 0, 0, TIME_STEP); |
|
||||||
|
|
||||||
REQUIRE(load_system_config(sim, "tests/test_orbital_period.toml")); |
|
||||||
|
|
||||||
OrbitTracker* tracker = create_orbit_tracker(2); |
|
||||||
|
|
||||||
double max_time = MAX_SIMULATION_DAYS * SECONDS_PER_DAY; |
|
||||||
while (sim->time < max_time && !tracker->orbit_completed) { |
|
||||||
update_simulation(sim); |
|
||||||
update_orbit_tracker(tracker, &sim->bodies[2], &sim->bodies[0], sim->time); |
|
||||||
} |
|
||||||
|
|
||||||
REQUIRE(tracker->orbit_completed); |
|
||||||
|
|
||||||
double measured_period_days = tracker->time_at_completion / SECONDS_PER_DAY; |
|
||||||
double period_error_days = fabs(measured_period_days - EXPECTED_PERIOD_DAYS); |
|
||||||
|
|
||||||
INFO("Expected period: " << EXPECTED_PERIOD_DAYS << " days"); |
|
||||||
INFO("Measured period: " << measured_period_days << " days"); |
|
||||||
INFO("Error: " << period_error_days << " days"); |
|
||||||
|
|
||||||
REQUIRE(period_error_days < 25.0); |
|
||||||
|
|
||||||
destroy_orbit_tracker(tracker); |
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Orbit direction - prograde for zero inclination", "[direction]") { |
|
||||||
const double TIME_STEP = 60.0; |
|
||||||
const double TEST_DURATION_DAYS = 1.0; |
|
||||||
const double SECONDS_PER_DAY = 86400.0; |
|
||||||
const int STEPS = (int)(TEST_DURATION_DAYS * SECONDS_PER_DAY / TIME_STEP); |
|
||||||
|
|
||||||
SimulationState* sim = create_simulation(2, 0, 0, TIME_STEP); |
|
||||||
REQUIRE(load_system_config(sim, "tests/test_energy.toml")); |
|
||||||
|
|
||||||
CelestialBody* sun = &sim->bodies[0]; |
|
||||||
CelestialBody* earth = &sim->bodies[1]; |
|
||||||
|
|
||||||
Vec3 initial_rel_pos = vec3_sub(earth->global_position, sun->global_position); |
|
||||||
double theta_start = atan2(initial_rel_pos.y, initial_rel_pos.x); |
|
||||||
|
|
||||||
for (int i = 0; i < STEPS; i++) { |
|
||||||
update_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
Vec3 final_rel_pos = vec3_sub(earth->global_position, sun->global_position); |
|
||||||
double theta_final = atan2(final_rel_pos.y, final_rel_pos.x); |
|
||||||
|
|
||||||
INFO("Initial angle: " << theta_start << " rad"); |
|
||||||
INFO("Final angle: " << theta_final << " rad"); |
|
||||||
|
|
||||||
REQUIRE(theta_final > theta_start); |
|
||||||
|
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
@ -1,39 +0,0 @@ |
|||||||
# Test Configuration: Sun + Earth + Mars (circular orbits) |
|
||||||
# Earth at 1 AU, Mars at 1.5 AU with circular orbits |
|
||||||
# Expected orbital periods: Earth ~365 days, Mars ~687 days |
|
||||||
|
|
||||||
[[bodies]] |
|
||||||
name = "Sun" |
|
||||||
mass = 1.989e30 |
|
||||||
radius = 6.96e8 |
|
||||||
parent_index = -1 |
|
||||||
color = { r = 1.0, g = 1.0, b = 0.0 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 0.0, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
|
|
||||||
[[bodies]] |
|
||||||
name = "Earth" |
|
||||||
mass = 5.972e24 |
|
||||||
radius = 6.371e6 |
|
||||||
parent_index = 0 |
|
||||||
color = { r = 0.0, g = 0.5, b = 1.0 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 1.496e11, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
|
|
||||||
[[bodies]] |
|
||||||
name = "Mars" |
|
||||||
mass = 6.39e23 |
|
||||||
radius = 3.3895e6 |
|
||||||
parent_index = 0 |
|
||||||
color = { r = 0.8, g = 0.3, b = 0.1 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 2.244e11, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
@ -1,133 +0,0 @@ |
|||||||
#include <catch2/catch_test_macros.hpp> |
|
||||||
#include <catch2/matchers/catch_matchers_floating_point.hpp> |
|
||||||
#include "../src/physics.h" |
|
||||||
#include "../src/orbital_mechanics.h" |
|
||||||
#include <cmath> |
|
||||||
|
|
||||||
TEST_CASE("True anomaly round-trip conversion at periapsis", "[orbital_elements][true_anomaly]") { |
|
||||||
double parent_mass = 5.972e24; |
|
||||||
|
|
||||||
OrbitalElements elements = {0}; |
|
||||||
elements.semi_major_axis = 7000e3; |
|
||||||
elements.eccentricity = 0.3; |
|
||||||
elements.true_anomaly = 0.0; |
|
||||||
elements.inclination = 0.0; |
|
||||||
elements.longitude_of_ascending_node = 0.0; |
|
||||||
elements.argument_of_periapsis = 0.0; |
|
||||||
|
|
||||||
Vec3 pos, vel; |
|
||||||
orbital_elements_to_cartesian(elements, parent_mass, &pos, &vel); |
|
||||||
OrbitalElements reconstructed = cartesian_to_orbital_elements(pos, vel, parent_mass); |
|
||||||
|
|
||||||
INFO("Original true_anomaly: " << elements.true_anomaly); |
|
||||||
INFO("Reconstructed true_anomaly: " << reconstructed.true_anomaly); |
|
||||||
|
|
||||||
REQUIRE_THAT(reconstructed.true_anomaly, |
|
||||||
Catch::Matchers::WithinAbs(elements.true_anomaly, 0.01)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("True anomaly round-trip conversion at apoapsis", "[orbital_elements][true_anomaly]") { |
|
||||||
double parent_mass = 5.972e24; |
|
||||||
|
|
||||||
OrbitalElements elements = {0}; |
|
||||||
elements.semi_major_axis = 7000e3; |
|
||||||
elements.eccentricity = 0.3; |
|
||||||
elements.true_anomaly = M_PI; |
|
||||||
elements.inclination = 0.0; |
|
||||||
elements.longitude_of_ascending_node = 0.0; |
|
||||||
elements.argument_of_periapsis = 0.0; |
|
||||||
|
|
||||||
Vec3 pos, vel; |
|
||||||
orbital_elements_to_cartesian(elements, parent_mass, &pos, &vel); |
|
||||||
OrbitalElements reconstructed = cartesian_to_orbital_elements(pos, vel, parent_mass); |
|
||||||
|
|
||||||
INFO("Original true_anomaly: " << elements.true_anomaly); |
|
||||||
INFO("Reconstructed true_anomaly: " << reconstructed.true_anomaly); |
|
||||||
|
|
||||||
REQUIRE_THAT(reconstructed.true_anomaly, |
|
||||||
Catch::Matchers::WithinAbs(elements.true_anomaly, 0.01)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("True anomaly round-trip conversion at 90 degrees", "[orbital_elements][true_anomaly]") { |
|
||||||
double parent_mass = 5.972e24; |
|
||||||
|
|
||||||
OrbitalElements elements = {0}; |
|
||||||
elements.semi_major_axis = 7000e3; |
|
||||||
elements.eccentricity = 0.3; |
|
||||||
elements.true_anomaly = M_PI / 2.0; |
|
||||||
elements.inclination = 0.0; |
|
||||||
elements.longitude_of_ascending_node = 0.0; |
|
||||||
elements.argument_of_periapsis = 0.0; |
|
||||||
|
|
||||||
Vec3 pos, vel; |
|
||||||
orbital_elements_to_cartesian(elements, parent_mass, &pos, &vel); |
|
||||||
OrbitalElements reconstructed = cartesian_to_orbital_elements(pos, vel, parent_mass); |
|
||||||
|
|
||||||
INFO("Original true_anomaly: " << elements.true_anomaly); |
|
||||||
INFO("Reconstructed true_anomaly: " << reconstructed.true_anomaly); |
|
||||||
|
|
||||||
REQUIRE_THAT(reconstructed.true_anomaly, |
|
||||||
Catch::Matchers::WithinAbs(elements.true_anomaly, 0.01)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("True anomaly round-trip conversion at 270 degrees", "[orbital_elements][true_anomaly]") { |
|
||||||
double parent_mass = 5.972e24; |
|
||||||
|
|
||||||
OrbitalElements elements = {0}; |
|
||||||
elements.semi_major_axis = 7000e3; |
|
||||||
elements.eccentricity = 0.3; |
|
||||||
elements.true_anomaly = 3.0 * M_PI / 2.0; |
|
||||||
elements.inclination = 0.0; |
|
||||||
elements.longitude_of_ascending_node = 0.0; |
|
||||||
elements.argument_of_periapsis = 0.0; |
|
||||||
|
|
||||||
Vec3 pos, vel; |
|
||||||
orbital_elements_to_cartesian(elements, parent_mass, &pos, &vel); |
|
||||||
OrbitalElements reconstructed = cartesian_to_orbital_elements(pos, vel, parent_mass); |
|
||||||
|
|
||||||
INFO("Original true_anomaly: " << elements.true_anomaly); |
|
||||||
INFO("Reconstructed true_anomaly: " << reconstructed.true_anomaly); |
|
||||||
|
|
||||||
REQUIRE_THAT(reconstructed.true_anomaly, |
|
||||||
Catch::Matchers::WithinAbs(elements.true_anomaly, 0.01)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Radius at periapsis matches expected value", "[orbital_elements][sanity]") { |
|
||||||
double parent_mass = 5.972e24; |
|
||||||
|
|
||||||
OrbitalElements peri = {0}; |
|
||||||
peri.semi_major_axis = 7000e3; |
|
||||||
peri.eccentricity = 0.3; |
|
||||||
peri.true_anomaly = 0.0; |
|
||||||
|
|
||||||
Vec3 pos, vel; |
|
||||||
orbital_elements_to_cartesian(peri, parent_mass, &pos, &vel); |
|
||||||
double r_peri = vec3_magnitude(pos); |
|
||||||
double expected_peri = peri.semi_major_axis * (1.0 - peri.eccentricity); |
|
||||||
|
|
||||||
INFO("At true_anomaly=0:"); |
|
||||||
INFO(" Calculated radius: " << r_peri); |
|
||||||
INFO(" Expected: " << expected_peri); |
|
||||||
|
|
||||||
REQUIRE_THAT(r_peri, Catch::Matchers::WithinAbs(expected_peri, 1.0)); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Radius at apoapsis matches expected value", "[orbital_elements][sanity]") { |
|
||||||
double parent_mass = 5.972e24; |
|
||||||
|
|
||||||
OrbitalElements apo = {0}; |
|
||||||
apo.semi_major_axis = 7000e3; |
|
||||||
apo.eccentricity = 0.3; |
|
||||||
apo.true_anomaly = M_PI; |
|
||||||
|
|
||||||
Vec3 pos, vel; |
|
||||||
orbital_elements_to_cartesian(apo, parent_mass, &pos, &vel); |
|
||||||
double r_apo = vec3_magnitude(pos); |
|
||||||
double expected_apo = apo.semi_major_axis * (1.0 + apo.eccentricity); |
|
||||||
|
|
||||||
INFO("At true_anomaly=pi:"); |
|
||||||
INFO(" Calculated radius: " << r_apo); |
|
||||||
INFO(" Expected: " << expected_apo); |
|
||||||
|
|
||||||
REQUIRE_THAT(r_apo, Catch::Matchers::WithinAbs(expected_apo, 1.0)); |
|
||||||
} |
|
||||||
Loading…
Reference in new issue