pressure-wave — ray-traced blast overpressure on body armor
pressure-wave is a Monte-Carlo-style simulator that estimates how much blast overpressure reaches a sensor at the center of a soldier’s torso under four armor configurations — full armor, helmet only, vest only, and bare — for a cube of candidate blast points around the body. It is a small research codebase aimed at quantifying the protective contribution of a MICH-style ballistic helmet and a SAPI plate vest against the impulse of a free-air blast, and at producing a tidy CSV that downstream R scripts can run an ANOVA and paired comparisons against.
Goals
- Estimate, for a given peak pressure, the total impulse delivered to a sensor at the head/chest origin from blast points scattered across a configurable cube.
- Compare four armor conditions on the same set of blast points so each point acts as its own paired “subject” for the statistical analysis.
- Stay simple enough to read end-to-end: a few hundred lines of
Python for the physics, a few dozen lines of R for the stats, and
two
.objmeshes for the geometry. - Be reproducible — every run records its parameters via the
CLI flags and writes a timestamped CSV next to the others in
data/.
How it works
The simulator (model.py) places
a spherical sensor at the origin and a blast point in space, then casts
a fan of rays from the blast point covering just the angular extents
needed to reach the sensor and any intervening geometry. Each ray is
followed by trace_ray: it advances to its next hit,
either the sensor sphere (via
line_sphere_intersect) or a triangle on the helmet or
vest mesh (via find_mesh_intersection, backed by trimesh
ray casting). On a geometry hit the ray reflects specularly with an
energy-loss factor of 0.8 per bounce and continues until it hits the
sensor, runs past --ray-max-length-meters, or exceeds
--ray-max-bounces.
When a ray finally strikes the sensor, the peak pressure at that
moment is taken to be the initial peak p0 attenuated
by an inverse-square distance law and the accumulated reflection
factor. The waveform is then modelled with a Friedlander pulse and
integrated over a 30 ms window to get an impulse contribution,
which is weighted by the differential solid angle of the ray and
summed into the total impulse for that armor configuration. The
driver (controller.py)
runs this four times per blast point — with both meshes, just
the helmet, just the vest, and no geometry — and the CLI
(simulate_blasts.py)
sweeps the blast point across a cube using
ProcessPoolExecutor for parallelism.
The R side (analysis.R) reads
a results CSV, splits blast points into near and intermediate-range
bands, pivots to long form with each blast point as its own subject,
and runs a repeated-measures ANOVA on impulse vs. armor condition.
When the omnibus test is significant, it follows up with paired
t-tests across the six condition pairs and applies a Holm
correction. It also reports the helmet’s effect size as both an
absolute impulse reduction (with 95% CI) and a Cohen’s
d on the per-subject percent reduction.
visualize.R contains the
box-, violin-, and 3D scatter plots used to eyeball the same data.
Inputs and outputs
- Geometry. Two
.objmeshes inassets/: a MICH-style helmet and a SAPI plate. They’re rotated and scaled into the project frame bygeometry.py— the helmet faces −Y, with Z up and X to the soldier’s left. - Parameters. Set per run on the CLI: the blast cube’s center, extent, and segmentation; the minimum stand-off distance; the initial peak pressure (default 50,000 kPa); ray resolution in arc-minutes; max bounces and max ray length; and the worker count.
- Output. A timestamped
blast-results-YYYY-MM-DD-HHMM.csvindata/with one row per blast point and columnsbp_x, bp_y, bp_z, full_armor, helmet_only, vest_only, no_armor. The repository ships several such files from prior runs at different resolutions and cube sizes.
License
AGPL-3.0-or-later · Copyright © 2025 SWGY, Inc. Source, assets, data, and the exploration notebook are all covered; see LICENSE for the full text.