Monte Carlo Ballistic Simulator

Build a C++ ballistic simulator with RK4 integration and CEP analysis.

Introduction

30 Second Summary

Every time a projectile is fired in the real world, small variations in launch conditions mean it never lands in exactly the same spot twice. Defense engineers quantify this uncertainty using statistical simulation to predict how accurate a weapon system really is.

In this project, you will build a Monte Carlo ballistic trajectory simulator in C++. You will fire a projectile through a realistic atmosphere model with aerodynamic drag, perturb launch parameters across 1000 runs, and compute the defense industry's standard accuracy metric.

What You'll Build

A modular C++ simulation engine that combines RK4 numerical integration, ISA atmosphere modeling, and Monte Carlo statistical analysis to predict projectile dispersion.

By the end of this project, you'll have:

  • A C++ physics engine with a 4th-order Runge-Kutta integrator and altitude-dependent aerodynamic drag using the ISA atmosphere model.
  • A Monte Carlo execution layer with config-driven, seed-controlled reproducibility across 1000 simulations.
  • A Python visualization pipeline producing trajectory family plots and impact dispersion histograms with CEP annotation.
  • Secret Mission: Add crosswind to produce 3D impact scatter with a true 2D dispersion ellipse and CEP circle.

Are there any prerequisites?

This project assumes you have C++ fundamentals (variables, functions, structs, loops) and basic physics knowledge (forces, acceleration, Newton's second law). A conceptual understanding of the Monte Carlo method (running many randomized trials to estimate a statistical outcome) will also help.

Set Up Project Structure and Build System

Before writing any simulation logic, you need a solid foundation. Defense simulation codebases separate concerns into layers: physics engines, execution controllers, configuration, and analysis.

In this step, you will build SALVO's modular directory structure, its build system, and the header file that declares every function the simulator will use. Think of the header as the contract between modules. Once the project compiles with empty stubs, you know every piece can talk to every other piece.

In this step, get ready to:
  • Create the directory structure and Makefile for the SALVO project.
  • Write the project header file declaring all structs and function prototypes.
  • Write skeleton source files and verify the project compiles.
Create the directory structure and Makefile

SALVO separates its code into five directories, each with a clear responsibility: headers, source files, configuration, analysis scripts, and output data.

  • In VS Code, open the integrated terminal by pressing Ctrl+` (backtick key, top-left of your keyboard).
  • Move to your Desktop by running this command:
cd ~/Desktop
  • Create the SALVO project directory and all subdirectories by running these commands:
mkdir salvo
cd salvo
mkdir include src config analysis output
  • Confirm the directories exist by running:
ls

You should see analysis, config, include, output, and src listed in the output.

Don't see the directories?

Make sure you ran cd salvo before running ls. If the salvo folder was not created, check that you were on the Desktop first.

Help me troubleshoot my SALVO directory structure

  • Now open this project folder in VS Code. Click File in the top menu bar.
  • Click Open Folder.
  • Navigate to your Desktop, select the salvo folder, and click Select Folder.

VS Code's Explorer sidebar now shows your five subdirectories. Next, you will create the Makefile that tells the compiler how to build your project.

  • Create a new file in the root of your salvo folder called Makefile (no extension). Right-click in the Explorer sidebar and select New File, then type Makefile.
  • Paste the following build configuration into Makefile:
CXX = g++
CXXFLAGS = -std=c++17 -Wall -Wextra -O2 -Iinclude
SRC = src/main.cpp src/physics.cpp src/monte_carlo.cpp src/config_parser.cpp
TARGET = salvo.exe

all: $(TARGET)

$(TARGET): $(SRC) include/salvo.h
	$(CXX) $(CXXFLAGS) -o $(TARGET) $(SRC)

run: $(TARGET)
	@if not exist output mkdir output
	.\$(TARGET)

clean:
	@if exist $(TARGET) del $(TARGET)

.PHONY: all clean run

What does this Makefile do?

  • CXX = g++ sets your compiler to MinGW's g++.
  • CXXFLAGS enables C++17, turns on all warnings, applies -O2 optimization, and tells the compiler where to find headers.
  • all compiles the four source files into salvo.exe.
  • run creates the output/ directory if needed, then executes the simulator.
  • clean removes the compiled executable so you can rebuild fresh.

Why Windows-specific commands in the Makefile?

MinGW's mingw32-make uses cmd.exe as its shell, not Bash. That is why the Makefile uses if not exist and del instead of mkdir -p and rm.

  • Save Makefile.

Makefile not saving correctly?

Make sure the file is named exactly Makefile with a capital M and no file extension. If VS Code added .txt, rename it.

Help me fix my Makefile setup

Write the project header file

The header file is SALVO's central contract. It declares every data structure and function that the modules will implement. Writing it first forces you to design the full interface before coding any logic.

  • In the Explorer sidebar, right-click the include folder and select New File.
  • Name it salvo.h and paste the following code:
#ifndef SALVO_H
#define SALVO_H

#include <vector>
#include <string>

// Projectile state vector: position (x, y) and velocity (vx, vy)
struct State {
    double x;
    double y;
    double vx;
    double vy;
};

// Simulation configuration loaded from config file
struct Config {
    double muzzle_velocity;
    double launch_angle;
    double projectile_mass;
    double projectile_diameter;
    double drag_coefficient;

    int num_simulations;
    unsigned int seed;
    double dt;

    double sigma_velocity;
    double sigma_angle;
    double sigma_drag;
};

// Result of a single trajectory simulation
struct ImpactResult {
    double range;
    double flight_time;
};

// Physics engine functions
double air_density(double altitude);

State compute_derivatives(const State& state, double mass, double diameter, double cd);
State rk4_step(const State& state, double dt, double mass, double diameter, double cd);
ImpactResult simulate_trajectory(double velocity, double angle_deg, double drag_cd,
                                  const Config& config,
                                  std::vector<State>* trajectory);

// Monte Carlo execution
std::vector<ImpactResult> run_monte_carlo(const Config& config,
                                           std::vector<std::vector<State>>* trajectories);
double compute_cep(const std::vector<ImpactResult>& results);

// Config parser
Config parse_config(const std::string& filepath);

// Output writers
void write_results_csv(const std::string& filepath,
                       const std::vector<ImpactResult>& results);
void write_trajectories_csv(const std::string& filepath,
                            const std::vector<std::vector<State>>& trajectories);

#endif

What do these declarations define?

  • State holds a projectile's position and velocity at any moment in flight. This is your state vector for the integrator.
  • Config stores every parameter the simulator needs: projectile properties, simulation settings, and the standard deviations (sigmas) that control how much randomness each Monte Carlo run adds.
  • ImpactResult captures where and when a projectile hits the ground.
  • The function prototypes define four groups: physics (atmosphere, derivatives, RK4, trajectory), Monte Carlo execution, config parsing, and CSV output.
  • Save include/salvo.h.

Seeing red squiggles in the header?

VS Code's IntelliSense may show errors because no .cpp files exist yet. These will disappear once you create the source files in the next substep.

Help me understand the errors in salvo.h

Write skeleton source files and verify compilation

Each source file needs stub function bodies so the linker can resolve every prototype declared in salvo.h. The stubs return placeholder values. You will replace them with real logic in later steps.

  • In the Explorer sidebar, right-click the src folder and select New File. Name it main.cpp and paste:
#include "salvo.h"
#include <iostream>

// Entry point - will be expanded to load config and run simulation
int main() {
    std::cout << "SALVO - build successful" << std::endl;
    return 0;
}
  • Create another file in src called physics.cpp and paste:
#include "salvo.h"

// Stub: returns sea-level density as placeholder
double air_density(double altitude) {
    (void)altitude;
    return 1.225;
}

// Stub: returns zero derivatives
State compute_derivatives(const State& state, double mass, double diameter, double cd) {
    (void)state; (void)mass; (void)diameter; (void)cd;
    return {};
}

// Stub: returns state unchanged
State rk4_step(const State& state, double dt, double mass, double diameter, double cd) {
    (void)dt; (void)mass; (void)diameter; (void)cd;
    return state;
}

// Stub: returns zero impact result
ImpactResult simulate_trajectory(double velocity, double angle_deg, double drag_cd,
                                  const Config& config,
                                  std::vector<State>* trajectory) {
    (void)velocity; (void)angle_deg; (void)drag_cd;
    (void)config; (void)trajectory;
    return {};
}
  • Create monte_carlo.cpp in src and paste:
#include "salvo.h"

// Stub: returns empty results vector
std::vector<ImpactResult> run_monte_carlo(const Config& config,
                                           std::vector<std::vector<State>>* trajectories) {
    (void)config; (void)trajectories;
    return {};
}

// Stub: returns zero CEP
double compute_cep(const std::vector<ImpactResult>& results) {
    (void)results;
    return 0.0;
}

// Stub: no-op CSV writer
void write_results_csv(const std::string& filepath,
                       const std::vector<ImpactResult>& results) {
    (void)filepath; (void)results;
}

// Stub: no-op trajectory writer
void write_trajectories_csv(const std::string& filepath,
                            const std::vector<std::vector<State>>& trajectories) {
    (void)filepath; (void)trajectories;
}
  • Create config_parser.cpp in src and paste:
#include "salvo.h"

// Stub: returns default config values
Config parse_config(const std::string& filepath) {
    (void)filepath;
    Config config = {};
    return config;
}

Why use (void) casts in the stubs?

The Makefile enables -Wall -Wextra which warns about unused parameters. Casting a parameter to (void) tells the compiler you intentionally ignored it. This keeps the build warning-free even before you write the real logic.

  • Save all four files.
  • In the VS Code terminal, make sure you are inside the salvo folder (you should already be there from earlier).
  • Compile the project by running:
mingw32-make

You should see the compiler command echoed and then a clean return to the prompt with no errors. The file salvo.exe now exists in your project root.

  • Confirm the executable was created by running:
ls salvo.exe

You should see salvo.exe listed.

  • Run the stub executable to verify it works:
.\salvo.exe

You should see SALVO - build successful printed to the terminal.

Seeing 'mingw32-make: command not found'?

This means MinGW's bin directory is not on your system PATH. Confirm the path to mingw32-make.exe and add it to your PATH environment variable.

If you see a linker error about missing functions, double-check that all four .cpp files are saved in the src/ folder and that each includes "salvo.h" at the top.

Help me fix my mingw32-make build errors

✔️ Awesome, I've got everything!

Great. Double check you have saved all files: Makefile, include/salvo.h, and all four .cpp files in src/.

ⓧ I'd like to double check the full code

Here is every file you should have at this point.

CXX = g++
CXXFLAGS = -std=c++17 -Wall -Wextra -O2 -Iinclude
SRC = src/main.cpp src/physics.cpp src/monte_carlo.cpp src/config_parser.cpp
TARGET = salvo.exe

all: $(TARGET)

$(TARGET): $(SRC) include/salvo.h
	$(CXX) $(CXXFLAGS) -o $(TARGET) $(SRC)

run: $(TARGET)
	@if not exist output mkdir output
	.\$(TARGET)

clean:
	@if exist $(TARGET) del $(TARGET)

.PHONY: all clean run
#ifndef SALVO_H
#define SALVO_H

#include <vector>
#include <string>

struct State {
    double x;
    double y;
    double vx;
    double vy;
};

struct Config {
    double muzzle_velocity;
    double launch_angle;
    double projectile_mass;
    double projectile_diameter;
    double drag_coefficient;

    int num_simulations;
    unsigned int seed;
    double dt;

    double sigma_velocity;
    double sigma_angle;
    double sigma_drag;
};

struct ImpactResult {
    double range;
    double flight_time;
};

double air_density(double altitude);

State compute_derivatives(const State& state, double mass, double diameter, double cd);
State rk4_step(const State& state, double dt, double mass, double diameter, double cd);
ImpactResult simulate_trajectory(double velocity, double angle_deg, double drag_cd,
                                  const Config& config,
                                  std::vector<State>* trajectory);

std::vector<ImpactResult> run_monte_carlo(const Config& config,
                                           std::vector<std::vector<State>>* trajectories);
double compute_cep(const std::vector<ImpactResult>& results);

Config parse_config(const std::string& filepath);

void write_results_csv(const std::string& filepath,
                       const std::vector<ImpactResult>& results);
void write_trajectories_csv(const std::string& filepath,
                            const std::vector<std::vector<State>>& trajectories);

#endif
#include "salvo.h"
#include <iostream>

int main() {
    std::cout << "SALVO - build successful" << std::endl;
    return 0;
}
#include "salvo.h"

double air_density(double altitude) {
    (void)altitude;
    return 1.225;
}

State compute_derivatives(const State& state, double mass, double diameter, double cd) {
    (void)state; (void)mass; (void)diameter; (void)cd;
    return {};
}

State rk4_step(const State& state, double dt, double mass, double diameter, double cd) {
    (void)dt; (void)mass; (void)diameter; (void)cd;
    return state;
}

ImpactResult simulate_trajectory(double velocity, double angle_deg, double drag_cd,
                                  const Config& config,
                                  std::vector<State>* trajectory) {
    (void)velocity; (void)angle_deg; (void)drag_cd;
    (void)config; (void)trajectory;
    return {};
}
#include "salvo.h"

std::vector<ImpactResult> run_monte_carlo(const Config& config,
                                           std::vector<std::vector<State>>* trajectories) {
    (void)config; (void)trajectories;
    return {};
}

double compute_cep(const std::vector<ImpactResult>& results) {
    (void)results;
    return 0.0;
}

void write_results_csv(const std::string& filepath,
                       const std::vector<ImpactResult>& results) {
    (void)filepath; (void)results;
}

void write_trajectories_csv(const std::string& filepath,
                            const std::vector<std::vector<State>>& trajectories) {
    (void)filepath; (void)trajectories;
}
#include "salvo.h"

Config parse_config(const std::string& filepath) {
    (void)filepath;
    Config config = {};
    return config;
}

Your project compiles cleanly with zero logic. Every module knows its interface. Next up, you will fill in the physics engine with a real atmosphere model and RK4 integrator that simulates a projectile tearing through altitude-dependent drag.

Build the Physics Engine with RK4 and ISA Atmosphere

Your project skeleton compiles cleanly, but every function is empty. Time to build the brain of the simulator: the physics engine.

A ballistic projectile in flight experiences two forces: gravity pulling it down, and aerodynamic drag slowing it along its velocity vector. The drag force depends on air density, which changes with altitude. To capture this correctly, you need an atmosphere model, a force calculator, and a numerical integrator that steps the projectile through time.

In this step, get ready to:
  • Implement the ISA troposphere model and quadratic drag computation.
  • Build a 4th-order Runge-Kutta integrator from scratch.
  • Wire up a trajectory simulation loop and verify it with a test shot.
Implement the ISA atmosphere model and drag computation

Air density at sea level is about 1.225 kg/m³, but it drops significantly at higher altitudes. A projectile flying through a 10 km arc encounters much less drag at the top than at launch. The ISA (International Standard Atmosphere) troposphere model captures this relationship with a power-law formula.

  • Open src/physics.cpp in VS Code.
  • Replace the stub comment with the following code to set up the file preamble and atmosphere model:
#include "salvo.h"
#include <cmath>

// Physical constants
static const double G = 9.80665;
static const double PI = 3.14159265358979323846;

double air_density(double altitude) {
    // ISA troposphere constants
    const double rho_0 = 1.225;      // sea-level density (kg/m^3)
    const double T_0 = 288.15;       // sea-level temperature (K)
    const double L = 0.0065;         // temperature lapse rate (K/m)
    const double exponent = 4.2561;  // derived from g*M/(R*L) - 1

    // Clamp altitude to troposphere bounds (0 - 11 km)
    if (altitude < 0.0) altitude = 0.0;
    if (altitude > 11000.0) altitude = 11000.0;

    double temp_ratio = 1.0 - (L * altitude) / T_0;
    return rho_0 * std::pow(temp_ratio, exponent);
}

What does this code do?

The ISA model says air density in the troposphere follows a power law: ρ = ρ₀ × (1 - L×h/T₀)^4.2561. The exponent 4.2561 comes from the physical constants g, molar mass of air, the gas constant, and the lapse rate.

Clamping altitude to 0-11 km keeps the model valid. Above 11 km the lapse rate changes (the stratosphere), but artillery projectiles rarely reach that high.

You won't see any output yet. The next chunk adds the force computation that uses this density function.

  • Below air_density(), add the drag force computation:
State compute_derivatives(const State& state, double mass, double diameter, double cd) {
    // Total speed of the projectile
    double speed = std::sqrt(state.vx * state.vx + state.vy * state.vy);

    State deriv;
    // Position derivatives are just velocities
    deriv.x = state.vx;
    deriv.y = state.vy;

    if (speed < 1e-10) {
        // Projectile essentially stopped, only gravity acts
        deriv.vx = 0.0;
        deriv.vy = -G;
    } else {
        // Quadratic aerodynamic drag: F = 0.5 * rho * Cd * A * v^2
        double area = PI * (diameter / 2.0) * (diameter / 2.0);
        double rho = air_density(state.y);
        double drag_factor = 0.5 * rho * cd * area * speed / mass;

        // Drag opposes velocity in each axis
        deriv.vx = -drag_factor * state.vx;
        deriv.vy = -G - drag_factor * state.vy;
    }

    return deriv;
}

What does this code do?

This function computes the rate of change of the state vector [x, y, vx, vy]. Position derivatives are simply the current velocities.

Velocity derivatives come from two forces. Gravity pulls down at 9.80665 m/s². Drag acts opposite to the velocity vector with magnitude proportional to v² (quadratic drag). The drag_factor term is the standard formulation: 0.5 × ρ × Cd × A × |v| / m.

  • Save src/physics.cpp.

You won't see output yet. The next substep builds the integrator that calls compute_derivatives() repeatedly to step the projectile forward in time.

Seeing a compile error?

If you see an error about std::sqrt or std::pow, make sure #include <cmath> is at the top of the file.

If you see errors about State or Config being undeclared, check that #include "salvo.h" is the very first line.

help me fix a compile error in physics.cpp

Build the RK4 numerical integrator

A simple Euler integrator just steps forward once per timestep. That works for straight lines, but a curved trajectory accumulates large errors. RK4 samples the derivatives four times per step (start, two midpoints, end) and weights them to capture curvature far more accurately.

  • Below compute_derivatives(), add the RK4 integrator function:
State rk4_step(const State& state, double dt, double mass, double diameter, double cd) {
    // Sample 1: derivatives at current state
    State k1 = compute_derivatives(state, mass, diameter, cd);

    // Sample 2: derivatives at midpoint using k1
    State s2;
    s2.x = state.x + 0.5 * dt * k1.x;
    s2.y = state.y + 0.5 * dt * k1.y;
    s2.vx = state.vx + 0.5 * dt * k1.vx;
    s2.vy = state.vy + 0.5 * dt * k1.vy;
    State k2 = compute_derivatives(s2, mass, diameter, cd);

    // Sample 3: derivatives at midpoint using k2
    State s3;
    s3.x = state.x + 0.5 * dt * k2.x;
    s3.y = state.y + 0.5 * dt * k2.y;
    s3.vx = state.vx + 0.5 * dt * k2.vx;
    s3.vy = state.vy + 0.5 * dt * k2.vy;
    State k3 = compute_derivatives(s3, mass, diameter, cd);

    // Sample 4: derivatives at full step using k3
    State s4;
    s4.x = state.x + dt * k3.x;
    s4.y = state.y + dt * k3.y;
    s4.vx = state.vx + dt * k3.vx;
    s4.vy = state.vy + dt * k3.vy;
    State k4 = compute_derivatives(s4, mass, diameter, cd);

    // Weighted combination: (k1 + 2*k2 + 2*k3 + k4) / 6
    State next;
    next.x = state.x + (dt / 6.0) * (k1.x + 2.0 * k2.x + 2.0 * k3.x + k4.x);
    next.y = state.y + (dt / 6.0) * (k1.y + 2.0 * k2.y + 2.0 * k3.y + k4.y);
    next.vx = state.vx + (dt / 6.0) * (k1.vx + 2.0 * k2.vx + 2.0 * k3.vx + k4.vx);
    next.vy = state.vy + (dt / 6.0) * (k1.vy + 2.0 * k2.vy + 2.0 * k3.vy + k4.vy);

    return next;
}

Why four samples instead of one?

Each k-sample evaluates the derivative at a slightly different point along the timestep. k1 measures the slope at the start. k2 and k3 measure it at the midpoint (using different estimates). k4 measures it at the end.

The weighted average (1/6)(k1 + 2k2 + 2k3 + k4) cancels out lower-order error terms. The result is 4th-order accuracy: doubling the timestep only introduces 1/16th the error of Euler.

  • Save src/physics.cpp.

You won't see output yet. The next chunk builds the simulation loop that calls rk4_step() repeatedly until the projectile hits the ground.

Wire up the trajectory simulation loop and verify with a test shot

The simulation loop initializes the projectile at the origin with a given velocity and angle, then repeatedly calls the RK4 integrator until the projectile returns to ground level.

  • Below rk4_step(), add the trajectory simulation function:
ImpactResult simulate_trajectory(double velocity, double angle_deg, double drag_cd,
                                  const Config& config,
                                  std::vector<State>* trajectory) {
    // Convert launch angle from degrees to radians
    double angle_rad = angle_deg * PI / 180.0;

    // Initialize state at origin with given velocity components
    State state;
    state.x = 0.0;
    state.y = 0.0;
    state.vx = velocity * std::cos(angle_rad);
    state.vy = velocity * std::sin(angle_rad);

    // Record initial state if trajectory storage provided
    if (trajectory) {
        trajectory->push_back(state);
    }

    double t = 0.0;

    // Step until projectile returns to ground (y < 0)
    while (state.y >= 0.0 || t < config.dt * 2.0) {
        state = rk4_step(state, config.dt, config.projectile_mass,
                         config.projectile_diameter, drag_cd);
        t += config.dt;

        if (trajectory) {
            trajectory->push_back(state);
        }

        // Safety cutoff for runaway simulations
        if (t > 1000.0) break;
    }

    ImpactResult result;
    result.range = state.x;
    result.flight_time = t;
    return result;
}

What does this code do?

The function decomposes the muzzle velocity into horizontal (vx) and vertical (vy) components using trigonometry. Then it runs the RK4 loop with a fixed timestep (config.dt) until the projectile's y-coordinate goes negative (ground impact).

The t < config.dt * 2.0 condition prevents the loop from exiting on the very first step (when y starts at 0). The 1000-second safety cutoff prevents infinite loops if something goes wrong with the parameters.

  • Save src/physics.cpp.

✔️ Awesome, I've got everything!

Great. Double check you've saved src/physics.cpp with all four functions in order: air_density, compute_derivatives, rk4_step, simulate_trajectory.

ⓧ I'd like to double check the full code

#include "salvo.h"
#include <cmath>

// Physical constants
static const double G = 9.80665;
static const double PI = 3.14159265358979323846;

double air_density(double altitude) {
    // ISA troposphere constants
    const double rho_0 = 1.225;      // sea-level density (kg/m^3)
    const double T_0 = 288.15;       // sea-level temperature (K)
    const double L = 0.0065;         // temperature lapse rate (K/m)
    const double exponent = 4.2561;  // derived from g*M/(R*L) - 1

    // Clamp altitude to troposphere bounds (0 - 11 km)
    if (altitude < 0.0) altitude = 0.0;
    if (altitude > 11000.0) altitude = 11000.0;

    double temp_ratio = 1.0 - (L * altitude) / T_0;
    return rho_0 * std::pow(temp_ratio, exponent);
}

State compute_derivatives(const State& state, double mass, double diameter, double cd) {
    // Total speed of the projectile
    double speed = std::sqrt(state.vx * state.vx + state.vy * state.vy);

    State deriv;
    // Position derivatives are just velocities
    deriv.x = state.vx;
    deriv.y = state.vy;

    if (speed < 1e-10) {
        // Projectile essentially stopped, only gravity acts
        deriv.vx = 0.0;
        deriv.vy = -G;
    } else {
        // Quadratic aerodynamic drag: F = 0.5 * rho * Cd * A * v^2
        double area = PI * (diameter / 2.0) * (diameter / 2.0);
        double rho = air_density(state.y);
        double drag_factor = 0.5 * rho * cd * area * speed / mass;

        // Drag opposes velocity in each axis
        deriv.vx = -drag_factor * state.vx;
        deriv.vy = -G - drag_factor * state.vy;
    }

    return deriv;
}

State rk4_step(const State& state, double dt, double mass, double diameter, double cd) {
    // Sample 1: derivatives at current state
    State k1 = compute_derivatives(state, mass, diameter, cd);

    // Sample 2: derivatives at midpoint using k1
    State s2;
    s2.x = state.x + 0.5 * dt * k1.x;
    s2.y = state.y + 0.5 * dt * k1.y;
    s2.vx = state.vx + 0.5 * dt * k1.vx;
    s2.vy = state.vy + 0.5 * dt * k1.vy;
    State k2 = compute_derivatives(s2, mass, diameter, cd);

    // Sample 3: derivatives at midpoint using k2
    State s3;
    s3.x = state.x + 0.5 * dt * k2.x;
    s3.y = state.y + 0.5 * dt * k2.y;
    s3.vx = state.vx + 0.5 * dt * k2.vx;
    s3.vy = state.vy + 0.5 * dt * k2.vy;
    State k3 = compute_derivatives(s3, mass, diameter, cd);

    // Sample 4: derivatives at full step using k3
    State s4;
    s4.x = state.x + dt * k3.x;
    s4.y = state.y + dt * k3.y;
    s4.vx = state.vx + dt * k3.vx;
    s4.vy = state.vy + dt * k3.vy;
    State k4 = compute_derivatives(s4, mass, diameter, cd);

    // Weighted combination: (k1 + 2*k2 + 2*k3 + k4) / 6
    State next;
    next.x = state.x + (dt / 6.0) * (k1.x + 2.0 * k2.x + 2.0 * k3.x + k4.x);
    next.y = state.y + (dt / 6.0) * (k1.y + 2.0 * k2.y + 2.0 * k3.y + k4.y);
    next.vx = state.vx + (dt / 6.0) * (k1.vx + 2.0 * k2.vx + 2.0 * k3.vx + k4.vx);
    next.vy = state.vy + (dt / 6.0) * (k1.vy + 2.0 * k2.vy + 2.0 * k3.vy + k4.vy);

    return next;
}

ImpactResult simulate_trajectory(double velocity, double angle_deg, double drag_cd,
                                  const Config& config,
                                  std::vector<State>* trajectory) {
    // Convert launch angle from degrees to radians
    double angle_rad = angle_deg * PI / 180.0;

    // Initialize state at origin with given velocity components
    State state;
    state.x = 0.0;
    state.y = 0.0;
    state.vx = velocity * std::cos(angle_rad);
    state.vy = velocity * std::sin(angle_rad);

    // Record initial state if trajectory storage provided
    if (trajectory) {
        trajectory->push_back(state);
    }

    double t = 0.0;

    // Step until projectile returns to ground (y < 0)
    while (state.y >= 0.0 || t < config.dt * 2.0) {
        state = rk4_step(state, config.dt, config.projectile_mass,
                         config.projectile_diameter, drag_cd);
        t += config.dt;

        if (trajectory) {
            trajectory->push_back(state);
        }

        // Safety cutoff for runaway simulations
        if (t > 1000.0) break;
    }

    ImpactResult result;
    result.range = state.x;
    result.flight_time = t;
    return result;
}

Now you need a way to actually fire the engine and see a result. You will add a temporary test in main.cpp that fires a single projectile and prints its range.

  • Open src/main.cpp and replace the stub with this temporary test code:
#include "salvo.h"
#include <iostream>

int main() {
    // Temporary test: fire a single projectile at 45 degrees, 900 m/s
    Config config;
    config.projectile_mass = 10.0;
    config.projectile_diameter = 0.155;
    config.dt = 0.01;

    // Simulate with muzzle velocity 900 m/s, 45 deg, drag coefficient 0.3
    ImpactResult result = simulate_trajectory(900.0, 45.0, 0.3, config, nullptr);

    std::cout << "=== Physics Engine Test ===" << std::endl;
    std::cout << "Muzzle velocity: 900 m/s" << std::endl;
    std::cout << "Launch angle: 45 deg" << std::endl;
    std::cout << "Range: " << result.range << " m ("
              << result.range / 1000.0 << " km)" << std::endl;
    std::cout << "Flight time: " << result.flight_time << " s" << std::endl;
    std::cout << "Vacuum range (no drag): ~82.5 km" << std::endl;

    return 0;
}

Why compare against the vacuum range?

In a vacuum (no drag), a 900 m/s projectile at 45 degrees has an analytical range of about 82.5 km (v²/g). With realistic aerodynamic drag on a 155mm projectile, the range should be dramatically shorter, roughly 20-30 km.

If your range is close to 82.5 km, drag is not being applied. If it is negative or near zero, the integrator has a sign error. A value in the 20-30 km range confirms the physics engine is working correctly.

  • Save src/main.cpp.
  • Compile and run the project by running this command:
mingw32-make clean && mingw32-make && .\salvo.exe

You should see output showing a range somewhere between 20,000 and 30,000 meters (20-30 km), along with a flight time of roughly 50-80 seconds.

Range showing ~82.5 km or higher?

This means drag is not being applied. Check that air_density() is returning a non-zero value by verifying the exponent calculation. Also confirm that compute_derivatives() is using the speed variable (not a hardcoded zero).

If you see a negative range or zero, check that your angle conversion uses PI / 180.0 (degrees to radians). A common mistake is passing radians to std::cos when the angle is already in degrees.

help me debug my physics engine output

Extra for Experts: Why not just use Euler integration?

Euler integration evaluates the derivative once and steps blindly forward. On a curved trajectory, this accumulates significant error because it ignores how the derivative changes mid-step.

RK4 evaluates the derivative four times per step and weights the midpoint estimates more heavily. For the same timestep, RK4 is orders of magnitude more accurate. This is why it remains the standard integrator in defense ballistic simulations.

Your physics engine is alive. A single projectile flies through a realistic atmosphere and lands at a physically correct distance. Next, you will build the Monte Carlo layer that fires this engine hundreds of times with random perturbations to compute statistical accuracy.

Build the Monte Carlo Execution Layer and CEP Calculation

Your physics engine can simulate a single projectile from launch to impact. But a single trajectory tells you nothing about how accurate the weapon system actually is.

Real-world launches have uncertainty in muzzle velocity, launch angle, and aerodynamic properties. A Monte Carlo simulation quantifies this uncertainty by running hundreds of perturbed trajectories and measuring the resulting dispersion. In this step, you will build the config parser, stochastic execution layer, and CEP calculation that turns SALVO into a statistical analysis tool.

In this step, get ready to:
  • Implement the config parser and write the default configuration file.
  • Build the Monte Carlo execution loop with seeded Gaussian perturbations.
  • Add CEP computation, CSV output, and wire everything together in main.
Implement the config parser and write the config file

Defense simulations are config-driven. You never hardcode parameters into source code because reproducibility requires exact control over every input.

  • Create the configuration file by opening config/salvo.cfg in VS Code and pasting this content:
# SALVO Configuration File
# Monte Carlo Ballistic Trajectory Simulator

# Projectile parameters
muzzle_velocity = 900.0
launch_angle = 45.0
projectile_mass = 10.0
projectile_diameter = 0.155
drag_coefficient = 0.3

# Simulation parameters
num_simulations = 1000
seed = 42
dt = 0.01

# Perturbation standard deviations (Gaussian noise)
sigma_velocity = 5.0
sigma_angle = 0.5
sigma_drag = 0.01

What do these parameters control?

The first group defines the nominal (ideal) launch conditions for a 155mm-class projectile.

The sigma values define how much each parameter varies between shots. A sigma_velocity of 5.0 means each shot's muzzle velocity is drawn from a Gaussian distribution centered on 900 m/s with a standard deviation of 5 m/s. The seed value ensures every run with the same seed produces identical results.

  • Save config/salvo.cfg.
  • Confirm the file exists by running:
ls config/

You should see salvo.cfg in the output.

Now you need a parser that reads this file and populates the Config struct. The parser uses hardcoded defaults so the simulation runs even if the file is missing.

  • Replace the contents of src/config_parser.cpp with this implementation:
#include "salvo.h"
#include <fstream>
#include <sstream>
#include <iostream>

Config parse_config(const std::string& filepath) {
    Config config;
    // Hardcoded defaults - simulation runs even without a config file
    config.muzzle_velocity = 900.0;
    config.launch_angle = 45.0;
    config.projectile_mass = 10.0;
    config.projectile_diameter = 0.155;
    config.drag_coefficient = 0.3;
    config.num_simulations = 1000;
    config.seed = 42;
    config.dt = 0.01;
    config.sigma_velocity = 5.0;
    config.sigma_angle = 0.5;
    config.sigma_drag = 0.01;

    std::ifstream file(filepath);
    if (!file.is_open()) {
        std::cerr << "Warning: Could not open " << filepath
                  << ", using defaults." << std::endl;
        return config;
    }

    std::string line;
    while (std::getline(file, line)) {
        // Skip empty lines and comments
        if (line.empty() || line[0] == '#') continue;

        // Parse "key = value" format
        std::istringstream iss(line);
        std::string key, equals, value;
        if (iss >> key >> equals >> value) {
            if (key == "muzzle_velocity") config.muzzle_velocity = std::stod(value);
            else if (key == "launch_angle") config.launch_angle = std::stod(value);
            else if (key == "projectile_mass") config.projectile_mass = std::stod(value);
            else if (key == "projectile_diameter") config.projectile_diameter = std::stod(value);
            else if (key == "drag_coefficient") config.drag_coefficient = std::stod(value);
            else if (key == "num_simulations") config.num_simulations = std::stoi(value);
            else if (key == "seed") config.seed = static_cast<unsigned int>(std::stoul(value));
            else if (key == "dt") config.dt = std::stod(value);
            else if (key == "sigma_velocity") config.sigma_velocity = std::stod(value);
            else if (key == "sigma_angle") config.sigma_angle = std::stod(value);
            else if (key == "sigma_drag") config.sigma_drag = std::stod(value);
        }
    }

    return config;
}

What does this code do?

The parser first sets all fields to safe defaults. Then it attempts to open the config file. If the file is missing, it warns and returns the defaults.

For each non-empty, non-comment line, it splits on whitespace to extract the key and value. The std::stod and std::stoi calls convert strings to the appropriate numeric types.

  • Save src/config_parser.cpp.
  • Verify the project still compiles by running:
mingw32-make

You should see salvo.exe built without errors.

Seeing a compile error?

Check that you have the #include <fstream>, #include <sstream>, and #include <iostream> headers at the top. Missing any of these causes undefined symbol errors for the stream classes.

Also verify that every else if branch ends with a semicolon after the function call.

Help me debug my config_parser.cpp compile error.

✔️ Awesome, I've got everything!

Great. Double check you have saved both config/salvo.cfg and src/config_parser.cpp.

ⓧ I'd like to double check the full code

Here is the complete config/salvo.cfg:

# SALVO Configuration File
# Monte Carlo Ballistic Trajectory Simulator

# Projectile parameters
muzzle_velocity = 900.0
launch_angle = 45.0
projectile_mass = 10.0
projectile_diameter = 0.155
drag_coefficient = 0.3

# Simulation parameters
num_simulations = 1000
seed = 42
dt = 0.01

# Perturbation standard deviations (Gaussian noise)
sigma_velocity = 5.0
sigma_angle = 0.5
sigma_drag = 0.01

Here is the complete src/config_parser.cpp:

#include "salvo.h"
#include <fstream>
#include <sstream>
#include <iostream>

Config parse_config(const std::string& filepath) {
    Config config;
    config.muzzle_velocity = 900.0;
    config.launch_angle = 45.0;
    config.projectile_mass = 10.0;
    config.projectile_diameter = 0.155;
    config.drag_coefficient = 0.3;
    config.num_simulations = 1000;
    config.seed = 42;
    config.dt = 0.01;
    config.sigma_velocity = 5.0;
    config.sigma_angle = 0.5;
    config.sigma_drag = 0.01;

    std::ifstream file(filepath);
    if (!file.is_open()) {
        std::cerr << "Warning: Could not open " << filepath
                  << ", using defaults." << std::endl;
        return config;
    }

    std::string line;
    while (std::getline(file, line)) {
        if (line.empty() || line[0] == '#') continue;

        std::istringstream iss(line);
        std::string key, equals, value;
        if (iss >> key >> equals >> value) {
            if (key == "muzzle_velocity") config.muzzle_velocity = std::stod(value);
            else if (key == "launch_angle") config.launch_angle = std::stod(value);
            else if (key == "projectile_mass") config.projectile_mass = std::stod(value);
            else if (key == "projectile_diameter") config.projectile_diameter = std::stod(value);
            else if (key == "drag_coefficient") config.drag_coefficient = std::stod(value);
            else if (key == "num_simulations") config.num_simulations = std::stoi(value);
            else if (key == "seed") config.seed = static_cast<unsigned int>(std::stoul(value));
            else if (key == "dt") config.dt = std::stod(value);
            else if (key == "sigma_velocity") config.sigma_velocity = std::stod(value);
            else if (key == "sigma_angle") config.sigma_angle = std::stod(value);
            else if (key == "sigma_drag") config.sigma_drag = std::stod(value);
        }
    }

    return config;
}
Implement the Monte Carlo execution loop

The Monte Carlo layer runs the physics engine hundreds of times, each time with slightly different parameters drawn from Gaussian distributions. This models the real-world variability that makes every shot land in a slightly different place.

  • Replace the contents of src/monte_carlo.cpp with the includes and run_monte_carlo() function:
#include "salvo.h"
#include <random>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <cstdio>

std::vector<ImpactResult> run_monte_carlo(const Config& config,
                                           std::vector<std::vector<State>>* trajectories) {
    // Seed the Mersenne Twister PRNG for deterministic results
    std::mt19937 rng(config.seed);
    std::normal_distribution<double> vel_dist(config.muzzle_velocity, config.sigma_velocity);
    std::normal_distribution<double> angle_dist(config.launch_angle, config.sigma_angle);
    std::normal_distribution<double> drag_dist(config.drag_coefficient, config.sigma_drag);

    std::vector<ImpactResult> results;
    results.reserve(config.num_simulations);

    // Only store first 50 trajectories (for visualization, not all 1000)
    int trajectory_samples = std::min(config.num_simulations, 50);

    for (int i = 0; i < config.num_simulations; i++) {
        // Draw perturbed parameters from Gaussian distributions
        double vel = vel_dist(rng);
        double angle = angle_dist(rng);
        double cd = std::max(0.01, drag_dist(rng));

        std::vector<State> traj;
        std::vector<State>* traj_ptr = nullptr;
        if (trajectories && i < trajectory_samples) {
            traj_ptr = &traj;
        }

        ImpactResult result = simulate_trajectory(vel, angle, cd, config, traj_ptr);
        results.push_back(result);

        if (trajectories && i < trajectory_samples) {
            trajectories->push_back(std::move(traj));
        }

        // Print progress every 100 simulations
        if ((i + 1) % 100 == 0) {
            std::cout << "\rSimulations: " << (i + 1)
                      << "/" << config.num_simulations << std::flush;
        }
    }

    return results;
}

What does this code do?

The function seeds std::mt19937 (a Mersenne Twister pseudo-random number generator) with the config seed. This guarantees identical output for the same seed value.

Each iteration draws velocity, angle, and drag coefficient from independent Gaussian distributions centered on the nominal values. The std::max(0.01, ...) on drag prevents physically impossible negative drag coefficients.

Only the first 50 trajectories are recorded point-by-point (for plotting later). All 1000 impact results are stored for statistical analysis.

Why seed control matters in defense simulation

Production simulation tools require deterministic reproducibility. If an analyst reports a concerning CEP value, their manager must be able to rerun the exact same simulation and get the exact same number. The seed makes this possible.

  • Save src/monte_carlo.cpp.
  • You won't see a change yet. The next chunk adds the remaining functions to this file.
Add CEP calculation, CSV output, and wire main.cpp

CEP (Circular Error Probable) is the defense industry's standard accuracy metric. It is the radius within which 50% of projectiles land. For range-only analysis (no crosswind), CEP reduces to the median absolute deviation from the mean impact range.

  • In src/monte_carlo.cpp, add the following three functions below run_monte_carlo():
double compute_cep(const std::vector<ImpactResult>& results) {
    // Calculate mean impact range
    double mean_range = 0.0;
    for (const auto& r : results) {
        mean_range += r.range;
    }
    mean_range /= results.size();

    // Compute absolute deviation from mean for each impact
    std::vector<double> distances;
    distances.reserve(results.size());
    for (const auto& r : results) {
        double d = std::abs(r.range - mean_range);
        distances.push_back(d);
    }

    // CEP = median absolute deviation (50th percentile)
    std::sort(distances.begin(), distances.end());

    size_t idx = distances.size() / 2;
    return distances[idx];
}

void write_results_csv(const std::string& filepath,
                       const std::vector<ImpactResult>& results) {
    FILE* f = std::fopen(filepath.c_str(), "w");
    if (!f) {
        std::cerr << "Error: Cannot write to " << filepath << std::endl;
        return;
    }

    std::fprintf(f, "range,flight_time\n");
    for (const auto& r : results) {
        std::fprintf(f, "%.4f,%.4f\n", r.range, r.flight_time);
    }
    std::fclose(f);
}

void write_trajectories_csv(const std::string& filepath,
                            const std::vector<std::vector<State>>& trajectories) {
    FILE* f = std::fopen(filepath.c_str(), "w");
    if (!f) {
        std::cerr << "Error: Cannot write to " << filepath << std::endl;
        return;
    }

    std::fprintf(f, "trajectory_id,x,y,vx,vy\n");
    for (size_t i = 0; i < trajectories.size(); i++) {
        for (const auto& s : trajectories[i]) {
            std::fprintf(f, "%zu,%.4f,%.4f,%.4f,%.4f\n", i, s.x, s.y, s.vx, s.vy);
        }
    }
    std::fclose(f);
}

What does this code do?

The compute_cep() function computes the mean impact range, then finds the absolute deviation of each shot from that mean. Sorting and taking the median (50th percentile) gives the CEP.

The two CSV writers use C-style fprintf for fast formatted output. The impacts file stores range and flight time for all 1000 runs. The trajectories file stores the full position and velocity history for the 50 sample trajectories.

  • Save src/monte_carlo.cpp.
  • Verify the file compiles by running:
mingw32-make

You should see the build succeed with no errors.

Undefined reference errors?

Make sure all three functions (compute_cep, write_results_csv, write_trajectories_csv) are placed BELOW the closing brace of run_monte_carlo() in the same file. They must not be nested inside another function.

Help me fix linker errors in my monte_carlo.cpp file.

✔️ Awesome, I've got everything!

Great. Make sure you have saved src/monte_carlo.cpp with all four functions.

ⓧ I'd like to double check the full code

Here is the complete src/monte_carlo.cpp:

#include "salvo.h"
#include <random>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <cstdio>

std::vector<ImpactResult> run_monte_carlo(const Config& config,
                                           std::vector<std::vector<State>>* trajectories) {
    std::mt19937 rng(config.seed);
    std::normal_distribution<double> vel_dist(config.muzzle_velocity, config.sigma_velocity);
    std::normal_distribution<double> angle_dist(config.launch_angle, config.sigma_angle);
    std::normal_distribution<double> drag_dist(config.drag_coefficient, config.sigma_drag);

    std::vector<ImpactResult> results;
    results.reserve(config.num_simulations);

    int trajectory_samples = std::min(config.num_simulations, 50);

    for (int i = 0; i < config.num_simulations; i++) {
        double vel = vel_dist(rng);
        double angle = angle_dist(rng);
        double cd = std::max(0.01, drag_dist(rng));

        std::vector<State> traj;
        std::vector<State>* traj_ptr = nullptr;
        if (trajectories && i < trajectory_samples) {
            traj_ptr = &traj;
        }

        ImpactResult result = simulate_trajectory(vel, angle, cd, config, traj_ptr);
        results.push_back(result);

        if (trajectories && i < trajectory_samples) {
            trajectories->push_back(std::move(traj));
        }

        if ((i + 1) % 100 == 0) {
            std::cout << "\rSimulations: " << (i + 1)
                      << "/" << config.num_simulations << std::flush;
        }
    }

    return results;
}

double compute_cep(const std::vector<ImpactResult>& results) {
    double mean_range = 0.0;
    for (const auto& r : results) {
        mean_range += r.range;
    }
    mean_range /= results.size();

    std::vector<double> distances;
    distances.reserve(results.size());
    for (const auto& r : results) {
        double d = std::abs(r.range - mean_range);
        distances.push_back(d);
    }

    std::sort(distances.begin(), distances.end());

    size_t idx = distances.size() / 2;
    return distances[idx];
}

void write_results_csv(const std::string& filepath,
                       const std::vector<ImpactResult>& results) {
    FILE* f = std::fopen(filepath.c_str(), "w");
    if (!f) {
        std::cerr << "Error: Cannot write to " << filepath << std::endl;
        return;
    }

    std::fprintf(f, "range,flight_time\n");
    for (const auto& r : results) {
        std::fprintf(f, "%.4f,%.4f\n", r.range, r.flight_time);
    }
    std::fclose(f);
}

void write_trajectories_csv(const std::string& filepath,
                            const std::vector<std::vector<State>>& trajectories) {
    FILE* f = std::fopen(filepath.c_str(), "w");
    if (!f) {
        std::cerr << "Error: Cannot write to " << filepath << std::endl;
        return;
    }

    std::fprintf(f, "trajectory_id,x,y,vx,vy\n");
    for (size_t i = 0; i < trajectories.size(); i++) {
        for (const auto& s : trajectories[i]) {
            std::fprintf(f, "%zu,%.4f,%.4f,%.4f,%.4f\n", i, s.x, s.y, s.vx, s.vy);
        }
    }
    std::fclose(f);
}

Now wire everything together in main.cpp. The new main function loads the config, runs the Monte Carlo batch, computes statistics, and writes output files.

  • Replace the entire contents of src/main.cpp with this implementation:
#include "salvo.h"
#include <iostream>
#include <string>

int main(int argc, char* argv[]) {
    // Allow config path override via command line argument
    std::string config_path = "config/salvo.cfg";
    if (argc > 1) {
        config_path = argv[1];
    }

    std::cout << "=== SALVO - Monte Carlo Ballistic Trajectory Simulator ===" << std::endl;
    std::cout << "Loading config: " << config_path << std::endl;

    Config config = parse_config(config_path);

    // Print loaded configuration for verification
    std::cout << "Muzzle velocity: " << config.muzzle_velocity << " m/s" << std::endl;
    std::cout << "Launch angle: " << config.launch_angle << " deg" << std::endl;
    std::cout << "Projectile mass: " << config.projectile_mass << " kg" << std::endl;
    std::cout << "Projectile diameter: " << config.projectile_diameter << " m" << std::endl;
    std::cout << "Drag coefficient: " << config.drag_coefficient << std::endl;
    std::cout << "Simulations: " << config.num_simulations << std::endl;
    std::cout << "Seed: " << config.seed << std::endl;
    std::cout << std::endl;

    // Run the Monte Carlo batch
    std::vector<std::vector<State>> trajectories;
    std::vector<ImpactResult> results = run_monte_carlo(config, &trajectories);

    // Compute statistics
    double cep = compute_cep(results);

    double mean_range = 0.0;
    for (const auto& r : results) {
        mean_range += r.range;
    }
    mean_range /= results.size();

    // Print results
    std::cout << std::endl;
    std::cout << "=== Results ===" << std::endl;
    std::cout << "Mean range: " << mean_range << " m" << std::endl;
    std::cout << "CEP (Circular Error Probable): " << cep << " m" << std::endl;

    // Write CSV output for visualization
    write_results_csv("output/impacts.csv", results);
    write_trajectories_csv("output/trajectories.csv", trajectories);

    std::cout << "Output written to output/" << std::endl;
    std::cout << "Run: python analysis/visualize.py" << std::endl;

    return 0;
}

What does this code do?

The main function orchestrates the full pipeline. It loads config from config/salvo.cfg (or a custom path passed as a command-line argument), prints the loaded parameters, runs the Monte Carlo batch, computes mean range and CEP, and writes two CSV files to the output/ directory.

The progress counter inside run_monte_carlo() uses \r (carriage return) to overwrite the same line as it counts up. This keeps the terminal clean during long runs.

  • Save src/main.cpp.
  • Build and run the full simulation by running:
mingw32-make run

You should see the config parameters printed, a progress counter reaching 1000/1000, and then the results with a mean range and CEP value.

  • Run it a second time with the same command:
mingw32-make run

The output should be identical to the first run. This confirms deterministic seed control is working.

  • Verify the CSV files were generated by running:
ls output/

You should see impacts.csv and trajectories.csv in the listing.

Seeing a linker error about undefined references?

This means one of the functions called in main.cpp is not defined in the corresponding .cpp file. Verify that src/monte_carlo.cpp defines all four functions: run_monte_carlo, compute_cep, write_results_csv, and write_trajectories_csv.

Also check that src/config_parser.cpp defines parse_config. The function signatures must match include/salvo.h exactly.

Help me fix undefined reference errors in my SALVO build.

✔️ Awesome, I've got everything!

Great. Confirm that running salvo.exe twice produces the same mean range and CEP values.

ⓧ I'd like to double check the full code

Here is the complete src/main.cpp:

#include "salvo.h"
#include <iostream>
#include <string>

int main(int argc, char* argv[]) {
    std::string config_path = "config/salvo.cfg";
    if (argc > 1) {
        config_path = argv[1];
    }

    std::cout << "=== SALVO - Monte Carlo Ballistic Trajectory Simulator ===" << std::endl;
    std::cout << "Loading config: " << config_path << std::endl;

    Config config = parse_config(config_path);

    std::cout << "Muzzle velocity: " << config.muzzle_velocity << " m/s" << std::endl;
    std::cout << "Launch angle: " << config.launch_angle << " deg" << std::endl;
    std::cout << "Projectile mass: " << config.projectile_mass << " kg" << std::endl;
    std::cout << "Projectile diameter: " << config.projectile_diameter << " m" << std::endl;
    std::cout << "Drag coefficient: " << config.drag_coefficient << std::endl;
    std::cout << "Simulations: " << config.num_simulations << std::endl;
    std::cout << "Seed: " << config.seed << std::endl;
    std::cout << std::endl;

    std::vector<std::vector<State>> trajectories;
    std::vector<ImpactResult> results = run_monte_carlo(config, &trajectories);

    double cep = compute_cep(results);

    double mean_range = 0.0;
    for (const auto& r : results) {
        mean_range += r.range;
    }
    mean_range /= results.size();

    std::cout << std::endl;
    std::cout << "=== Results ===" << std::endl;
    std::cout << "Mean range: " << mean_range << " m" << std::endl;
    std::cout << "CEP (Circular Error Probable): " << cep << " m" << std::endl;

    write_results_csv("output/impacts.csv", results);
    write_trajectories_csv("output/trajectories.csv", trajectories);

    std::cout << "Output written to output/" << std::endl;
    std::cout << "Run: python analysis/visualize.py" << std::endl;

    return 0;
}

Your Monte Carlo pipeline is complete. You now have 1000 simulated impacts with statistical dispersion analysis. Next up, you will build the Python visualization layer that turns those CSV files into publication-quality trajectory and histogram plots.

Build the Python Visualization Layer

Your C++ simulator is generating thousands of data points across 1000 Monte Carlo runs. Right now those results live in CSV files that are hard to interpret at a glance.

In this step, you will build a Python visualization script that reads the CSV output and produces publication-quality plots. One plot overlays all 50 sample trajectories to show the family of flight paths. The other shows a histogram of impact ranges with statistical markers for mean and CEP.

In this step, get ready to:
  • Write the Python visualization script with data loaders and matplotlib plotting.
  • Run the full simulation-to-visualization pipeline and verify the output figure.
Write the CSV data loaders

The visualization script needs to parse two CSV files: output/impacts.csv (1000 rows of impact range and flight time) and output/trajectories.csv (position data for 50 sample trajectories).

  • Create the file analysis/visualize.py and add the imports and both data-loading functions by pasting this code:
import matplotlib.pyplot as plt
import numpy as np
import csv
import os

def load_impacts(filepath):
    # Read impact range and flight time from each simulation run
    ranges = []
    flight_times = []
    with open(filepath, 'r') as f:
        reader = csv.DictReader(f)
        for row in reader:
            ranges.append(float(row['range']))
            flight_times.append(float(row['flight_time']))
    return np.array(ranges), np.array(flight_times)

def load_trajectories(filepath):
    # Group trajectory points by ID into a dictionary of x/y lists
    trajectories = {}
    with open(filepath, 'r') as f:
        reader = csv.DictReader(f)
        for row in reader:
            tid = int(row['trajectory_id'])
            if tid not in trajectories:
                trajectories[tid] = {'x': [], 'y': []}
            trajectories[tid]['x'].append(float(row['x']))
            trajectories[tid]['y'].append(float(row['y']))
    return trajectories

What does this code do?

These two functions handle all the CSV parsing for the visualization layer.

  • load_impacts() reads each row from the impacts file and returns NumPy arrays of ranges and flight times. NumPy arrays enable vectorized math (like computing means and medians across all 1000 values at once).
  • load_trajectories() groups position data by trajectory ID into a dictionary. Each key maps to a dict with x and y lists, giving you one plottable series per trajectory.

You won't see a change yet. The next chunk adds the plotting logic that uses these loaders.

Write the plotting main function

The main() function computes summary statistics, creates a 1x2 figure with a trajectory family plot and an impact histogram, then saves the result as a PNG.

  • Add the following main() function and entry-point guard below load_trajectories() in analysis/visualize.py:
def main():
    # Load simulation output data
    ranges, flight_times = load_impacts('output/impacts.csv')
    trajectories = load_trajectories('output/trajectories.csv')

    # Compute statistics: mean impact range and CEP (median absolute deviation)
    mean_range = np.mean(ranges)
    cep = np.median(np.abs(ranges - mean_range))

    # Create a 1x2 figure for trajectory family and impact histogram
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    fig.suptitle('SALVO - Monte Carlo Ballistic Trajectory Analysis',
                 fontsize=14, fontweight='bold')

    # Left plot: overlay all sample trajectories (range vs altitude)
    ax1 = axes[0]
    for tid, traj in trajectories.items():
        x_km = np.array(traj['x']) / 1000.0
        y_km = np.array(traj['y']) / 1000.0
        ax1.plot(x_km, y_km, alpha=0.3, linewidth=0.8, color='steelblue')
    ax1.set_xlabel('Range (km)')
    ax1.set_ylabel('Altitude (km)')
    ax1.set_title('Trajectory Family (N={})'.format(len(trajectories)))
    ax1.set_ylim(bottom=0)
    ax1.grid(True, alpha=0.3)

    # Right plot: impact range histogram with mean and CEP markers
    ax2 = axes[1]
    ax2.hist(ranges / 1000.0, bins=40, color='steelblue', edgecolor='navy', alpha=0.7)
    ax2.axvline(mean_range / 1000.0, color='red', linestyle='--', linewidth=2,
                label='Mean: {:.2f} km'.format(mean_range / 1000.0))
    ax2.axvline((mean_range - cep) / 1000.0, color='orange', linestyle=':', linewidth=1.5)
    ax2.axvline((mean_range + cep) / 1000.0, color='orange', linestyle=':', linewidth=1.5,
                label='CEP: {:.1f} m'.format(cep))
    ax2.set_xlabel('Impact Range (km)')
    ax2.set_ylabel('Count')
    ax2.set_title('Impact Dispersion (N={})'.format(len(ranges)))
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    # Save the figure to disk
    plt.tight_layout()
    os.makedirs('output', exist_ok=True)
    plt.savefig('output/salvo_analysis.png', dpi=150, bbox_inches='tight')
    plt.show()
    print('Analysis saved to output/salvo_analysis.png')

if __name__ == '__main__':
    main()

What does this code do?

  • The CEP is computed as the median absolute deviation from the mean range. This mirrors the C++ computation and gives you a Python-side annotation value.
  • The left subplot loops over all 50 trajectory dictionaries and plots each as a semi-transparent line. Overlaying them reveals the dispersion envelope.
  • The right subplot draws a 40-bin histogram of impact ranges. A red dashed line marks the mean. Two orange dotted lines mark one CEP distance above and below the mean.
  • plt.savefig() writes the figure at 150 DPI before plt.show() displays it interactively. This ensures the file is always saved even if you close the window early.
  • Save analysis/visualize.py.

Seeing a syntax error when you save?

  • Check that the main() function is at the top level (no accidental indentation under load_trajectories).
  • Verify you have a blank line between the end of load_trajectories() and the start of def main():.
  • help me fix a syntax error in my visualize.py

✔️ Awesome, I've got everything!

Great. Make sure you have saved the file before moving on.

ⓧ I'd like to double check the full code

Here is the complete analysis/visualize.py for reference:

import matplotlib.pyplot as plt
import numpy as np
import csv
import os

def load_impacts(filepath):
    ranges = []
    flight_times = []
    with open(filepath, 'r') as f:
        reader = csv.DictReader(f)
        for row in reader:
            ranges.append(float(row['range']))
            flight_times.append(float(row['flight_time']))
    return np.array(ranges), np.array(flight_times)

def load_trajectories(filepath):
    trajectories = {}
    with open(filepath, 'r') as f:
        reader = csv.DictReader(f)
        for row in reader:
            tid = int(row['trajectory_id'])
            if tid not in trajectories:
                trajectories[tid] = {'x': [], 'y': []}
            trajectories[tid]['x'].append(float(row['x']))
            trajectories[tid]['y'].append(float(row['y']))
    return trajectories

def main():
    ranges, flight_times = load_impacts('output/impacts.csv')
    trajectories = load_trajectories('output/trajectories.csv')

    mean_range = np.mean(ranges)
    cep = np.median(np.abs(ranges - mean_range))

    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    fig.suptitle('SALVO - Monte Carlo Ballistic Trajectory Analysis',
                 fontsize=14, fontweight='bold')

    ax1 = axes[0]
    for tid, traj in trajectories.items():
        x_km = np.array(traj['x']) / 1000.0
        y_km = np.array(traj['y']) / 1000.0
        ax1.plot(x_km, y_km, alpha=0.3, linewidth=0.8, color='steelblue')
    ax1.set_xlabel('Range (km)')
    ax1.set_ylabel('Altitude (km)')
    ax1.set_title('Trajectory Family (N={})'.format(len(trajectories)))
    ax1.set_ylim(bottom=0)
    ax1.grid(True, alpha=0.3)

    ax2 = axes[1]
    ax2.hist(ranges / 1000.0, bins=40, color='steelblue', edgecolor='navy', alpha=0.7)
    ax2.axvline(mean_range / 1000.0, color='red', linestyle='--', linewidth=2,
                label='Mean: {:.2f} km'.format(mean_range / 1000.0))
    ax2.axvline((mean_range - cep) / 1000.0, color='orange', linestyle=':', linewidth=1.5)
    ax2.axvline((mean_range + cep) / 1000.0, color='orange', linestyle=':', linewidth=1.5,
                label='CEP: {:.1f} m'.format(cep))
    ax2.set_xlabel('Impact Range (km)')
    ax2.set_ylabel('Count')
    ax2.set_title('Impact Dispersion (N={})'.format(len(ranges)))
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    os.makedirs('output', exist_ok=True)
    plt.savefig('output/salvo_analysis.png', dpi=150, bbox_inches='tight')
    plt.show()
    print('Analysis saved to output/salvo_analysis.png')

if __name__ == '__main__':
    main()
Run the full simulation-to-visualization pipeline

Now you will run the C++ simulator to regenerate fresh CSV data, then feed that data into the Python visualization script.

  • Rebuild and run the simulator by running this command:
mingw32-make run
  • You should see the simulation progress counter reach 1000/1000, followed by the mean range and CEP values printed to the console.
  • Run the visualization script by running this command:
python analysis/visualize.py
  • A matplotlib window opens showing two side-by-side plots. Close the window to continue.
  • Verify that output/salvo_analysis.png was generated by running:
ls output/salvo_analysis.png

You should see the file listed. The left plot shows 50 parabolic arcs with visible spread between them. The right plot shows a roughly Gaussian histogram with a red dashed mean line and orange dotted CEP boundaries.

What should I see in the plots?

The trajectory family plot shows all 50 sample arcs overlaid. They fan out slightly because each run perturbed the muzzle velocity, launch angle, and drag coefficient.

The histogram shows where all 1000 projectiles landed on the range axis. The orange CEP markers bound the region that contains roughly 50% of impacts around the mean. This is the same CEP value the C++ code prints to the console.

Seeing a ModuleNotFoundError?

  • If the error mentions matplotlib or numpy, install them by running pip install matplotlib numpy and try again.
  • If you see a FileNotFoundError for the CSV files, make sure you ran mingw32-make run first. The script expects output/impacts.csv and output/trajectories.csv to already exist.
  • help me debug my SALVO visualization pipeline

Your simulation pipeline is complete. You have a C++ physics engine that fires projectiles through a realistic atmosphere, a Monte Carlo layer that quantifies dispersion, and a Python visualization that makes the results interpretable. Next up, the Secret Mission extends this to 3D with crosswind.

Secret mission

Add Crosswind and 2D Impact Scatter

Your simulator models range dispersion, but real projectiles drift sideways too. Crosswind pushes the projectile laterally, turning your 1D impact line into a 2D impact cloud. Expand SALVO to three dimensions, compute true circular CEP on the impact plane, and visualize the dispersion ellipse with a CEP circle overlay.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs or cloud resources to worry about. Decide whether to keep your project files as-is, remove the compiled executable, or delete everything.

Resources you used:

  • Local salvo project directory on your Desktop (source code, config, analysis scripts).
  • Compiled executable salvo.exe in the project root.
  • Generated output files in output/ (CSV data and PNG visualization).

Keep everything running

No action needed. Choose this if you want to keep the project as a portfolio piece or continue experimenting with different config parameters.

  • The project is entirely local. Your source code, config files, and output data stay on your machine.
  • You can re-run the simulation any time with mingw32-make run.
  • Tweak values in config/salvo.cfg to explore different scenarios.

Pause - I'll come back to this later

Remove the compiled executable to free up a small amount of disk space while keeping all your source code and outputs intact.

  • Navigate to your salvo project directory and run:
mingw32-make clean
  • This deletes salvo.exe only. Your source code, config, and output files remain.
  • To rebuild later, run mingw32-make from the project directory.

Delete - I don't want to use this again

Remove all project files from your machine entirely. This deletes your source code, compiled executable, config, and all output data.

  • Delete the entire project directory by running the command for your operating system.

On macOS or Linux:

rm -rf ~/Desktop/salvo

On Windows:

rmdir /s /q %USERPROFILE%\Desktop\salvo
  • This removes all source files, the compiled executable, config files, and generated output.
  • This action cannot be undone. Make sure you have copied anything you want to keep.

Nice Work!

Nice Work!

You just built SALVO, a Monte Carlo ballistic trajectory simulator from the ground up in C++. This is the kind of simulation infrastructure that defense and aerospace teams actually build to quantify weapon system accuracy.

You've learned how to:

  • Build a 3DOF point-mass ballistic physics engine with a 4th-order Runge-Kutta integrator and altitude-dependent aerodynamic drag using the ISA atmosphere model.
  • Design a Monte Carlo simulation layer with Gaussian perturbations, config-driven execution, and seed-controlled reproducibility for deterministic results.
  • Compute Circular Error Probable (CEP), the defense industry's standard accuracy metric, and produce trajectory family and impact dispersion visualizations with matplotlib.
  • Secret Mission: Extended the simulator to 3D with crosswind, producing true 2D impact scatter with a CEP circle visualization.

Ready to quiz yourself?