Ballista — A Monte-Carlo Trajectory Simulation and Analysis

Goal

Most student ballistic sims stop at "plot a parabola." Ballista's goal was to go further in two specific directions:

1. Physical realism — model the things that actually determine where a real projectile lands: Mach-dependent drag (not a single flat drag coefficient), altitude-dependent air density, and wind that grows with altitude and couples correctly into the drag calculation rather than being bolted on as an afterthought.
2. Simulation engineering rigor — validate the numerical method against known closed-form answers instead of just trusting the output, and check that the Monte Carlo layer has actually converged (i.e. that running it again with a different random seed gives you the same answer) instead of reporting whatever one random draw happened to produce.

The end product answers a real ballistics question: given a launch setup with some amount of uncertainty in muzzle velocity, angle, drag, and wind, where does the projectile land, and how tightly clustered are the impacts? That clustering is expressed as CEP (Circular Error Probable) — the radius of a circle, centered on the mean impact point, that contains half of all impacts. It's the standard way accuracy is quantified in real exterior ballistics.

Screenshots

Trajectory family, range dispersion, and impact scatter

This is the core output of a single simulation run (2000 randomized launches from one nominal setup). Three views of the same data:

• Left — Trajectory Family. 50 sampled trajectories (out of the 2000 total impacts) plotted as altitude vs. range, colored by how far that particular shot ended up landing. This is what "dispersion" looks like in flight, before it's collapsed down to a single impact point — you can see the trajectories fan out mid-flight as their randomized velocity, angle, drag, and wind pull them apart.
• Middle — Range Dispersion. A histogram of where all 2000 shots landed along the range axis. The dashed red line is the mean impact range; the dotted orange lines mark one CEP radius on either side of the mean. In this run: mean range ≈ 7.54 km.
• Right — Impact Scatter with CEP Circle. The same 2000 impacts, but now plotted as (range error, lateral error) relative to the mean impact point — this is the actual "shot group" a real ballistics analyst would look at. The red circle has radius equal to the CEP (≈165 m here), meaning roughly half of the 2000 points fall inside it and half fall outside. A tighter cluster means a more precise weapon system; CEP is the single number that summarizes how tight that cluster is.

Input sensitivity analysis

This answers a different question than the plot above: not "how dispersed are the impacts," but "which of the four randomized inputs (velocity, angle, form factor, wind) is actually responsible for that dispersion?" Each Monte Carlo run already logs its randomly sampled inputs alongside its result, so this is a linear regression of range (and separately, lateral deviation) against each input, scaled by that input's own standard deviation — so a velocity change (hundreds of m/s) and a form-factor change (a small dimensionless number) can be compared on the same footing.

• Left — Range Sensitivity. Form factor dominates: a one-standard-deviation change in form factor moves range by about 119 m, roughly 4x the effect of a one-standard-deviation change in velocity. This makes physical sense — form factor directly scales drag, and range is very sensitive to drag over a multi-kilometer flight.
• Middle — Lateral Deviation Sensitivity. Wind dominates almost completely here, which also makes sense: velocity, angle, and form factor mostly affect how far the shot goes, not how far sideways it drifts. Crosswind is the one input that has a direct, sustained sideways effect over the whole flight time.
• Right — Strongest Driver. A direct scatter of range vs. form factor (the strongest range driver from the left panel), with a fitted trend line, so the tornado-chart bar is backed by the raw data rather than just a summary number.

Monte Carlo stability across 100 seeds

The two plots above are both from a single run with a single random seed. This one asks: if you'd used a different seed, would you have gotten a meaningfully different answer? analysis/seed_stability.py reruns the full 2000-impact simulation 100 times, once per seed, and looks at how much the headline numbers (mean range, CEP) move around.

• Left & middle — histograms of per-seed mean range and per-seed CEP. Both are tightly clustered: mean range varies by only 0.04% of its value across 100 independent runs, and CEP by about 1.5%. That's the confirmation that num_simulations = 2000 is a large enough sample for these numbers to be trustworthy, not an artifact of one lucky/unlucky draw.
• Right — CEP vs. seed value, a flat scatter with a dashed mean line. The point is the absence of a pattern: if CEP trended upward with seed value, or clustered at certain seed values, that would suggest a bug (e.g. a seeding or RNG-state issue). It doesn't — it's random noise around a stable mean, which is exactly what you want to see. One seed (22) landed at CEP ≈ 169.7 m, about 3.4 standard deviations out; rerunning that single seed reproduces the same number exactly, confirming it's genuine sampling variance (unsurprising at 100 draws of a median-based statistic) rather than nondeterminism.

Trajectory animation

A single sampled trajectory played back frame-by-frame rather than shown as a static arc — mostly to make the flight readable as motion (speed near launch and near apex, acceleration on the way down) rather than just a shape.

Engineering decisions

RK4 over Euler. Projectile motion under drag is stiff enough that Euler integration accumulates visible energy error over a multi-second flight. RK4 gives 4th-order convergence for a small extra cost per step — confirmed empirically below, not just assumed.

Mach-dependent drag instead of constant Cd. Real projectile drag isn't a single number; it rises sharply through the transonic region around Mach 1 and decays afterward. drag_coefficient_g1() implements a simplified approximation of the G1 standard drag curve, and the config exposes a form_factor (the ballistic term "i") that scales it — this is how real exterior ballistics separates shape (the standard curve) from this specific projectile (the form factor), rather than hand-tuning one flat Cd to match an expected range.

Relative-velocity drag coupling for crosswind. Wind isn't a separate additive term on top of drag — it changes the relative velocity the projectile experiences, which is what drag actually acts on. compute_derivatives() subtracts the local wind from velocity before computing drag, then applies the resulting force to the projectile's absolute state. A naive implementation that adds wind as a separate force gets this interaction wrong.

Altitude-dependent air density and wind shear. Air density falls off with altitude via the standard barometric formula, and crosswind speed grows with altitude via a power-law boundary-layer profile (the Hellmann exponent, wind_shear_exponent in the config). Both feed into the same drag calculation, so a trajectory that goes higher genuinely experiences different drag and different wind than one that stays low.

Impact-point interpolation. The integrator advances in fixed steps, so the raw last state after the projectile crosses y=0 has already overshot the ground. interpolate_impact() linearly interpolates between the last state above ground and the first state below it to recover the true crossing point. This directly affects the CEP number — without it, CEP is biased by however coarse dt happens to be.

Assert-based tests instead of a test framework. No Catch2/GoogleTest dependency, on purpose — for a project this size, a header-only assert macro keeps the build simple and portable (including across the MinGW toolchain) without sacrificing what the tests actually check.

CMake instead of a bare Makefile. Lets ballista_core be built once and linked into both the main binary and the test binary, and makes ctest a first-class way to run validation. On MinGW, the build also statically links the C++ runtime (-static -static-libgcc -static-libstdc++) so the resulting .exe doesn't depend on whichever runtime DLLs happen to be on PATH at launch time — this bit us during development when a conda environment's own bundled DLLs shadowed the ones the binary was actually linked against, and static linking removes that whole failure mode.

Validation

Four automated tests run against the physics core (tests/test_physics.cpp, run via ctest or the ballista_tests binary directly):

Test	What it checks	Result
Zero-drag range	Simulated range with drag disabled matches the closed-form v²sin(2θ)/g	Error: 0.00000 m
RK4 convergence order	Halving dt should shrink integration error by ~2⁴=16× for a 4th-order method	Ratio: 16.01
CEP against known circle	compute_cep() recovers the exact radius of a synthetic circular impact distribution	Error: 0.00000 m
Impact interpolation	Linear interpolation lands exactly on y=0, not past it	Error: < 1e-9 m

The RK4 test deliberately runs in a subsonic, wind-free regime, so the drag model's Mach-table kinks and the wind profile's altitude clamp don't locally break smoothness — those are properties of the simplified drag/wind models, not of the integrator, and isolating them keeps the convergence claim honest.

Monte Carlo convergence, swept across 100 independent RNG seeds (analysis/seed_stability.py, N=2000 impacts per seed):

Metric	Mean	Std dev	Std as % of mean
Mean range	7540.4 m	3.0 m	0.04%
CEP	161.5 m	2.4 m	1.48%

Run it yourself:

ctest --output-on-failure
python analysis/seed_stability.py --exe build/ballista --n 100

Project structure

ballista/
├── CMakeLists.txt
├── config/
│   └── ballista.cfg          # simulation parameters
├── include/
│   └── ballista.h
├── src/
│   ├── main.cpp
│   ├── physics.cpp           # RK4 integrator, drag model, wind model
│   ├── monte_carlo.cpp       # randomized sampling, CEP, CSV output
│   └── config_parser.cpp
├── tests/
│   └── test_physics.cpp      # 4 assert-based validation tests
├── analysis/
│   ├── visualize.py          # trajectory family / dispersion / CEP scatter
│   ├── sensitivity.py        # standardized sensitivity of range & lateral error
│   ├── seed_stability.py     # Monte Carlo convergence across seeds
│   └── animate.py            # single-trajectory GIF
└── docs/images/               # curated screenshots used in this README

Build & run

Linux / macOS

mkdir build && cd build
cmake .. -DCMAKE_BUILD_TYPE=Release
make -j4
ctest --output-on-failure
cd ..
./build/ballista config/ballista.cfg

Windows (MinGW)

mkdir build
cd build
cmake .. -G "MinGW Makefiles" -DCMAKE_BUILD_TYPE=Release
mingw32-make -j4
./ballista_tests.exe
cd ..
./build/ballista.exe config/ballista.cfg

Analysis (either platform)

pip install matplotlib numpy pillow
python analysis/visualize.py
python analysis/sensitivity.py
python analysis/animate.py 0
python analysis/seed_stability.py --exe build/ballista --n 100

Configuration reference

config/ballista.cfg controls the simulation:

Key	Meaning
muzzle_velocity	Nominal launch speed, m/s
launch_angle	Nominal launch angle, degrees
projectile_mass, projectile_diameter	Physical properties feeding the drag calculation
form_factor	Ballistic form factor (i) scaling the G1 drag curve — the main "how draggy is this specific projectile" knob
num_simulations	Monte Carlo sample count per run
seed	RNG seed — same seed always reproduces the same impacts
dt	Integration timestep, seconds
sigma_velocity, sigma_angle, sigma_form_factor	Per-run randomization std devs around the nominal values above
crosswind_speed, sigma_crosswind	Nominal crosswind at wind_ref_altitude, and its randomization
wind_ref_altitude	Reference altitude for the wind profile, m
wind_shear_exponent	Power-law exponent for wind growth with altitude

Unrecognized keys produce a warning and are ignored; malformed values keep their default and print which line failed, rather than crashing.

Possible extensions

• Sweep num_simulations itself (500 / 2000 / 8000) and plot CEP std-dev vs. N to show the √N Monte Carlo convergence rate directly, rather than only fixed-N seed stability.
• Compare CEP output against a published reference dispersion figure for a real weapon system, and discuss where the simplified drag/wind models over- or under-estimate.
• Replace the piecewise-linear G1 approximation with a full tabulated drag function for a smoother RHS (would also let the RK4 convergence test run in the supersonic regime without losing order).

License

MIT — see LICENSE.