diff --git a/old_tests/test_barkers_equation.cpp b/old_tests/test_barkers_equation.cpp deleted file mode 100644 index fcee8cc..0000000 --- a/old_tests/test_barkers_equation.cpp +++ /dev/null @@ -1,155 +0,0 @@ -#include -#include -#include "../src/orbital_mechanics.h" -#include - -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 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)); -} diff --git a/old_tests/test_cartesian_to_elements_advanced.cpp b/old_tests/test_cartesian_to_elements_advanced.cpp deleted file mode 100644 index 7338c4f..0000000 --- a/old_tests/test_cartesian_to_elements_advanced.cpp +++ /dev/null @@ -1,508 +0,0 @@ -#include -#include -#include -#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)); - } -} diff --git a/old_tests/test_energy.cpp b/old_tests/test_energy.cpp deleted file mode 100644 index 9601109..0000000 --- a/old_tests/test_energy.cpp +++ /dev/null @@ -1,34 +0,0 @@ -#include -#include "../src/physics.h" -#include "../src/simulation.h" -#include "../src/config_loader.h" -#include "../src/test_utilities.h" -#include - -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); -} diff --git a/old_tests/test_energy.toml b/old_tests/test_energy.toml deleted file mode 100644 index 1b7dc41..0000000 --- a/old_tests/test_energy.toml +++ /dev/null @@ -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 -} diff --git a/old_tests/test_inclined_orbits.cpp b/old_tests/test_inclined_orbits.cpp deleted file mode 100644 index 87dfaa4..0000000 --- a/old_tests/test_inclined_orbits.cpp +++ /dev/null @@ -1,205 +0,0 @@ -#include -#include -#include "../src/physics.h" -#include "../src/simulation.h" -#include "../src/config_loader.h" -#include "../src/test_utilities.h" -#include - -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); -} diff --git a/old_tests/test_inclined_orbits.toml b/old_tests/test_inclined_orbits.toml deleted file mode 100644 index eba1ec1..0000000 --- a/old_tests/test_inclined_orbits.toml +++ /dev/null @@ -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 -} diff --git a/old_tests/test_orbital_period.cpp b/old_tests/test_orbital_period.cpp deleted file mode 100644 index dab10aa..0000000 --- a/old_tests/test_orbital_period.cpp +++ /dev/null @@ -1,102 +0,0 @@ -#include -#include "../src/physics.h" -#include "../src/simulation.h" -#include "../src/config_loader.h" -#include "../src/test_utilities.h" -#include - -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); -} diff --git a/old_tests/test_orbital_period.toml b/old_tests/test_orbital_period.toml deleted file mode 100644 index f8fa430..0000000 --- a/old_tests/test_orbital_period.toml +++ /dev/null @@ -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 -} diff --git a/old_tests/test_true_anomaly_roundtrip.cpp b/old_tests/test_true_anomaly_roundtrip.cpp deleted file mode 100644 index 26f282e..0000000 --- a/old_tests/test_true_anomaly_roundtrip.cpp +++ /dev/null @@ -1,133 +0,0 @@ -#include -#include -#include "../src/physics.h" -#include "../src/orbital_mechanics.h" -#include - -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)); -}