projectile-drag-integrator

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/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

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

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

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

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

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

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

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

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

1	#ifndef BALLISTICS_IO_HPP
2	#define BALLISTICS_IO_HPP
3	
4	#include <cstddef>
5	#include <iosfwd>
6	#include <string>
7	
8	#include "ballistics/integrator.hpp"
9	#include "ballistics/model.hpp"
10	#include "ballistics/types.hpp"
11	
12	namespace ballistics {
13	
14	/// A complete scenario parsed from the `.shot` text format.
15	struct Scenario {
16	    Projectile projectile;
17	    State launch;
18	    double dt{0.0};
19	    std::size_t max_steps{0};
20	    double z_impact{0.0};
21	};
22	
23	/// Parse a scenario from the engine's plain-text `.shot` format.
24	///
25	/// The format is line-oriented. Blank lines and lines beginning with '#' are
26	/// ignored. Recognized records:
27	///   mass     <m>                    # projectile mass [kg]            (required)
28	///   drag     <Cd> <area>            # drag coefficient and area [m^2] (required)
29	///   gravity  <g>                    # gravity [m/s^2] (default 9.80665)
30	///   atmosphere <rho0> <H>           # sea-level density, lower scale height
31	///                                   #   (default 1.225 8500); rho0=0 vacuum
32	///   wind     <wx> <wy> <wz>         # constant wind vector [m/s] (default 0 0 0)
33	///   spin     <sx> <sy> <sz>         # spin vector [rad/s], sets Magnus axis
34	///                                   #   (default 0 0 0)
35	///   magnus   <Cl>                   # Magnus lift coefficient (default 0)
36	///   omega    <ox> <oy> <oz>         # rotating-frame angular velocity [rad/s]
37	///                                   #   (default 0 0 0)
38	///   position <x> <y> <z>            # launch position [m] (default 0 0 0)
39	///   launch   <speed> <elev_deg> <azim_deg>
40	///                                   # launch speed [m/s], elevation above the
41	///                                   #   horizon [deg], azimuth from +x toward
42	///                                   #   +y [deg]                      (required)
43	///   velocity <vx> <vy> <vz>         # explicit launch velocity [m/s] (this OR
44	///                                   #   launch)
45	///   ground   <z_impact>             # impact plane altitude [m] (default 0)
46	///   dt       <value>                # integration step [s]            (required)
47	///   maxsteps <n>                    # step cap                        (required)
48	///
49	/// Throws std::runtime_error on malformed input or a missing required record.
50	Scenario parse_scenario(std::istream& in);
51	
52	/// Convenience overload that parses from a file path.
53	Scenario parse_scenario_file(const std::string& path);
54	
55	/// Write a human-readable flight report (root-found apex / impact / range /
56	/// time of flight, lateral deflection, energy dissipated) to `out`. Fully
57	/// implemented; used by the demo CLI.
58	void write_report(std::ostream& out, const Projectile& projectile,
59	                  const Flight& flight);
60	
61	} // namespace ballistics
62	
63	#endif // BALLISTICS_IO_HPP
64

1	#include <cstdlib>
2	#include <iostream>
3	
4	#include "ballistics/integrator.hpp"
5	#include "ballistics/io.hpp"
6	#include "ballistics/model.hpp"
7	
8	// Command-line demo driver: read a .shot scenario, fire the projectile through
9	// the 3D exterior-ballistics integrator, and print
10	// the root-found apex / impact / range / time-of-flight and energy dissipated.
11	//   ballistics_cli <scenario.shot>
12	int main(int argc, char** argv) {
13	    if (argc != 2) {
14	        std::cerr << "usage: " << argv[0] << " <scenario.shot>\n";
15	        return 2;
16	    }
17	    try {
18	        ballistics::Scenario sc = ballistics::parse_scenario_file(argv[1]);
19	        ballistics::Integrator integrator(sc.projectile);
20	        ballistics::Flight flight =
21	            integrator.fire(sc.launch, sc.dt, sc.max_steps, sc.z_impact);
22	        ballistics::write_report(std::cout, sc.projectile, flight);
23	    } catch (const std::exception& ex) {
24	        std::cerr << "error: " << ex.what() << "\n";
25	        return 1;
26	    }
27	    return 0;
28	}
29

[ 25%] Building CXX object CMakeFiles/ballistics.dir/src/solver.cpp.o
[ 37%] Building CXX object CMakeFiles/ballistics.dir/src/io.cpp.o
[ 12%] Building CXX object CMakeFiles/ballistics.dir/src/model.cpp.o
[ 50%] Linking CXX static library libballistics.a
[ 50%] Built target ballistics
[ 62%] Building CXX object CMakeFiles/ballistics_cli.dir/src/main.cpp.o
[ 75%] Building CXX object tests/CMakeFiles/test_smoke.dir/test_smoke.cpp.o
[ 87%] Linking CXX executable ballistics_cli
[ 87%] Built target ballistics_cli
[100%] Linking CXX executable test_smoke
[100%] Built target test_smoke
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

[stdout]
[ 25%] Building CXX object CMakeFiles/ballistics.dir/src/solver.cpp.o
[ 37%] Building CXX object CMakeFiles/ballistics.dir/src/io.cpp.o
[ 12%] Building CXX object CMakeFiles/ballistics.dir/src/model.cpp.o
[ 50%] Linking CXX static library libballistics.a
[ 50%] Built target ballistics
[ 62%] Building CXX object CMakeFiles/ballistics_cli.dir/src/main.cpp.o
[ 75%] Building CXX object tests/CMakeFiles/test_smoke.dir/test_smoke.cpp.o
[ 87%] Linking CXX executable ballistics_cli
[ 87%] Built target ballistics_cli
[100%] Linking CXX executable test_smoke
[100%] Built target test_smoke
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

apex_t  got 4.588723  exact 4.588723
apex_z  got 103.246267  exact 103.246267
imp_t   got 9.177446  exact 9.177446
range   got 367.556423  exact 367.556423
imp_z   got 1.387779e-17 (want 0)
steps 918 states 919 time 9.1800

[stdout]
apex_t  got 4.588723  exact 4.588723
apex_z  got 103.246267  exact 103.246267
imp_t   got 9.177446  exact 9.177446
range   got 367.556423  exact 367.556423
imp_z   got 1.387779e-17 (want 0)
steps 918 states 919 time 9.1800

# A spin-stabilized round fired at 250 m/s, 35 deg elevation, due +x (azimuth 0),
# through an exponential atmosphere (sea-level rho0 = 1.225 kg/m^3, scale height
# 8500 m). A steady crosswind blows in +y and the round carries right-hand spin
# about the +z axis, so the Magnus force deflects it laterally out of the firing
# plane. Drag plus the thinning air aloft make the impact noticeably shorter
# than the drag-free range.
mass       8.0
drag       0.30 0.0078     # Cd, reference area [m^2] (~ 10 cm diameter)
gravity    9.80665
atmosphere 1.225 8500.0
wind       0.0 6.0 0.0     # 6 m/s crosswind toward +y
spin       0.0 0.0 220.0   # spin about +z (rad/s); only the axis matters
magnus     0.18            # Magnus lift coefficient
position   0.0 0.0 0.0
launch     250.0 35.0 0.0  # speed, elevation [deg], azimuth [deg]
ground     0.0
dt         0.002
maxsteps   2000000
=== REPORT ===
# projectile 3D exterior-ballistics flight report
mass=8.000000e+00 Cd=3.000000e-01 area=7.800000e-03 g=9.806650e+00 Cl=1.800000e-01
rho0=1.225000e+00 H1=8.500000e+03 zt=4.000000e+02 H2=6.200000e+03
steps=12564 grid_time=2.512800e+01 impacted=1
apex_time=1.186335e+01 apex_height=7.804948e+02
impact_time=2.512667e+01 range=3.628516e+03 downrange=3.546984e+03 crossrange=7.648704e+02
energy_lost=1.941071e+04

[stdout]
# A spin-stabilized round fired at 250 m/s, 35 deg elevation, due +x (azimuth 0),
# through an exponential atmosphere (sea-level rho0 = 1.225 kg/m^3, scale height
# 8500 m). A steady crosswind blows in +y and the round carries right-hand spin
# about the +z axis, so the Magnus force deflects it laterally out of the firing
# plane. Drag plus the thinning air aloft make the impact noticeably shorter
# than the drag-free range.
mass       8.0
drag       0.30 0.0078     # Cd, reference area [m^2] (~ 10 cm diameter)
gravity    9.80665
atmosphere 1.225 8500.0
wind       0.0 6.0 0.0     # 6 m/s crosswind toward +y
spin       0.0 0.0 220.0   # spin about +z (rad/s); only the axis matters
magnus     0.18            # Magnus lift coefficient
position   0.0 0.0 0.0
launch     250.0 35.0 0.0  # speed, elevation [deg], azimuth [deg]
ground     0.0
dt         0.002
maxsteps   2000000
=== REPORT ===
# projectile 3D exterior-ballistics flight report
mass=8.000000e+00 Cd=3.000000e-01 area=7.800000e-03 g=9.806650e+00 Cl=1.800000e-01
rho0=1.225000e+00 H1=8.500000e+03 zt=4.000000e+02 H2=6.200000e+03
steps=12564 grid_time=2.512800e+01 impacted=1
apex_time=1.186335e+01 apex_height=7.804948e+02
impact_time=2.512667e+01 range=3.628516e+03 downrange=3.546984e+03 crossrange=7.648704e+02
energy_lost=1.941071e+04