tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it./opt/ballistics/src/solver.cpp
1 #include "ballistics/integrator.hpp"
2
3 #include <cstddef>
4 #include <stdexcept>
5 #include <vector>
6
7 #include "ballistics/model.hpp"
8 #include "ballistics/types.hpp"
9
10 namespace ballistics {
11
12 // ============================================================================
13 // STUB IMPLEMENTATION -- replace the bodies below with the real 3D high-order
14 // fixed-step integrator and sub-step event reporting.
15 //
16 // As shipped this compiles cleanly (so the Docker image builds) but is WRONG:
17 // * step() returns the state unchanged -- the projectile is frozen;
18 // * simulate() records only the launch state, never integrates, and reports
19 // no apex and no impact.
20 // Every flight test therefore fails. Implement `step` as one fourth-order step
21 // of the six-state system, assembling the acceleration from the projectile
22 // parameters and using axpy for state algebra. Implement `simulate` by marching
23 // from the launch state, reporting apex and ground impact as sub-step events
24 // rather than grid samples, and filling EVERY Flight field exactly per the
25 // header contract, including range as horizontal distance from launch, time as
26 // steps*dt, and the impacted / has_apex flags. Honor the dt > 0 contract.
27 // ============================================================================
28
29 State Integrator::step(const State& s, double dt) const {
30 if (!(dt > 0.0)) {
31 throw std::invalid_argument("Integrator::step: dt must be positive");
32 }
33 // TODO: replace with one fourth-order step over [0, dt].
34 (void)projectile_;
35 return s;
36 }
37
38 Flight Integrator::simulate(const State& initial, double dt,
39 std::size_t max_steps, double z_impact) const {
40 if (!(dt > 0.0)) {
41 throw std::invalid_argument(
42 "Integrator::simulate: dt must be positive");
43 }
44 (void)max_steps;
45 (void)z_impact;
46 // TODO: march, locate the apex and impact, and fill the Flight.
47 Flight f;
48 f.states.push_back(initial);
49 f.times.push_back(0.0);
50 f.steps = 0;
51 f.time = 0.0;
52 f.impacted = false;
53 f.has_apex = false;
54 return f;
55 }
56
57 } // namespace ballistics
58
/opt/ballistics/include/ballistics/types.hpp
1 #ifndef BALLISTICS_TYPES_HPP
2 #define BALLISTICS_TYPES_HPP
3
4 #include <cmath>
5
6 namespace ballistics {
7
8 /// A 3D Cartesian vector [m] / [m/s] / [m/s^2] depending on context.
9 /// Plain value type with the small amount of vector algebra the dynamics need.
10 struct Vec3 {
11 double x{0.0};
12 double y{0.0};
13 double z{0.0};
14 };
15
16 inline Vec3 operator+(const Vec3& a, const Vec3& b) {
17 return Vec3{a.x + b.x, a.y + b.y, a.z + b.z};
18 }
19 inline Vec3 operator-(const Vec3& a, const Vec3& b) {
20 return Vec3{a.x - b.x, a.y - b.y, a.z - b.z};
21 }
22 inline Vec3 operator*(double s, const Vec3& a) {
23 return Vec3{s * a.x, s * a.y, s * a.z};
24 }
25 inline double dot(const Vec3& a, const Vec3& b) {
26 return a.x * b.x + a.y * b.y + a.z * b.z;
27 }
28 /// Right-handed cross product a x b.
29 inline Vec3 cross(const Vec3& a, const Vec3& b) {
30 return Vec3{a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z,
31 a.x * b.y - a.y * b.x};
32 }
33 inline double norm(const Vec3& a) { return std::sqrt(dot(a, a)); }
34
35 /// Phase-space state of a point-mass projectile in 3D: position r = (x, y, z)
36 /// [m] with z measured upward, and velocity v = (vx, vy, vz) [m/s]. Plain value
37 /// type used throughout. The dynamics are autonomous (time does not appear
38 /// explicitly), so the integrator carries the simulation clock separately.
39 struct State {
40 Vec3 r{};
41 Vec3 v{};
42 };
43
44 /// Element-wise scaled-add over the full 6-component state: returns a + h * b.
45 /// Keep this helper here so the integrator stays a thin numerical layer.
46 /// (Support code.)
47 inline State axpy(const State& a, double h, const State& b) {
48 return State{Vec3{a.r.x + h * b.r.x, a.r.y + h * b.r.y, a.r.z + h * b.r.z},
49 Vec3{a.v.x + h * b.v.x, a.v.y + h * b.v.y, a.v.z + h * b.v.z}};
50 }
51
52 } // namespace ballistics
53
54 #endif // BALLISTICS_TYPES_HPP
55
/opt/ballistics/include/ballistics/integrator.hpp
1 #ifndef BALLISTICS_INTEGRATOR_HPP
2 #define BALLISTICS_INTEGRATOR_HPP
3
4 #include <cstddef>
5 #include <vector>
6
7 #include "ballistics/model.hpp"
8 #include "ballistics/types.hpp"
9
10 namespace ballistics {
11
12 /// Result of a 3D projectile flight simulation.
13 ///
14 /// `states[i]` is the fixed-step sample at grid time `times[i] = i * dt`, with
15 /// states[0] the launch state. The apex and impact are sub-step events (not the
16 /// nearest grid sample), reported to high accuracy.
17 struct Flight {
18 /// Grid samples; states[0] is the launch state.
19 std::vector<State> states;
20 /// Simulation time [s] of each sample (times[i] = i * dt). Parallel to
21 /// `states`, supplied so energy budgets can be integrated directly.
22 std::vector<double> times;
23
24 /// Number of grid steps actually taken (== states.size() - 1).
25 std::size_t steps{0};
26 /// Final grid time [s] = steps * dt (NOT the root-found impact time).
27 double time{0.0};
28
29 /// True iff the run terminated at a ground-plane impact (z crossing
30 /// z_impact from above). False if it stopped at the step cap instead, in
31 /// which case the impact fields below are not meaningful.
32 bool impacted{false};
33
34 /// Root-found apex (highest point: vz crosses zero downward). Valid once at
35 /// least one apex crossing has been seen.
36 State apex{};
37 double apex_time{0.0};
38 bool has_apex{false};
39
40 /// Root-found impact state and time (z == z_impact, descending). Valid iff
41 /// `impacted` is true.
42 State impact{};
43 double impact_time{0.0};
44
45 /// Horizontal range [m] from the launch point to the impact point, i.e.
46 /// sqrt(dx^2 + dy^2) measured in the ground (x, y) plane. Zero unless the
47 /// projectile impacted.
48 double range{0.0};
49 };
50
51 /// Fixed-step fourth-order integrator for the 3D projectile, augmented with
52 /// sub-step event reporting for the apex and ground impact.
53 class Integrator {
54 public:
55 explicit Integrator(const Projectile& projectile)
56 : projectile_(projectile) {}
57
58 /// Advance a single state by one fourth-order step of size `dt`. `dt` must be > 0
59 /// (throws std::invalid_argument otherwise).
60 State step(const State& s, double dt) const;
61
62 /// Simulate the flight from `initial`, taking fixed steps of size `dt`, until
63 /// EITHER the projectile crosses the impact plane z = `z_impact` while
64 /// descending OR `max_steps` steps have been taken, whichever comes first.
65 ///
66 /// The launch state is recorded as states[0]; each stepped grid state is
67 /// appended. The apex (vz: + -> -) and the impact (z crossing z_impact from
68 /// above) are located as sub-step events and stored in `apex` / `impact`
69 /// with their times; `range` is the horizontal
70 /// launch-to-impact distance. `dt` must be > 0 (throws
71 /// std::invalid_argument otherwise). Returns a Flight with
72 /// steps == states.size() - 1, time == steps * dt, and `impacted` /
73 /// `has_apex` set appropriately.
74 Flight simulate(const State& initial, double dt, std::size_t max_steps,
75 double z_impact = 0.0) const;
76
77 /// Alias for simulate(), matching exterior-ballistics "fire a shot"
78 /// terminology used by the CLI and the .shot format.
79 Flight fire(const State& initial, double dt, std::size_t max_steps,
80 double z_impact = 0.0) const {
81 return simulate(initial, dt, max_steps, z_impact);
82 }
83
84 private:
85 const Projectile& projectile_;
86 };
87
88 } // namespace ballistics
89
90 #endif // BALLISTICS_INTEGRATOR_HPP
91
/opt/ballistics/include/ballistics/model.hpp
1 #ifndef BALLISTICS_MODEL_HPP
2 #define BALLISTICS_MODEL_HPP
3
4 #include "ballistics/types.hpp"
5
6 namespace ballistics {
7
8 /// Two-layer model atmosphere. Air density decays with altitude using one scale
9 /// height below a transition altitude and a second scale height above it, with
10 /// the density kept continuous at the transition. Fully implemented support code.
11 class Atmosphere {
12 public:
13 /// Construct from sea-level density and lower-layer scale height. The
14 /// default transition is 400 m, and ordinary atmospheric cases use a shorter
15 /// upper-layer scale height; near-uniform atmospheres keep the supplied scale
16 /// height in both layers. Throws std::invalid_argument on invalid parameters.
17 Atmosphere(double rho0, double scale_height);
18
19 /// Construct with explicit lower-layer scale height, transition altitude, and
20 /// upper-layer scale height.
21 Atmosphere(double rho0, double lower_scale_height,
22 double transition_altitude, double upper_scale_height);
23
24 /// A vacuum (rho0 = 0): density is identically zero at every altitude.
25 static Atmosphere vacuum();
26
27 double rho0() const { return rho0_; }
28 double scale_height() const { return lower_scale_height_; }
29 double lower_scale_height() const { return lower_scale_height_; }
30 double transition_altitude() const { return transition_altitude_; }
31 double upper_scale_height() const { return upper_scale_height_; }
32
33 /// Air density at altitude z [m]. For rho0 = 0 this is identically 0.
34 double density(double z) const;
35
36 private:
37 double rho0_;
38 double lower_scale_height_;
39 double transition_altitude_;
40 double upper_scale_height_;
41 };
42
43 /// A point-mass projectile flying through a model atmosphere under uniform
44 /// gravity, quadratic aerodynamic drag, spin/Magnus deflection, and an optional
45 /// rotating-frame Coriolis term. The state is the full 3D phase vector s = (r, v)
46 /// with z up.
47 ///
48 /// Drag opposes the air-relative velocity, where air-relative velocity means the
49 /// projectile velocity minus the constant wind. The drag magnitude scales with
50 /// air density and the square of the air-relative speed. The spin/Magnus force is
51 /// perpendicular to the air-relative velocity and is directed by the right-handed
52 /// spin axis, so it curves the path without doing mechanical work. The Coriolis
53 /// convention is the rotating-frame acceleration from the stored omega vector.
54 ///
55 /// This class is a fully-implemented data container plus energy diagnostics; the
56 /// numerical integrator assembles the dynamics from these parameters.
57 class Projectile {
58 public:
59 /// Construct from mass [kg], drag coefficient Cd [-], reference area A
60 /// [m^2], gravity g [m/s^2], the atmosphere, the constant wind vector w
61 /// [m/s], the spin vector S [rad/s] (only its direction matters; sets the
62 /// Magnus axis), the Magnus lift coefficient Cl [-], and optionally the
63 /// rotating-frame angular-velocity vector omega [rad/s].
64 ///
65 /// Throws std::invalid_argument if mass <= 0, Cd < 0, area < 0, g < 0, or
66 /// Cl < 0.
67 Projectile(double mass, double Cd, double area, double gravity,
68 const Atmosphere& atmosphere, const Vec3& wind, const Vec3& spin,
69 double Cl, const Vec3& omega = Vec3{});
70
71 double mass() const { return mass_; }
72 double Cd() const { return Cd_; }
73 double area() const { return area_; }
74 double gravity() const { return gravity_; }
75 double Cl() const { return Cl_; }
76 const Vec3& wind() const { return wind_; }
77 const Vec3& spin() const { return spin_; }
78 const Vec3& omega() const { return omega_; }
79 const Atmosphere& atmosphere() const { return atmosphere_; }
80
81 /// Air density at altitude z [m] (delegates to the atmosphere).
82 double density(double z) const { return atmosphere_.density(z); }
83
84 /// Ground speed magnitude |v| [m/s]. Fully implemented.
85 double speed(const State& s) const { return norm(s.v); }
86
87 /// Total mechanical energy per unit mass E/m = 1/2 |v|^2 + g z. For the
88 /// no-wind case this decreases monotonically with aerodynamic dissipation.
89 /// Fully implemented.
90 double specific_energy(const State& s) const;
91
92 /// Instantaneous mechanical-energy dissipation rate per unit mass. In the
93 /// no-wind case this equals -d/dt(specific_energy). Fully implemented.
94 double drag_power(const State& s) const;
95
96 private:
97 double mass_;
98 double Cd_;
99 double area_;
100 double gravity_;
101 Atmosphere atmosphere_;
102 Vec3 wind_;
103 Vec3 spin_;
104 Vec3 omega_;
105 double Cl_;
106 };
107
108 } // namespace ballistics
109
110 #endif // BALLISTICS_MODEL_HPP
111
find /opt/ballistics -type f -not -path '*/build/*' -not -path '*/.git/*' | head -100
/opt/ballistics/CMakeLists.txt /opt/ballistics/examples/cannonball.shot /opt/ballistics/include/ballistics/integrator.hpp /opt/ballistics/include/ballistics/io.hpp /opt/ballistics/include/ballistics/model.hpp /opt/ballistics/include/ballistics/types.hpp /opt/ballistics/src/io.cpp /opt/ballistics/src/main.cpp /opt/ballistics/src/model.cpp /opt/ballistics/src/solver.cpp /opt/ballistics/tests/CMakeLists.txt /opt/ballistics/tests/test_smoke.cpp /opt/ballistics/tests/test_util.hpp [stdout] /opt/ballistics/CMakeLists.txt /opt/ballistics/examples/cannonball.shot /opt/ballistics/include/ballistics/integrator.hpp /opt/ballistics/include/ballistics/io.hpp /opt/ballistics/include/ballistics/model.hpp /opt/ballistics/include/ballistics/types.hpp /opt/ballistics/src/io.cpp /opt/ballistics/src/main.cpp /opt/ballistics/src/model.cpp /opt/ballistics/src/solver.cpp /opt/ballistics/tests/CMakeLists.txt /opt/ballistics/tests/test_smoke.cpp /opt/ballistics/tests/test_util.hpp
/opt/ballistics/src/model.cpp
1 #include "ballistics/model.hpp"
2
3 #include <cmath>
4 #include <stdexcept>
5
6 namespace ballistics {
7
8 namespace {
9
10 double default_upper_scale_height(double lower_scale_height) {
11 return lower_scale_height > 1.0e9 ? lower_scale_height : 6200.0;
12 }
13
14 } // namespace
15
16 Atmosphere::Atmosphere(double rho0, double scale_height)
17 : Atmosphere(rho0, scale_height, 400.0,
18 default_upper_scale_height(scale_height)) {}
19
20 Atmosphere::Atmosphere(double rho0, double lower_scale_height,
21 double transition_altitude, double upper_scale_height)
22 : rho0_(rho0),
23 lower_scale_height_(lower_scale_height),
24 transition_altitude_(transition_altitude),
25 upper_scale_height_(upper_scale_height) {
26 if (rho0_ < 0.0) {
27 throw std::invalid_argument("Atmosphere: rho0 must be non-negative");
28 }
29 if (!(lower_scale_height_ > 0.0)) {
30 throw std::invalid_argument(
31 "Atmosphere: lower scale height must be positive");
32 }
33 if (!(upper_scale_height_ > 0.0)) {
34 throw std::invalid_argument(
35 "Atmosphere: upper scale height must be positive");
36 }
37 if (transition_altitude_ < 0.0) {
38 throw std::invalid_argument(
39 "Atmosphere: transition altitude must be non-negative");
40 }
41 }
42
43 Atmosphere Atmosphere::vacuum() { return Atmosphere(0.0, 1.0); }
44
45 double Atmosphere::density(double z) const {
46 if (rho0_ == 0.0) return 0.0;
47 if (z <= transition_altitude_) {
48 return rho0_ * std::exp(-z / lower_scale_height_);
49 }
50 const double rho_transition =
51 rho0_ * std::exp(-transition_altitude_ / lower_scale_height_);
52 return rho_transition *
53 std::exp(-(z - transition_altitude_) / upper_scale_height_);
54 }
55
56 Projectile::Projectile(double mass, double Cd, double area, double gravity,
57 const Atmosphere& atmosphere, const Vec3& wind,
58 const Vec3& spin, double Cl, const Vec3& omega)
59 : mass_(mass),
60 Cd_(Cd),
61 area_(area),
62 gravity_(gravity),
63 atmosphere_(atmosphere),
64 wind_(wind),
65 spin_(spin),
66 omega_(omega),
67 Cl_(Cl) {
68 if (!(mass_ > 0.0)) {
69 throw std::invalid_argument("Projectile: mass must be positive");
70 }
71 if (Cd_ < 0.0) {
72 throw std::invalid_argument("Projectile: Cd must be non-negative");
73 }
74 if (area_ < 0.0) {
75 throw std::invalid_argument("Projectile: area must be non-negative");
76 }
77 if (gravity_ < 0.0) {
78 throw std::invalid_argument("Projectile: gravity must be non-negative");
79 }
80 if (Cl_ < 0.0) {
81 throw std::invalid_argument("Projectile: Cl must be non-negative");
82 }
83 }
84
85 double Projectile::specific_energy(const State& s) const {
86 const double v = speed(s);
87 return 0.5 * v * v + gravity_ * s.r.z;
88 }
89
90 double Projectile::drag_power(const State& s) const {
91 const double rho = density(s.r.z);
92 const Vec3 vrel = s.v - wind_;
93 const double q = norm(vrel);
94 const double kd = rho * Cd_ * area_ / (2.0 * mass_);
95 return kd * q * dot(s.v, vrel);
96 }
97
98 } // namespace ballistics
99
/opt/ballistics/tests/test_smoke.cpp
1 #include <stdexcept>
2
3 #include "ballistics/integrator.hpp"
4 #include "ballistics/model.hpp"
5 #include "ballistics/types.hpp"
6 #include "test_util.hpp"
7
8 using namespace ballistics;
9
10 namespace {
11
12 Projectile vacuum_projectile() {
13 return Projectile(1.0, 0.0, 0.0, 9.80665, Atmosphere::vacuum(),
14 Vec3{0.0, 0.0, 0.0}, Vec3{0.0, 0.0, 0.0}, 0.0);
15 }
16
17 } // namespace
18
19 BALLISTICS_TEST("vacuum_flight_smoke") {
20 const Projectile projectile = vacuum_projectile();
21 Integrator integrator(projectile);
22 const State launch{Vec3{0.0, 0.0, 0.0}, Vec3{40.0, 2.0, 45.0}};
23 const Flight flight = integrator.simulate(launch, 0.01, 20000, 0.0);
24
25 ballistics_test::expect(flight.impacted, "flight reports an impact");
26 ballistics_test::expect(flight.has_apex, "flight reports an apex");
27 ballistics_test::expect(flight.states.size() > 100, "many samples recorded");
28 ballistics_test::expect(flight.steps == flight.states.size() - 1,
29 "step count matches samples");
30 ballistics_test::expect(flight.times.size() == flight.states.size(),
31 "times parallel states");
32 ballistics_test::expect(flight.apex.r.z > launch.r.z, "apex is above launch");
33 ballistics_test::expect(flight.range > 0.0, "impact range is positive");
34 ballistics_test::expect(flight.impact.v.z < 0.0, "impact is descending");
35 }
36
37 BALLISTICS_TEST("step_cap_smoke") {
38 const Projectile projectile = vacuum_projectile();
39 Integrator integrator(projectile);
40 const State launch{Vec3{0.0, 0.0, 0.0}, Vec3{20.0, 0.0, 80.0}};
41 const std::size_t cap = 10;
42 const Flight flight = integrator.simulate(launch, 0.001, cap, -1.0e9);
43 ballistics_test::expect(!flight.impacted, "cap run does not impact");
44 ballistics_test::expect(flight.steps == cap, "cap run honors max_steps");
45 ballistics_test::expect(flight.states.back().r.z > flight.states.front().r.z,
46 "cap run advances upward");
47 }
48
49 BALLISTICS_TEST("nonpositive_dt_throws") {
50 const Projectile projectile = vacuum_projectile();
51 Integrator integrator(projectile);
52 const State launch{};
53 bool threw_step = false;
54 bool threw_sim = false;
55 try {
56 (void)integrator.step(launch, 0.0);
57 } catch (const std::invalid_argument&) {
58 threw_step = true;
59 }
60 try {
61 (void)integrator.simulate(launch, -0.01, 1, 0.0);
62 } catch (const std::invalid_argument&) {
63 threw_sim = true;
64 }
65 ballistics_test::expect(threw_step, "step rejects dt <= 0");
66 ballistics_test::expect(threw_sim, "simulate rejects dt <= 0");
67 }
68
69 BALLISTICS_TEST_MAIN()
70
/opt/ballistics/tests/test_util.hpp
1 #ifndef BALLISTICS_TEST_UTIL_HPP
2 #define BALLISTICS_TEST_UTIL_HPP
3
4 // Minimal dependency-free test harness for the ballistics suite.
5 //
6 // Each test executable defines test functions and registers them with
7 // BALLISTICS_TEST(name) { ... }. main() (provided by BALLISTICS_TEST_MAIN) runs
8 // every registered test and returns non-zero if any assertion fails, which
9 // CTest interprets as a failure.
10
11 #include <cmath>
12 #include <cstddef>
13 #include <functional>
14 #include <iostream>
15 #include <string>
16 #include <vector>
17
18 namespace ballistics_test {
19
20 struct Case {
21 std::string name;
22 std::function<void()> fn;
23 };
24
25 inline std::vector<Case>& registry() {
26 static std::vector<Case> cases;
27 return cases;
28 }
29
30 struct Registrar {
31 Registrar(const std::string& name, std::function<void()> fn) {
32 registry().push_back(Case{name, std::move(fn)});
33 }
34 };
35
36 // Thrown by an assertion failure to abort the current test case.
37 struct AssertionError {
38 std::string message;
39 };
40
41 inline void expect(bool cond, const std::string& what) {
42 if (!cond) {
43 throw AssertionError{what};
44 }
45 }
46
47 inline void expect_near(double got, double want, double tol, const std::string& what) {
48 if (std::fabs(got - want) > tol) {
49 throw AssertionError{what + " (got " + std::to_string(got) +
50 ", want " + std::to_string(want) +
51 ", tol " + std::to_string(tol) + ")"};
52 }
53 }
54
55 // Relative-error check (with an absolute floor) for quantities whose scale
56 // varies widely across cases.
57 inline void expect_rel(double got, double want, double rel_tol,
58 const std::string& what) {
59 const double denom = std::fabs(want) > 1e-300 ? std::fabs(want) : 1.0;
60 if (std::fabs(got - want) / denom > rel_tol) {
61 throw AssertionError{what + " (got " + std::to_string(got) +
62 ", want " + std::to_string(want) +
63 ", rel_tol " + std::to_string(rel_tol) + ")"};
64 }
65 }
66
67 inline int run_all() {
68 int failures = 0;
69 for (const auto& c : registry()) {
70 try {
71 c.fn();
72 std::cout << "[ PASS ] " << c.name << "\n";
73 } catch (const AssertionError& e) {
74 std::cout << "[ FAIL ] " << c.name << ": " << e.message << "\n";
75 ++failures;
76 } catch (const std::exception& e) {
77 std::cout << "[ FAIL ] " << c.name
78 << ": unexpected exception: " << e.what() << "\n";
79 ++failures;
80 } catch (...) {
81 std::cout << "[ FAIL ] " << c.name << ": unknown exception\n";
82 ++failures;
83 }
84 }
85 std::cout << "----\n"
86 << (registry().size() - failures) << "/" << registry().size()
87 << " tests passed\n";
88 return failures == 0 ? 0 : 1;
89 }
90
91 } // namespace ballistics_test
92
93 #define BALLISTICS_CONCAT_INNER(a, b) a##b
94 #define BALLISTICS_CONCAT(a, b) BALLISTICS_CONCAT_INNER(a, b)
95
96 #define BALLISTICS_TEST(NAME) \
97 static void BALLISTICS_CONCAT(ballistics_test_fn_, __LINE__)(); \
98 static ::ballistics_test::Registrar BALLISTICS_CONCAT( \
99 ballistics_test_reg_, __LINE__){ \
100 NAME, &BALLISTICS_CONCAT(ballistics_test_fn_, __LINE__)}; \
101 static void BALLISTICS_CONCAT(ballistics_test_fn_, __LINE__)()
102
103 #define BALLISTICS_TEST_MAIN() \
104 int main() { return ::ballistics_test::run_all(); }
105
106 #endif // BALLISTICS_TEST_UTIL_HPP
107
/opt/ballistics/tests/CMakeLists.txt
1 # Smoke test only. The hidden verifier supplies the authoritative grading suite.
2
3 add_executable(test_smoke test_smoke.cpp)
4 target_link_libraries(test_smoke PRIVATE ballistics)
5 target_include_directories(test_smoke PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
6 add_test(NAME smoke COMMAND test_smoke)
7
/opt/ballistics/src/io.cpp
1 #include "ballistics/io.hpp"
2
3 #include <cmath>
4 #include <fstream>
5 #include <iomanip>
6 #include <optional>
7 #include <sstream>
8 #include <stdexcept>
9 #include <string>
10
11 namespace ballistics {
12
13 namespace {
14
15 [[noreturn]] void fail(std::size_t line_no, const std::string& msg) {
16 std::ostringstream oss;
17 oss << "parse_scenario: line " << line_no << ": " << msg;
18 throw std::runtime_error(oss.str());
19 }
20
21 constexpr double kPi = 3.14159265358979323846;
22
23 } // namespace
24
25 Scenario parse_scenario(std::istream& in) {
26 std::optional<double> mass, dt, gravity, Cl;
27 std::optional<double> Cd, area;
28 std::optional<double> rho0, scale_height;
29 std::optional<double> launch_speed, elev_deg, azim_deg;
30 Vec3 wind{}, spin{}, omega{}, position{};
31 std::optional<Vec3> velocity;
32 double z_impact = 0.0;
33 std::optional<std::size_t> max_steps;
34
35 std::string line;
36 std::size_t line_no = 0;
37 while (std::getline(in, line)) {
38 ++line_no;
39 const auto hash = line.find('#');
40 if (hash != std::string::npos) line.erase(hash);
41 std::istringstream ls(line);
42 std::string tag;
43 if (!(ls >> tag)) continue;
44
45 if (tag == "mass") {
46 double v;
47 if (!(ls >> v)) fail(line_no, "expected: mass <m>");
48 mass = v;
49 } else if (tag == "drag") {
50 double c, a;
51 if (!(ls >> c >> a)) fail(line_no, "expected: drag <Cd> <area>");
52 Cd = c;
53 area = a;
54 } else if (tag == "gravity") {
55 double v;
56 if (!(ls >> v)) fail(line_no, "expected: gravity <g>");
57 gravity = v;
58 } else if (tag == "atmosphere") {
59 double r, h;
60 if (!(ls >> r >> h))
61 fail(line_no, "expected: atmosphere <rho0> <H>");
62 rho0 = r;
63 scale_height = h;
64 } else if (tag == "wind") {
65 if (!(ls >> wind.x >> wind.y >> wind.z))
66 fail(line_no, "expected: wind <wx> <wy> <wz>");
67 } else if (tag == "spin") {
68 if (!(ls >> spin.x >> spin.y >> spin.z))
69 fail(line_no, "expected: spin <sx> <sy> <sz>");
70 } else if (tag == "magnus") {
71 double v;
72 if (!(ls >> v)) fail(line_no, "expected: magnus <Cl>");
73 Cl = v;
74 } else if (tag == "omega") {
75 if (!(ls >> omega.x >> omega.y >> omega.z))
76 fail(line_no, "expected: omega <ox> <oy> <oz>");
77 } else if (tag == "position") {
78 if (!(ls >> position.x >> position.y >> position.z))
79 fail(line_no, "expected: position <x> <y> <z>");
80 } else if (tag == "launch") {
81 double sp, el, az;
82 if (!(ls >> sp >> el >> az))
83 fail(line_no,
84 "expected: launch <speed> <elev_deg> <azim_deg>");
85 launch_speed = sp;
86 elev_deg = el;
87 azim_deg = az;
88 } else if (tag == "velocity") {
89 Vec3 v;
90 if (!(ls >> v.x >> v.y >> v.z))
91 fail(line_no, "expected: velocity <vx> <vy> <vz>");
92 velocity = v;
93 } else if (tag == "ground") {
94 double v;
95 if (!(ls >> v)) fail(line_no, "expected: ground <z_impact>");
96 z_impact = v;
97 } else if (tag == "dt") {
98 double v;
99 if (!(ls >> v)) fail(line_no, "expected: dt <value>");
100 dt = v;
101 } else if (tag == "maxsteps") {
102 std::size_t v;
103 if (!(ls >> v)) fail(line_no, "expected: maxsteps <n>");
104 max_steps = v;
105 } else {
106 fail(line_no, "unknown record '" + tag + "'");
107 }
108 }
109
110 if (!mass) fail(line_no, "missing required record: mass");
111 if (!Cd || !area) fail(line_no, "missing required record: drag");
112 if (!velocity && !launch_speed)
113 fail(line_no, "missing required record: launch (or velocity)");
114 if (!dt) fail(line_no, "missing required record: dt");
115 if (!max_steps) fail(line_no, "missing required record: maxsteps");
116
117 const double g = gravity.value_or(9.80665);
118 const Atmosphere atmosphere(rho0.value_or(1.225),
119 scale_height.value_or(8500.0));
120
121 Projectile projectile(*mass, *Cd, *area, g, atmosphere, wind, spin,
122 Cl.value_or(0.0), omega);
123
124 Vec3 v0;
125 if (velocity) {
126 v0 = *velocity;
127 } else {
128 const double el = *elev_deg * kPi / 180.0;
129 const double az = *azim_deg * kPi / 180.0;
130 v0 = Vec3{*launch_speed * std::cos(el) * std::cos(az),
131 *launch_speed * std::cos(el) * std::sin(az),
132 *launch_speed * std::sin(el)};
133 }
134
135 State launch{position, v0};
136 return Scenario{projectile, launch, *dt, *max_steps, z_impact};
137 }
138
139 Scenario parse_scenario_file(const std::string& path) {
140 std::ifstream in(path);
141 if (!in) {
142 throw std::runtime_error("parse_scenario_file: cannot open '" + path +
143 "'");
144 }
145 return parse_scenario(in);
146 }
147
148 void write_report(std::ostream& out, const Projectile& projectile,
149 const Flight& flight) {
150 out << std::scientific << std::setprecision(6);
151 out << "# projectile 3D exterior-ballistics flight report\n";
152 out << "mass=" << projectile.mass() << " Cd=" << projectile.Cd()
153 << " area=" << projectile.area() << " g=" << projectile.gravity()
154 << " Cl=" << projectile.Cl() << "\n";
155 out << "rho0=" << projectile.atmosphere().rho0()
156 << " H1=" << projectile.atmosphere().lower_scale_height()
157 << " zt=" << projectile.atmosphere().transition_altitude()
158 << " H2=" << projectile.atmosphere().upper_scale_height() << "\n";
159 out << "steps=" << flight.steps << " grid_time=" << flight.time
160 << " impacted=" << (flight.impacted ? 1 : 0) << "\n";
161
162 if (flight.has_apex) {
163 out << "apex_time=" << flight.apex_time
164 << " apex_height=" << flight.apex.r.z << "\n";
165 }
166 if (flight.impacted) {
167 out << "impact_time=" << flight.impact_time
168 << " range=" << flight.range
169 << " downrange=" << flight.impact.r.x
170 << " crossrange=" << flight.impact.r.y << "\n";
171 }
172 if (!flight.states.empty()) {
173 const double e0 = projectile.specific_energy(flight.states.front());
174 const double e1 = projectile.specific_energy(flight.states.back());
175 out << "energy_lost=" << (e0 - e1) << "\n";
176 }
177 }
178
179 } // namespace ballistics
180
#include "ballistics/integrator.hpp"
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <vector>
#include "ballistics/model.hpp"
#include "ballistics/types.hpp"
namespace ballistics {
namespace {
// Right-hand side of the autonomous 6-state system s' = f(s):
// r' = v
// v' = g + drag + Magnus + Coriolis
//
// Drag opposes the air-relative velocity vrel = v - wind with the quadratic law
// a_drag = -(rho Cd A / 2m) |vrel| vrel,
// matching Projectile::drag_power (kd = rho Cd A / 2m). The Magnus/spin term is
// built the same way but rotated 90 deg about the unit spin axis s_hat,
// a_magnus = (rho Cl A / 2m) |vrel| (s_hat x vrel),
// so it is perpendicular to vrel (no mechanical work) and curves the path; only
// the direction of the spin vector matters. Coriolis is the rotating-frame term
// a_cor = -2 omega x v.
State derivative(const Projectile& p, const State& s) {
const Vec3 vrel = s.v - p.wind();
const double q = norm(vrel);
const double rho = p.density(s.r.z);
const double half_rho_A_over_m = rho * p.area() / (2.0 * p.mass());
// Gravity (z up).
Vec3 a{0.0, 0.0, -p.gravity()};
// Quadratic aerodynamic drag.
const double kd = half_rho_A_over_m * p.Cd();
a = a - (kd * q) * vrel;
// Spin/Magnus deflection (direction-only spin axis), work-free.
const double spin_mag = norm(p.spin());
if (p.Cl() > 0.0 && spin_mag > 0.0 && q > 0.0) {
const Vec3 s_hat = (1.0 / spin_mag) * p.spin();
const double kl = half_rho_A_over_m * p.Cl();
a = a + (kl * q) * cross(s_hat, vrel);
}
// Rotating-frame Coriolis acceleration.
const Vec3 om = p.omega();
if (om.x != 0.0 || om.y != 0.0 || om.z != 0.0) {
a = a - 2.0 * cross(om, s.v);
}
return State{s.v, a};
}
} // namespace
State Integrator::step(const State& s, double dt) const {
if (!(dt > 0.0)) {
throw std::invalid_argument("Integrator::step: dt must be positive");
}
// Classic fourth-order Runge-Kutta over [0, dt].
const State k1 = derivative(projectile_, s);
const State k2 = derivative(projectile_, axpy(s, 0.5 * dt, k1));
const State k3 = derivative(projectile_, axpy(s, 0.5 * dt, k2));
const State k4 = derivative(projectile_, axpy(s, dt, k3));
// increment = k1 + 2 k2 + 2 k3 + k4
State increment = axpy(k1, 2.0, k2);
increment = axpy(increment, 2.0, k3);
increment = axpy(increment, 1.0, k4);
return axpy(s, dt / 6.0, increment);
}
Flight Integrator::simulate(const State& initial, double dt,
std::size_t max_steps, double z_impact) const {
if (!(dt > 0.0)) {
throw std::invalid_argument(
"Integrator::simulate: dt must be positive");
}
Flight f;
f.states.push_back(initial);
f.times.push_back(0.0);
std::size_t taken = 0;
while (taken < max_steps) {
const State prev = f.states.back();
const State next = step(prev, dt);
const double t_prev = static_cast<double>(taken) * dt;
++taken;
f.states.push_back(next);
f.times.push_back(static_cast<double>(taken) * dt);
// Apex: vz crosses zero downward. Locate the sub-step root by bisecting
// a single RK4 step of size tau in (0, dt].
if (!f.has_apex && prev.v.z > 0.0 && next.v.z <= 0.0) {
double lo = 0.0, hi = dt;
for (int it = 0; it < 100; ++it) {
const double mid = 0.5 * (lo + hi);
if (step(prev, mid).v.z > 0.0) {
lo = mid;
} else {
hi = mid;
}
}
const double tau = 0.5 * (lo + hi);
f.apex = step(prev, tau);
f.apex_time = t_prev + tau;
f.has_apex = true;
}
// Impact: z crosses z_impact from above (descending through the plane).
if (prev.r.z > z_impact && next.r.z <= z_impact) {
double lo = 0.0, hi = dt;
for (int it = 0; it < 100; ++it) {
const double mid = 0.5 * (lo + hi);
if (step(prev, mid).r.z > z_impact) {
lo = mid;
} else {
hi = mid;
}
}
const double tau = 0.5 * (lo + hi);
f.impact = step(prev, tau);
f.impact_time = t_prev + tau;
f.impacted = true;
const double dx = f.impact.r.x - initial.r.x;
const double dy = f.impact.r.y - initial.r.y;
f.range = std::sqrt(dx * dx + dy * dy);
break;
}
}
f.steps = taken;
f.time = static_cast<double>(taken) * dt;
return f;
}
} // namespace ballistics
cmake -S /opt/ballistics -B /opt/ballistics/build -DCMAKE_BUILD_TYPE=Release >/tmp/cmake.log 2>&1 && cmake --build /opt/ballistics/build --parallel >/tmp/build.log 2>&1 && cd /opt/ballistics/build && ctest --output-on-failure; echo "EXIT=$?"
Test project /opt/ballistics/build
Start 1: smoke
1/1 Test #1: smoke ............................ Passed 0.00 sec
100% tests passed, 0 tests failed out of 1
Total Test time (real) = 0.00 sec
EXIT=0
[stdout]
Test project /opt/ballistics/build
Start 1: smoke
1/1 Test #1: smoke ............................ Passed 0.00 sec
100% tests passed, 0 tests failed out of 1
Total Test time (real) = 0.00 sec
EXIT=0cat > /tmp/check.cpp <<'EOF'
#include <cstdio>
#include <cmath>
#include "ballistics/integrator.hpp"
#include "ballistics/model.hpp"
using namespace ballistics;
int main(){
Projectile p(1.0,0.0,0.0,9.80665,Atmosphere::vacuum(),Vec3{},Vec3{},0.0);
Integrator in(p);
State l{Vec3{0,0,0},Vec3{40,2,45}};
Flight f=in.simulate(l,0.01,20000,0.0);
double g=9.80665;
printf("apex_t got=%.6f want=%.6f\n",f.apex_time,45.0/g);
printf("apex_z got=%.6f want=%.6f\n",f.apex.r.z,45.0*45.0/(2*g));
printf("impact_t got=%.6f want=%.6f\n",f.impact_time,2*45.0/g);
double T=2*45.0/g, hs=std::sqrt(40.0*40+2*2);
printf("range got=%.6f want=%.6f\n",f.range,hs*T);
printf("impact vz got=%.6f want=%.6f impacted=%d has_apex=%d\n",f.impact.v.z,-45.0,f.impacted,f.has_apex);
return 0;
}
EOF
g++ -std=c++17 -I/opt/ballistics/include /tmp/check.cpp /opt/ballistics/src/model.cpp /opt/ballistics/src/solver.cpp -o /tmp/check && /tmp/checkapex_t got=4.588723 want=4.588723 apex_z got=103.246267 want=103.246267 impact_t got=9.177446 want=9.177446 range got=367.556423 want=367.556423 impact vz got=-45.000000 want=-45.000000 impacted=1 has_apex=1 [stdout] apex_t got=4.588723 want=4.588723 apex_z got=103.246267 want=103.246267 impact_t got=9.177446 want=9.177446 range got=367.556423 want=367.556423 impact vz got=-45.000000 want=-45.000000 impacted=1 has_apex=1
No file changes, the agent left the workspace untouched.
== Restoring pristine project ==
== Staging candidate solver ==
== Injecting hidden grading tests ==
== Configuring (cmake) ==
-- The CXX compiler identification is GNU 11.4.0
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Configuring done
-- Generating done
-- Build files have been written to: /tmp/tmp.MhAAARTXwd/ballistics/build_grade
== Building ==
[ 4%] Building CXX object CMakeFiles/ballistics.dir/src/model.cpp.o
[ 8%] Building CXX object CMakeFiles/ballistics.dir/src/solver.cpp.o
[ 12%] Building CXX object CMakeFiles/ballistics.dir/src/io.cpp.o
[ 16%] Linking CXX static library libballistics.a
[ 16%] Built target ballistics
[ 20%] Building CXX object CMakeFiles/ballistics_cli.dir/src/main.cpp.o
[ 25%] Building CXX object tests/CMakeFiles/test_step.dir/test_step.cpp.o
[ 29%] Building CXX object tests/CMakeFiles/test_vacuum.dir/test_vacuum.cpp.o
[ 33%] Building CXX object tests/CMakeFiles/test_drag.dir/test_drag.cpp.o
[ 37%] Building CXX object tests/CMakeFiles/test_events.dir/test_events.cpp.o
[ 41%] Building CXX object tests/CMakeFiles/test_events_dp.dir/test_events_dp.cpp.o
[ 45%] Building CXX object tests/CMakeFiles/test_lateral.dir/test_lateral.cpp.o
[ 50%] Building CXX object tests/CMakeFiles/test_book.dir/test_book.cpp.o
[ 54%] Building CXX object tests/CMakeFiles/test_consistency.dir/test_consistency.cpp.o
[ 58%] Building CXX object tests/CMakeFiles/test_energy.dir/test_energy.cpp.o
[ 62%] Linking CXX executable ballistics_cli
[ 62%] Built target ballistics_cli
[ 66%] Linking CXX executable test_events
[ 70%] Linking CXX executable test_energy
[ 75%] Linking CXX executable test_consistency
[ 79%] Linking CXX executable test_vacuum
[ 79%] Built target test_events
[ 79%] Built target test_energy
[ 79%] Built target test_consistency
[ 79%] Built target test_vacuum
[ 83%] Linking CXX executable test_drag
[ 87%] Linking CXX executable test_lateral
[ 91%] Linking CXX executable test_book
[ 91%] Built target test_drag
[ 95%] Linking CXX executable test_step
[ 95%] Built target test_lateral
[100%] Linking CXX executable test_events_dp
[100%] Built target test_book
[100%] Built target test_events_dp
[100%] Built target test_step
== Running hidden tests ==
Test project /tmp/tmp.MhAAARTXwd/ballistics/build_grade
Start 1: test_step
1/9 Test #1: test_step ........................ Passed 0.00 sec
Start 2: test_vacuum
2/9 Test #2: test_vacuum ...................... Passed 0.01 sec
Start 3: test_events
3/9 Test #3: test_events ......................***Failed 0.00 sec
[ PASS ] apex_residual_coarse_dt
[ PASS ] impact_residual_coarse_dt
[ FAIL ] degenerate_immediate_impact: immediate impact detected
----
2/3 tests passed
Start 4: test_events_dp
4/9 Test #4: test_events_dp ................... Passed 0.01 sec
Start 5: test_drag
5/9 Test #5: test_drag ........................ Passed 0.01 sec
Start 6: test_lateral
6/9 Test #6: test_lateral ..................... Passed 0.02 sec
Start 7: test_energy
7/9 Test #7: test_energy ...................... Passed 0.01 sec
Start 8: test_book
8/9 Test #8: test_book ........................***Failed 0.02 sec
[ PASS ] step_cap_bookkeeping
[ PASS ] range_is_horizontal_from_launch
[ PASS ] grid_time_not_impact_time
[ PASS ] offaxis_3d_impact_bookkeeping
[ FAIL ] immediate_impact_nonorigin_bookkeeping: immediate impact detected
[ PASS ] apex_strictly_before_impact
[ PASS ] descending_launch_has_no_apex
[ PASS ] step_cap_no_apex_no_impact
[ PASS ] simulate_rejects_nonpositive_dt
----
8/9 tests passed
Start 9: test_consistency
9/9 Test #9: test_consistency ................. Passed 0.01 sec
78% tests passed, 2 tests failed out of 9
Label Time Summary:
hidden = 0.08 sec*proc (9 tests)
Total Test time (real) = 0.08 sec
The following tests FAILED:
3 - test_events (Failed)
8 - test_book (Failed)
Errors while running CTest
FAIL: hidden tests failedReproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_d433732ead154e3a. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_d433732ead154e3a · verifier authoritative; classifier explanatory.