use crate::cluster_bc::ClusterBCDegradation;
use crate::mc_stats::{
wilson_interval, BernoulliConfidenceSequence, ConfidenceLevel, Welford,
};
use crate::pitch_damping::{calculate_pitch_damping_coefficient, PitchDampingCoefficients};
use crate::precession_nutation::{
calculate_combined_angular_motion, projectile_moments_of_inertia, AngularState,
PrecessionNutationParams,
};
use crate::trajectory_sampling::{
projected_sample_count, sample_trajectory, TrajectoryData, TrajectoryOutputs,
TrajectorySample,
};
use crate::trajectory_observation::{bracket_param, Bracket, TrajectoryTermination};
use crate::wind_shear::WindShearModel;
use crate::DragModel;
use nalgebra::{Vector3, Vector6};
use std::error::Error;
use std::fmt;
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "cli", derive(clap::ValueEnum))]
pub enum UnitSystem {
Metric,
Imperial,
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum OutputFormat {
Table,
Json,
Csv,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
#[cfg_attr(feature = "cli", derive(clap::ValueEnum))]
pub enum BcReferenceStandard {
#[default]
Icao,
ArmyStandardMetro,
}
pub const BC_REFERENCE_STANDARD_INERT_WARNING: &str =
"warning: --bc-reference army-standard-metro has no effect together with a custom drag \
table (--drag-table): the deck's Cd is divided by sectional density, not a BC value, so \
no BC-reference conversion applies";
#[derive(Debug)]
pub struct BallisticsError {
message: String,
}
impl fmt::Display for BallisticsError {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
write!(f, "{}", self.message)
}
}
impl Error for BallisticsError {}
impl From<String> for BallisticsError {
fn from(msg: String) -> Self {
BallisticsError { message: msg }
}
}
impl From<&str> for BallisticsError {
fn from(msg: &str) -> Self {
BallisticsError {
message: msg.to_string(),
}
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum DropsReference {
#[default]
Los,
Target,
}
#[derive(Debug, Clone)]
pub struct BallisticInputs {
pub bc_value: f64, pub bc_type: DragModel, pub bc_reference_standard: BcReferenceStandard,
pub bullet_mass: f64, pub muzzle_velocity: f64, pub bullet_diameter: f64, pub bullet_length: f64,
pub muzzle_angle: f64, pub target_distance: f64, pub azimuth_angle: f64, pub shot_azimuth: f64,
pub shooting_angle: f64, pub cant_angle: f64,
pub sight_height: f64, pub sight_offset_lateral_m: f64,
pub muzzle_height: f64, pub target_height: f64, pub zero_poi_vertical_m: f64,
pub zero_poi_horizontal_m: f64,
pub ground_threshold: f64,
pub altitude: f64, pub temperature: f64, pub pressure: f64, pub humidity: f64,
pub latitude: Option<f64>,
pub wind_speed: f64, pub wind_angle: f64,
pub twist_rate: f64, pub is_twist_right: bool, pub caliber_inches: f64, pub weight_grains: f64, pub manufacturer: Option<String>, pub bullet_model: Option<String>, pub bullet_id: Option<String>, pub bullet_cluster: Option<usize>,
pub use_rk4: bool, pub use_adaptive_rk45: bool,
pub enable_advanced_effects: bool,
pub enable_magnus: bool, pub enable_coriolis: bool, pub use_powder_sensitivity: bool,
pub powder_temp_sensitivity: f64, pub powder_temp: f64, pub powder_temp_curve: Option<Vec<(f64, f64)>>,
pub powder_curve_temp_c: Option<f64>,
pub tipoff_yaw: f64, pub tipoff_decay_distance: f64, pub use_bc_segments: bool,
pub bc_segments: Option<Vec<(f64, f64)>>, pub bc_segments_data: Option<Vec<crate::BCSegmentData>>, pub use_enhanced_spin_drift: bool,
pub use_form_factor: bool,
pub enable_wind_shear: bool,
pub wind_shear_model: String,
pub enable_trajectory_sampling: bool,
pub sample_interval: f64, pub drops_reference: DropsReference,
pub enable_pitch_damping: bool,
pub enable_precession_nutation: bool,
pub enable_aerodynamic_jump: bool,
pub use_cluster_bc: bool,
pub custom_drag_table: Option<crate::drag::DragTable>,
pub cd_scale: f64,
pub bc_type_str: Option<String>,
}
impl BallisticInputs {
pub fn humidity_percent(&self) -> f64 {
(self.humidity * 100.0).clamp(0.0, 100.0)
}
pub fn windage_zero_bias_rad(&self, zero_distance_m: f64) -> f64 {
if zero_distance_m > 0.0 {
(self.zero_poi_horizontal_m + self.sight_offset_lateral_m) / zero_distance_m
} else {
0.0
}
}
pub fn sectional_density_lb_in2(&self) -> Option<f64> {
let weight_gr = if self.weight_grains > 0.0 {
self.weight_grains
} else {
self.bullet_mass / crate::constants::GRAINS_TO_KG };
let diameter_in = if self.caliber_inches > 0.0 {
self.caliber_inches
} else {
self.bullet_diameter / 0.0254 };
if weight_gr > 0.0 && diameter_in > 0.0 {
Some(weight_gr / 7000.0 / (diameter_in * diameter_in))
} else {
None
}
}
pub fn custom_drag_denominator(&self, fallback_bc: f64) -> f64 {
match self.sectional_density_lb_in2() {
Some(sd) => sd,
None => {
static WARN_ONCE: std::sync::Once = std::sync::Once::new();
WARN_ONCE.call_once(|| {
eprintln!(
"Warning: custom drag table active but bullet mass/diameter are \
unavailable; falling back to bc_value for the retardation denominator"
);
});
fallback_bc
}
}
}
pub fn bc_reference_standard_inert_warning(&self) -> Option<&'static str> {
if self.custom_drag_table.is_some()
&& matches!(self.bc_reference_standard, BcReferenceStandard::ArmyStandardMetro)
{
Some(BC_REFERENCE_STANDARD_INERT_WARNING)
} else {
None
}
}
pub fn normalize_for_solve(&mut self) {
if matches!(
self.bc_reference_standard,
BcReferenceStandard::ArmyStandardMetro
) {
self.bc_value *= crate::constants::ASM_TO_ICAO_BC;
if let Some(segments) = self.bc_segments.as_mut() {
for (_mach, bc) in segments.iter_mut() {
*bc *= crate::constants::ASM_TO_ICAO_BC;
}
}
if let Some(segments) = self.bc_segments_data.as_mut() {
for segment in segments.iter_mut() {
segment.bc_value *= crate::constants::ASM_TO_ICAO_BC;
}
}
self.bc_reference_standard = BcReferenceStandard::Icao;
}
self.caliber_inches = self.bullet_diameter / 0.0254;
self.weight_grains = self.bullet_mass / crate::constants::GRAINS_TO_KG;
self.muzzle_velocity = resolve_powder_adjusted_velocity(
self.muzzle_velocity,
self.temperature,
self.use_powder_sensitivity,
self.powder_temp_sensitivity,
self.powder_temp,
self.powder_temp_curve.as_deref(),
self.powder_curve_temp_c,
);
}
}
impl Default for BallisticInputs {
fn default() -> Self {
let mass_kg = 0.01;
let diameter_m = 0.00762;
let bc = 0.5;
let muzzle_angle_rad = 0.0;
let bc_type = DragModel::G1;
Self {
bc_value: bc,
bc_type,
bc_reference_standard: BcReferenceStandard::Icao,
bullet_mass: mass_kg,
muzzle_velocity: 800.0,
bullet_diameter: diameter_m,
bullet_length: crate::stability::estimate_bullet_length_m(diameter_m, mass_kg),
muzzle_angle: muzzle_angle_rad,
target_distance: 100.0,
azimuth_angle: 0.0,
shot_azimuth: 0.0,
shooting_angle: 0.0,
cant_angle: 0.0,
sight_height: 0.05,
sight_offset_lateral_m: 0.0, muzzle_height: 0.0, target_height: 0.0, zero_poi_vertical_m: 0.0, zero_poi_horizontal_m: 0.0,
ground_threshold: -100.0,
altitude: 0.0,
temperature: 15.0,
pressure: 1013.25, humidity: 0.5, latitude: None,
wind_speed: 0.0,
wind_angle: 0.0,
twist_rate: 12.0, is_twist_right: true,
caliber_inches: diameter_m / 0.0254, weight_grains: mass_kg / crate::constants::GRAINS_TO_KG, manufacturer: None,
bullet_model: None,
bullet_id: None,
bullet_cluster: None,
use_rk4: true, use_adaptive_rk45: true,
enable_advanced_effects: false,
enable_magnus: false,
enable_coriolis: false,
use_powder_sensitivity: false,
powder_temp_sensitivity: 0.0,
powder_temp: 15.0,
powder_temp_curve: None,
powder_curve_temp_c: None,
tipoff_yaw: 0.0,
tipoff_decay_distance: 50.0,
use_bc_segments: false,
bc_segments: None,
bc_segments_data: None,
use_enhanced_spin_drift: false,
use_form_factor: false,
enable_wind_shear: false,
wind_shear_model: "none".to_string(),
enable_trajectory_sampling: false,
sample_interval: 10.0, drops_reference: DropsReference::Los, enable_pitch_damping: false,
enable_precession_nutation: false,
enable_aerodynamic_jump: false,
use_cluster_bc: false,
custom_drag_table: None,
cd_scale: 1.0,
bc_type_str: None,
}
}
}
pub fn interpolate_powder_temp_curve(curve: &[(f64, f64)], temp_c: f64) -> f64 {
debug_assert!(!curve.is_empty());
if curve.is_empty() {
return 0.0;
}
let mut sorted;
let pts: &[(f64, f64)] = if curve.windows(2).all(|w| w[0].0 <= w[1].0) {
curve
} else {
sorted = curve.to_vec();
sorted.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
&sorted
};
let n = pts.len();
if temp_c <= pts[0].0 {
return pts[0].1; }
if temp_c >= pts[n - 1].0 {
return pts[n - 1].1; }
for i in 1..n {
let (t0, v0) = pts[i - 1];
let (t1, v1) = pts[i];
if temp_c <= t1 {
let span = t1 - t0;
if span.abs() < f64::EPSILON {
return v1; }
let f = (temp_c - t0) / span;
return v0 + f * (v1 - v0);
}
}
pts[n - 1].1
}
pub fn parse_powder_sweep(s: &str) -> Result<Vec<f64>, String> {
const MAX_SWEEP_ROWS: usize = 500;
let parts: Vec<&str> = s.split(':').collect();
if parts.len() != 3 {
return Err(format!(
"Invalid --sweep '{}': expected START:END:STEP (e.g. \"20:110:10\")",
s
));
}
let parse = |p: &str, name: &str| -> Result<f64, String> {
p.trim()
.parse::<f64>()
.map_err(|_| format!("Invalid --sweep {}: '{}' is not a number", name, p.trim()))
};
let start = parse(parts[0], "START")?;
let end = parse(parts[1], "END")?;
let step = parse(parts[2], "STEP")?;
if !step.is_finite() || step <= 0.0 {
return Err(format!("Invalid --sweep STEP {}: must be positive", step));
}
if !start.is_finite() || !end.is_finite() || end < start {
return Err(format!(
"Invalid --sweep range {}:{}: END must be >= START",
start, end
));
}
let n_f = ((end - start) / step + 1e-9).floor();
if !n_f.is_finite() || n_f + 1.0 > MAX_SWEEP_ROWS as f64 {
return Err(format!(
"--sweep would produce more than {} rows; use a larger STEP",
MAX_SWEEP_ROWS
));
}
let n = n_f as usize + 1;
Ok((0..n).map(|i| start + step * i as f64).collect())
}
pub fn resolve_powder_adjusted_velocity(
nominal_velocity_mps: f64,
ambient_temperature_c: f64,
use_powder_sensitivity: bool,
powder_temp_sensitivity_mps_per_c: f64,
powder_reference_temp_c: f64,
powder_temp_curve: Option<&[(f64, f64)]>,
powder_curve_temp_c: Option<f64>,
) -> f64 {
if let Some(curve) = powder_temp_curve {
if !curve.is_empty() {
let lookup_c = powder_curve_temp_c.unwrap_or(ambient_temperature_c);
return interpolate_powder_temp_curve(curve, lookup_c);
}
return nominal_velocity_mps;
}
if use_powder_sensitivity {
let temp_delta_c = ambient_temperature_c - powder_reference_temp_c;
return nominal_velocity_mps + powder_temp_sensitivity_mps_per_c * temp_delta_c;
}
nominal_velocity_mps
}
#[derive(Debug, Clone)]
pub struct WindConditions {
pub speed: f64, pub direction: f64,
pub vertical_speed: f64,
}
impl Default for WindConditions {
fn default() -> Self {
Self {
speed: 0.0,
direction: 0.0,
vertical_speed: 0.0,
}
}
}
#[derive(Debug, Clone)]
pub struct AtmosphericConditions {
pub temperature: f64, pub pressure: f64, pub humidity: f64,
pub altitude: f64, }
impl Default for AtmosphericConditions {
fn default() -> Self {
Self {
temperature: 15.0,
pressure: 1013.25,
humidity: 50.0,
altitude: 0.0,
}
}
}
#[derive(Debug, Clone)]
pub struct TrajectoryPoint {
pub time: f64,
pub position: Vector3<f64>,
pub velocity_magnitude: f64,
pub kinetic_energy: f64,
pub drag_coefficient: Option<f64>,
}
impl TrajectoryPoint {
pub fn drag_coefficient_json_value(&self, with_drag_coefficient: bool) -> Option<f64> {
if with_drag_coefficient {
self.drag_coefficient
} else {
None
}
}
}
#[derive(Debug, Clone)]
pub struct TrajectoryResult {
pub max_range: f64,
pub max_height: f64,
pub time_of_flight: f64,
pub impact_velocity: f64,
pub impact_energy: f64,
pub projectile_mass_kg: f64,
pub line_of_sight_height_m: f64,
pub station_speed_of_sound_mps: f64,
pub termination: TrajectoryTermination,
pub points: Vec<TrajectoryPoint>,
pub sampled_points: Option<Vec<TrajectorySample>>, pub min_pitch_damping: Option<f64>, pub transonic_mach: Option<f64>, pub angular_state: Option<AngularState>, pub max_yaw_angle: Option<f64>, pub max_precession_angle: Option<f64>, pub aerodynamic_jump: Option<crate::aerodynamic_jump::AerodynamicJumpComponents>,
pub mach_1_2_distance_m: Option<f64>,
pub mach_1_0_distance_m: Option<f64>,
pub mach_0_9_distance_m: Option<f64>,
}
const RK45_TOLERANCE: f64 = 1e-6;
const RK45_SAFETY_FACTOR: f64 = 0.9;
const RK45_MAX_DT: f64 = 0.01;
const RK45_MIN_DT: f64 = 1e-6;
const TRAJECTORY_TIME_LIMIT_S: f64 = 100.0;
pub const MAX_TRAJECTORY_POINTS: usize = 250_000;
fn cli_rk45_error_norm(
position: &Vector3<f64>,
velocity: &Vector3<f64>,
fifth_position: &Vector3<f64>,
fifth_velocity: &Vector3<f64>,
fourth_position: &Vector3<f64>,
fourth_velocity: &Vector3<f64>,
) -> f64 {
let pack_state = |position: &Vector3<f64>, velocity: &Vector3<f64>| {
Vector6::new(
position.x, position.y, position.z, velocity.x, velocity.y, velocity.z,
)
};
let state = pack_state(position, velocity);
let fifth_order = pack_state(fifth_position, fifth_velocity);
let fourth_order = pack_state(fourth_position, fourth_velocity);
crate::trajectory_integration::rk45_error_norm(&state, &fifth_order, &fourth_order)
}
struct Rk45Trial {
position: Vector3<f64>,
velocity: Vector3<f64>,
suggested_dt: f64,
error: f64,
}
struct Rk45AcceptedStep {
position: Vector3<f64>,
velocity: Vector3<f64>,
used_dt: f64,
next_dt: f64,
error: f64,
}
#[derive(Default)]
struct MachTransitionTracker {
previous_mach: Option<f64>,
crossed_transonic: bool,
crossed_subsonic: bool,
crossed_narrow: bool,
mach_1_2_distance_m: Option<f64>,
mach_1_0_distance_m: Option<f64>,
mach_0_9_distance_m: Option<f64>,
}
impl MachTransitionTracker {
fn record_downward_crossings(&mut self, mach: f64, downrange_m: f64, distances: &mut Vec<f64>) {
if !mach.is_finite() {
self.previous_mach = None;
return;
}
if let Some(previous_mach) = self.previous_mach {
if !self.crossed_transonic && previous_mach >= 1.2 && mach < 1.2 {
self.crossed_transonic = true;
distances.push(downrange_m);
self.mach_1_2_distance_m = Some(downrange_m);
}
if !self.crossed_subsonic && previous_mach >= 1.0 && mach < 1.0 {
self.crossed_subsonic = true;
distances.push(downrange_m);
self.mach_1_0_distance_m = Some(downrange_m);
}
if !self.crossed_narrow && previous_mach >= 0.9 && mach < 0.9 {
self.crossed_narrow = true;
self.mach_0_9_distance_m = Some(downrange_m);
}
}
self.previous_mach = Some(mach);
}
}
impl TrajectoryResult {
pub fn position_at_range(&self, target_range: f64) -> Option<Vector3<f64>> {
if self.points.is_empty() {
return None;
}
for i in 0..self.points.len() - 1 {
let p1 = &self.points[i];
let p2 = &self.points[i + 1];
if p1.position.x <= target_range && p2.position.x >= target_range {
let dx = p2.position.x - p1.position.x;
if dx.abs() < 1e-10 {
return Some(p1.position);
}
let t = (target_range - p1.position.x) / dx;
return Some(Vector3::new(
target_range,
p1.position.y + t * (p2.position.y - p1.position.y),
p1.position.z + t * (p2.position.z - p1.position.z),
));
}
}
self.points.last().map(|p| p.position)
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum StationAtmosphereResolution {
LegacyDefaultSentinels,
Authoritative,
}
#[derive(Clone)]
pub struct TrajectorySolver {
inputs: BallisticInputs,
wind: WindConditions,
atmosphere: AtmosphericConditions,
station_atmosphere_resolution: StationAtmosphereResolution,
max_range: f64,
time_step: f64,
max_trajectory_points: usize,
cluster_bc: Option<ClusterBCDegradation>,
precession_nutation_inertias: (f64, f64),
wind_sock: Option<crate::wind::WindSock>,
atmo_sock: Option<crate::atmosphere::AtmoSock>,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(crate) enum ZeroTargetFrame {
SightLine,
WorldVertical,
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct ZeroCrossings {
pub near_m: Option<f64>,
pub far_m: Option<f64>,
}
impl TrajectorySolver {
pub fn new(
inputs: BallisticInputs,
wind: WindConditions,
atmosphere: AtmosphericConditions,
) -> Self {
Self::new_with_station_atmosphere_resolution(
inputs,
wind,
atmosphere,
StationAtmosphereResolution::LegacyDefaultSentinels,
)
}
pub fn new_with_resolved_station_atmosphere(
inputs: BallisticInputs,
wind: WindConditions,
atmosphere: AtmosphericConditions,
) -> Self {
Self::new_with_station_atmosphere_resolution(
inputs,
wind,
atmosphere,
StationAtmosphereResolution::Authoritative,
)
}
fn new_with_station_atmosphere_resolution(
mut inputs: BallisticInputs,
wind: WindConditions,
atmosphere: AtmosphericConditions,
station_atmosphere_resolution: StationAtmosphereResolution,
) -> Self {
inputs.normalize_for_solve();
let cluster_bc = if inputs.use_cluster_bc {
Some(ClusterBCDegradation::new())
} else {
None
};
let precession_nutation_inertias = projectile_moments_of_inertia(
inputs.bullet_mass,
inputs.bullet_diameter,
inputs.bullet_length,
);
Self {
inputs,
wind,
atmosphere,
station_atmosphere_resolution,
max_range: 1000.0,
time_step: 0.001,
max_trajectory_points: MAX_TRAJECTORY_POINTS,
cluster_bc,
precession_nutation_inertias,
wind_sock: None,
atmo_sock: None,
}
}
pub fn set_max_range(&mut self, range: f64) {
self.max_range = range;
}
pub fn set_time_step(&mut self, step: f64) {
self.time_step = step;
}
pub(crate) fn calculate_and_set_zero_angle(
&mut self,
target_distance_m: f64,
target_height_m: f64,
frame: ZeroTargetFrame,
) -> Result<f64, BallisticsError> {
let angle = self.find_zero_angle(target_distance_m, target_height_m, frame)?;
let angle = if target_distance_m > 0.0 {
angle + self.inputs.zero_poi_vertical_m / target_distance_m
} else {
angle
};
self.inputs.muzzle_angle = angle;
self.apply_windage_zero_bias(target_distance_m);
Ok(angle)
}
pub(crate) fn apply_windage_zero_bias(&mut self, target_distance_m: f64) {
self.inputs.azimuth_angle += self.inputs.windage_zero_bias_rad(target_distance_m);
}
fn find_zero_angle(
&self,
target_distance_m: f64,
target_height_m: f64,
frame: ZeroTargetFrame,
) -> Result<f64, BallisticsError> {
let mut low_angle = 0.0;
let mut high_angle = 0.2; let tolerance = 1e-7;
let max_iterations = 60;
let low_height = self.zero_trial_height_at(low_angle, target_distance_m, frame)?;
let high_height = self.zero_trial_height_at(high_angle, target_distance_m, frame)?;
match (low_height, high_height) {
(Some(low_height), Some(high_height)) => {
let low_error = low_height - target_height_m;
let high_error = high_height - target_height_m;
if low_error > 0.0 && high_error > 0.0 {
} else if low_error < 0.0 && high_error < 0.0 {
let mut expanded = false;
for multiplier in [2.0, 3.0, 4.0] {
let new_high = (high_angle * multiplier).min(0.785);
if let Ok(Some(height)) =
self.zero_trial_height_at(new_high, target_distance_m, frame)
{
if height - target_height_m > 0.0 {
high_angle = new_high;
expanded = true;
break;
}
}
if new_high >= 0.785 {
break;
}
}
if !expanded {
return Err("Cannot find zero angle: target beyond effective range even at maximum angle".into());
}
}
}
(None, Some(_)) => {
}
(Some(_), None) => {
return Err(
"Cannot find zero angle: high angle trajectory doesn't reach target distance"
.into(),
);
}
(None, None) => {
return Err(
"Cannot find zero angle: trajectory cannot reach target distance at any angle"
.into(),
);
}
}
for _ in 0..max_iterations {
let mid_angle = (low_angle + high_angle) / 2.0;
match self.zero_trial_height_at(mid_angle, target_distance_m, frame)? {
Some(height) => {
let error = height - target_height_m;
if error.abs() < 0.0001 {
return Ok(mid_angle);
}
if (high_angle - low_angle).abs() < tolerance {
if error.abs() < 0.01 {
return Ok(mid_angle);
}
return Err("Zero angle did not converge: residual height error too large (target not reachable / not bracketed)".into());
}
if error > 0.0 {
high_angle = mid_angle;
} else {
low_angle = mid_angle;
}
}
None => {
low_angle = mid_angle;
if (high_angle - low_angle).abs() < tolerance {
return Err("Trajectory cannot reach target distance - angle converged without valid solution".into());
}
}
}
}
Err("Failed to find zero angle".into())
}
fn zero_trial_height_at(
&self,
angle_rad: f64,
target_distance_m: f64,
frame: ZeroTargetFrame,
) -> Result<Option<f64>, BallisticsError> {
let mut trial = self.clone();
trial.inputs.muzzle_angle = angle_rad;
trial.inputs.enable_aerodynamic_jump = false;
trial.inputs.cant_angle = 0.0;
if frame == ZeroTargetFrame::SightLine {
trial.inputs.shooting_angle = 0.0;
}
trial.set_max_range(target_distance_m * 2.0);
let result = trial.solve()?;
for (index, point) in result.points.iter().enumerate() {
if point.position.x >= target_distance_m {
let shot_y_m = if index == 0 {
point.position.y
} else {
let previous = &result.points[index - 1];
let span = point.position.x - previous.position.x;
let fraction = (target_distance_m - previous.position.x) / span;
previous.position.y + fraction * (point.position.y - previous.position.y)
};
return Ok(Some(crate::atmosphere::shot_frame_altitude(
0.0,
target_distance_m,
shot_y_m,
trial.inputs.shooting_angle,
)));
}
}
Ok(None)
}
fn find_zero_range(
&self,
angle_rad: f64,
target_height_m: f64,
frame: ZeroTargetFrame,
) -> Result<ZeroCrossings, BallisticsError> {
let mut trial = self.clone();
trial.inputs.muzzle_angle = angle_rad;
trial.inputs.enable_aerodynamic_jump = false;
trial.inputs.cant_angle = 0.0;
if frame == ZeroTargetFrame::SightLine {
trial.inputs.shooting_angle = 0.0;
}
let result = trial.solve()?;
let mut near_crossing: Option<f64> = None;
let mut far_crossing: Option<f64> = None;
let mut previous: Option<(f64, f64)> = None; for point in &result.points {
let height = crate::atmosphere::shot_frame_altitude(
0.0,
point.position.x,
point.position.y,
trial.inputs.shooting_angle,
);
let error = height - target_height_m;
if let Some((prev_x, prev_error)) = previous {
if prev_error == 0.0 {
if near_crossing.is_none() {
near_crossing = Some(prev_x);
} else {
far_crossing = Some(prev_x);
}
}
if prev_error * error < 0.0 {
let fraction = prev_error / (prev_error - error);
let crossing = prev_x + fraction * (point.position.x - prev_x);
if prev_error < 0.0 && error > 0.0 {
if near_crossing.is_none() {
near_crossing = Some(crossing);
}
} else {
far_crossing = Some(crossing);
}
}
}
previous = Some((point.position.x, error));
}
if let Some((last_x, last_error)) = previous {
if last_error == 0.0 {
if near_crossing.is_none() {
near_crossing = Some(last_x);
} else {
far_crossing = Some(last_x);
}
}
}
if near_crossing.is_none() && far_crossing.is_none() {
return Err(BallisticsError::from(
"Cannot find zero range: this angle never crosses the target height within the \
solved range (angle too shallow to reach it, or both crossings lie beyond \
the solver's max range)."
.to_string(),
));
}
Ok(ZeroCrossings {
near_m: near_crossing,
far_m: far_crossing,
})
}
pub fn equivalent_horizontal_range(
&self,
target_range_m: f64,
zero_distance_m: f64,
) -> Option<f64> {
if !target_range_m.is_finite() || !zero_distance_m.is_finite() {
return None;
}
if target_range_m <= zero_distance_m || target_range_m <= 0.0 {
return None;
}
fn path_y_at(points: &[TrajectoryPoint], distance_m: f64) -> Option<f64> {
match bracket_param(points.len(), |i| points[i].position.x, distance_m) {
Bracket::Below => Some(points[0].position.y),
Bracket::Above | Bracket::Degenerate => None,
Bracket::Inside { lo, t } => {
let hi = lo + 1;
Some(
points[lo].position.y
+ t * (points[hi].position.y - points[lo].position.y),
)
}
}
}
let mut inclined = self.clone();
inclined.inputs.enable_trajectory_sampling = false;
let inclined_result = inclined.solve().ok()?;
let los_height = inclined_result.line_of_sight_height_m;
let inclined_drop = los_height - path_y_at(&inclined_result.points, target_range_m)?;
let correction = inclined_drop / target_range_m;
if correction <= 0.0 {
return None;
}
let mut flat = self.clone();
flat.inputs.enable_trajectory_sampling = false;
flat.inputs.shooting_angle = 0.0;
let flat_result = flat.solve().ok()?;
let flat_correction_at = |range_m: f64| -> Option<f64> {
Some((los_height - path_y_at(&flat_result.points, range_m)?) / range_m)
};
let flat_terminal_m = flat_result.points.last().map(|p| p.position.x)?;
let mut low = zero_distance_m.max(1.0);
let mut high = target_range_m.min(flat_terminal_m);
if high <= low {
return None;
}
if flat_correction_at(low)? - correction > 0.0 {
return None; }
if flat_correction_at(high)? - correction < 0.0 {
return None; }
for _ in 0..60 {
let mid = 0.5 * (low + high);
let error = flat_correction_at(mid)? - correction;
if error.abs() == 0.0 {
return Some(mid);
}
if error < 0.0 {
low = mid;
} else {
high = mid;
}
if high - low < 0.01 {
break;
}
}
Some(0.5 * (low + high))
}
fn validate_for_solve(&self) -> Result<(), BallisticsError> {
let require_finite = |name: &str, value: f64| {
if value.is_finite() {
Ok(())
} else {
Err(BallisticsError::from(format!("{name} must be finite")))
}
};
let require_positive = |name: &str, value: f64| {
if value.is_finite() && value > 0.0 {
Ok(())
} else {
Err(BallisticsError::from(format!(
"{name} must be finite and greater than zero"
)))
}
};
if self.inputs.custom_drag_table.is_none() {
require_positive("bc_value", self.inputs.bc_value)?;
}
require_positive("bullet_mass", self.inputs.bullet_mass)?;
require_positive("bullet_diameter", self.inputs.bullet_diameter)?;
require_positive("muzzle_velocity", self.inputs.muzzle_velocity)?;
require_positive("cd_scale", self.inputs.cd_scale)?;
require_finite("muzzle_angle", self.inputs.muzzle_angle)?;
require_finite("azimuth_angle", self.inputs.azimuth_angle)?;
require_finite("shooting_angle", self.inputs.shooting_angle)?;
require_finite("cant_angle", self.inputs.cant_angle)?;
require_finite("muzzle_height", self.inputs.muzzle_height)?;
for (name, value) in [
("zero_poi_vertical_m", self.inputs.zero_poi_vertical_m),
("zero_poi_horizontal_m", self.inputs.zero_poi_horizontal_m),
] {
require_finite(name, value)?;
if value.abs() >= 1.0 {
return Err(BallisticsError::from(format!(
"{name} must be smaller than 1.0 m in magnitude (it is a linear POI \
offset at the zero range, in meters)"
)));
}
}
require_finite(
"sight_offset_lateral_m",
self.inputs.sight_offset_lateral_m,
)?;
if self.inputs.sight_offset_lateral_m.abs() >= 0.5 {
return Err(BallisticsError::from(
"sight_offset_lateral_m must be smaller than 0.5 m in magnitude (it is \
the lateral sight-to-bore mount offset, in meters)",
));
}
if !(self.inputs.ground_threshold.is_finite()
|| self.inputs.ground_threshold == f64::NEG_INFINITY)
{
return Err(BallisticsError::from(
"ground_threshold must be finite or negative infinity",
));
}
match &self.wind_sock {
Some(wind_sock) => wind_sock
.validate_segments()
.map_err(BallisticsError::from)?,
None => {
require_finite("wind.speed", self.wind.speed)?;
require_finite("wind.direction", self.wind.direction)?;
require_finite("wind.vertical_speed", self.wind.vertical_speed)?;
}
}
require_finite("atmosphere.temperature", self.atmosphere.temperature)?;
require_finite("atmosphere.pressure", self.atmosphere.pressure)?;
require_finite("atmosphere.humidity", self.atmosphere.humidity)?;
require_finite("atmosphere.altitude", self.atmosphere.altitude)?;
require_positive("max_range", self.max_range)?;
if !self.inputs.use_rk4 || !self.inputs.use_adaptive_rk45 {
require_positive("time_step", self.time_step)?;
}
if self.inputs.enable_trajectory_sampling {
require_finite("sight_height", self.inputs.sight_height)?;
require_positive("sample_interval", self.inputs.sample_interval)?;
projected_sample_count(self.max_range, self.inputs.sample_interval)?;
}
if self.inputs.drops_reference == DropsReference::Target {
require_finite("target_height", self.inputs.target_height)?;
if self.inputs.shooting_angle.cos() <= 1e-9 {
return Err(BallisticsError::from(
"drops reference 'target' is undefined for shooting angles at or beyond 90 degrees",
));
}
}
if self.inputs.enable_coriolis {
require_finite("shot_azimuth", self.inputs.shot_azimuth)?;
if let Some(latitude) = self.inputs.latitude {
require_finite("latitude", latitude)?;
}
}
Ok(())
}
fn validate_result_sanity(&self, result: &TrajectoryResult) -> Result<(), BallisticsError> {
let require_finite = |name: &str, value: f64| {
if value.is_finite() {
Ok(())
} else {
Err(BallisticsError::from(format!(
"trajectory result contains non-finite {name}"
)))
}
};
let require_non_negative = |name: &str, value: f64| {
if value >= 0.0 {
Ok(())
} else {
Err(BallisticsError::from(format!(
"trajectory result contains non-physical negative {name} ({value})"
)))
}
};
let require_indexed_finite = |collection: &str, index: usize, field: &str, value: f64| {
if value.is_finite() {
Ok(())
} else {
Err(BallisticsError::from(format!(
"trajectory result contains non-finite {collection}[{index}].{field}"
)))
}
};
let require_indexed_non_negative =
|collection: &str, index: usize, field: &str, value: f64| {
if value >= 0.0 {
Ok(())
} else {
Err(BallisticsError::from(format!(
"trajectory result contains non-physical negative {collection}[{index}].{field} ({value})"
)))
}
};
require_finite("max_range", result.max_range)?;
require_finite("max_height", result.max_height)?;
require_finite("time_of_flight", result.time_of_flight)?;
require_finite("impact_velocity", result.impact_velocity)?;
require_finite("impact_energy", result.impact_energy)?;
require_finite("projectile_mass_kg", result.projectile_mass_kg)?;
require_finite(
"line_of_sight_height_m",
result.line_of_sight_height_m,
)?;
require_finite(
"station_speed_of_sound_mps",
result.station_speed_of_sound_mps,
)?;
require_non_negative("max_range", result.max_range)?;
require_non_negative("time_of_flight", result.time_of_flight)?;
require_non_negative("impact_velocity", result.impact_velocity)?;
require_non_negative("impact_energy", result.impact_energy)?;
require_non_negative("projectile_mass_kg", result.projectile_mass_kg)?;
require_non_negative(
"station_speed_of_sound_mps",
result.station_speed_of_sound_mps,
)?;
for (index, point) in result.points.iter().enumerate() {
require_indexed_finite("points", index, "time", point.time)?;
require_indexed_finite("points", index, "position.x", point.position.x)?;
require_indexed_finite("points", index, "position.y", point.position.y)?;
require_indexed_finite("points", index, "position.z", point.position.z)?;
require_indexed_finite(
"points",
index,
"velocity_magnitude",
point.velocity_magnitude,
)?;
require_indexed_finite("points", index, "kinetic_energy", point.kinetic_energy)?;
require_indexed_non_negative("points", index, "time", point.time)?;
require_indexed_non_negative(
"points",
index,
"velocity_magnitude",
point.velocity_magnitude,
)?;
require_indexed_non_negative("points", index, "kinetic_energy", point.kinetic_energy)?;
}
if let Some(samples) = &result.sampled_points {
for (index, sample) in samples.iter().enumerate() {
require_indexed_finite("sampled_points", index, "distance_m", sample.distance_m)?;
require_indexed_finite("sampled_points", index, "drop_m", sample.drop_m)?;
require_indexed_finite(
"sampled_points",
index,
"wind_drift_m",
sample.wind_drift_m,
)?;
require_indexed_finite(
"sampled_points",
index,
"velocity_mps",
sample.velocity_mps,
)?;
require_indexed_finite("sampled_points", index, "energy_j", sample.energy_j)?;
require_indexed_finite("sampled_points", index, "time_s", sample.time_s)?;
}
}
for (name, value) in [
("min_pitch_damping", result.min_pitch_damping),
("transonic_mach", result.transonic_mach),
("max_yaw_angle", result.max_yaw_angle),
("max_precession_angle", result.max_precession_angle),
] {
if let Some(value) = value {
require_finite(name, value)?;
}
}
if let Some(state) = result.angular_state {
for (name, value) in [
("angular_state.pitch_angle", state.pitch_angle),
("angular_state.yaw_angle", state.yaw_angle),
("angular_state.pitch_rate", state.pitch_rate),
("angular_state.yaw_rate", state.yaw_rate),
("angular_state.precession_angle", state.precession_angle),
("angular_state.nutation_phase", state.nutation_phase),
] {
require_finite(name, value)?;
}
}
if let Some(jump) = result.aerodynamic_jump {
for (name, value) in [
("aerodynamic_jump.vertical_jump_moa", jump.vertical_jump_moa),
(
"aerodynamic_jump.horizontal_jump_moa",
jump.horizontal_jump_moa,
),
("aerodynamic_jump.jump_angle_rad", jump.jump_angle_rad),
(
"aerodynamic_jump.magnus_component_moa",
jump.magnus_component_moa,
),
("aerodynamic_jump.yaw_component_moa", jump.yaw_component_moa),
(
"aerodynamic_jump.stabilization_factor",
jump.stabilization_factor,
),
] {
require_finite(name, value)?;
}
}
Ok(())
}
fn validate_integration_state(
&self,
position: &Vector3<f64>,
velocity: &Vector3<f64>,
time: f64,
) -> Result<(), BallisticsError> {
if !(position.iter().all(|value| value.is_finite())
&& velocity.iter().all(|value| value.is_finite())
&& time.is_finite())
{
return Err(BallisticsError::from(
"trajectory integration produced a non-finite state (often from physically \
extreme inputs — e.g. an absurd bore/muzzle height placing the launch far \
from sea level, or a degenerate atmosphere; check those inputs, or set \
--altitude explicitly)",
));
}
let speed = velocity.magnitude();
let budget = self.speed_budget(time);
if speed > budget {
return Err(BallisticsError::from(format!(
"trajectory integration diverged: speed {speed:.3e} m/s at t={time:.6}s exceeds \
the physical budget of {budget:.3e} m/s"
)));
}
Ok(())
}
fn speed_budget(&self, time: f64) -> f64 {
let scalar_wind = self.wind.speed.abs() + self.wind.vertical_speed.abs();
let wind_bound = match &self.wind_sock {
Some(sock) => scalar_wind.max(sock.max_speed_mps()),
None => scalar_wind,
};
2.0 * (self.inputs.muzzle_velocity + wind_bound + 10.0)
+ crate::constants::G_ACCEL_MPS2 * time
}
fn push_trajectory_point(
&self,
points: &mut Vec<TrajectoryPoint>,
point: TrajectoryPoint,
) -> Result<(), BallisticsError> {
if points.len() >= self.max_trajectory_points {
return Err(BallisticsError::from(format!(
"trajectory point limit of {} exceeded",
self.max_trajectory_points
)));
}
points.push(point);
Ok(())
}
pub fn set_wind_segments(&mut self, segments: Vec<crate::wind::WindSegment>) {
self.wind_sock = if segments.is_empty() {
None
} else {
Some(crate::wind::WindSock::new(segments))
};
}
pub fn set_atmo_segments(&mut self, segments: Vec<crate::atmosphere::AtmoSegment>) {
self.atmo_sock = if segments.is_empty() {
None
} else {
Some(crate::atmosphere::AtmoSock::new(segments))
};
}
fn launch_angles_from(
&self,
aj: Option<&crate::aerodynamic_jump::AerodynamicJumpComponents>,
) -> (f64, f64) {
let (mut elev, mut azim) = (self.inputs.muzzle_angle, self.inputs.azimuth_angle);
if self.inputs.cant_angle != 0.0 {
let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
let (e0, a0) = (elev, azim);
elev = e0 * cos_c - a0 * sin_c;
azim = a0 * cos_c + e0 * sin_c;
}
match aj {
Some(c) => {
const MOA_PER_RAD: f64 = 3437.7467707849;
(
elev + c.vertical_jump_moa / MOA_PER_RAD,
azim + c.horizontal_jump_moa / MOA_PER_RAD,
)
}
None => (elev, azim),
}
}
fn aerodynamic_jump_components(
&self,
) -> Option<crate::aerodynamic_jump::AerodynamicJumpComponents> {
if !self.inputs.enable_aerodynamic_jump {
return None;
}
let diameter_m = self.inputs.bullet_diameter;
if !(self.inputs.twist_rate.is_finite()
&& self.inputs.twist_rate != 0.0
&& diameter_m.is_finite()
&& diameter_m > 0.0
&& self.inputs.bullet_length.is_finite()
&& self.inputs.bullet_length > 0.0
&& self.inputs.muzzle_velocity.is_finite())
{
return None;
}
let (_, _, temp_c, pressure_hpa) = self.resolved_atmosphere();
let sg = crate::stability::compute_stability_coefficient(
&self.inputs,
(self.atmosphere.altitude, temp_c, pressure_hpa, 0.0),
);
if !(sg.is_finite() && sg > 0.0) {
return None;
}
let length_calibers = self.inputs.bullet_length / diameter_m;
const MS_TO_MPH: f64 = 2.236_936_292_054_4;
let crosswind_from_right_mps = if let Some(sock) = &self.wind_sock {
-sock.vector_for_range_stateless(0.0)[2]
} else {
self.wind.speed * self.wind.direction.sin()
};
let crosswind_from_right_mph = crosswind_from_right_mps * MS_TO_MPH;
let vertical_jump_moa = crate::aerodynamic_jump::litz_crosswind_jump_moa(
sg,
length_calibers,
crosswind_from_right_mph,
self.inputs.is_twist_right,
);
if !vertical_jump_moa.is_finite() {
return None;
}
const MOA_PER_RAD: f64 = 3437.7467707849;
Some(crate::aerodynamic_jump::AerodynamicJumpComponents {
vertical_jump_moa,
horizontal_jump_moa: 0.0,
jump_angle_rad: vertical_jump_moa.abs() / MOA_PER_RAD,
magnus_component_moa: 0.0,
yaw_component_moa: 0.0,
stabilization_factor: (sg / 1.5).clamp(0.0, 1.0),
})
}
fn resolved_atmosphere(&self) -> (f64, f64, f64, f64) {
let (temp_c, pressure_hpa) = match self.station_atmosphere_resolution {
StationAtmosphereResolution::LegacyDefaultSentinels => {
crate::atmosphere::resolve_station_conditions(
self.atmosphere.temperature,
self.atmosphere.pressure,
self.atmosphere.altitude,
)
}
StationAtmosphereResolution::Authoritative => {
(self.atmosphere.temperature, self.atmosphere.pressure)
}
};
let (density, speed_of_sound) = crate::atmosphere::calculate_atmosphere(
self.atmosphere.altitude,
Some(temp_c),
Some(pressure_hpa),
self.atmosphere.humidity,
);
(density, speed_of_sound, temp_c, pressure_hpa)
}
fn precession_nutation_params(
&self,
velocity_mps: f64,
air_density_kg_m3: f64,
speed_of_sound_mps: f64,
) -> PrecessionNutationParams {
let (spin_inertia, transverse_inertia) = self.precession_nutation_inertias;
let spin_rate_rad_s = if self.inputs.twist_rate > 0.0 {
let velocity_fps = velocity_mps * 3.28084;
let twist_rate_ft = self.inputs.twist_rate / 12.0;
(velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI
} else {
0.0
};
PrecessionNutationParams {
mass_kg: self.inputs.bullet_mass,
caliber_m: self.inputs.bullet_diameter,
length_m: self.inputs.bullet_length,
spin_rate_rad_s,
spin_inertia,
transverse_inertia,
velocity_mps,
air_density_kg_m3,
mach: velocity_mps / speed_of_sound_mps,
pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
nutation_damping_factor: 0.05,
}
}
fn append_terminal_endpoint(
&self,
points: &mut Vec<TrajectoryPoint>,
post_position: Vector3<f64>,
post_velocity: Vector3<f64>,
post_time: f64,
max_height: &mut f64,
) -> Result<TrajectoryTermination, BallisticsError> {
let previous = points
.last()
.cloned()
.ok_or_else(|| BallisticsError::from("No trajectory points generated"))?;
let mut crossings = Vec::with_capacity(3);
if previous.position.x < self.max_range && post_position.x >= self.max_range {
let span = post_position.x - previous.position.x;
if span.is_finite() && span > 0.0 {
crossings.push((
(self.max_range - previous.position.x) / span,
TrajectoryTermination::MaxRange,
));
}
}
if self.inputs.ground_threshold.is_finite()
&& previous.position.y > self.inputs.ground_threshold
&& post_position.y <= self.inputs.ground_threshold
{
let span = post_position.y - previous.position.y;
if span.is_finite() && span < 0.0 {
crossings.push((
(self.inputs.ground_threshold - previous.position.y) / span,
TrajectoryTermination::GroundThreshold,
));
}
}
if previous.time < TRAJECTORY_TIME_LIMIT_S && post_time >= TRAJECTORY_TIME_LIMIT_S {
let span = post_time - previous.time;
if span.is_finite() && span > 0.0 {
crossings.push((
(TRAJECTORY_TIME_LIMIT_S - previous.time) / span,
TrajectoryTermination::TimeLimit,
));
}
}
let (fraction, termination) = crossings
.into_iter()
.filter(|(fraction, _)| fraction.is_finite() && (0.0..=1.0).contains(fraction))
.min_by(|left, right| {
let priority = |termination: TrajectoryTermination| match termination {
TrajectoryTermination::GroundThreshold => 0,
TrajectoryTermination::MaxRange => 1,
TrajectoryTermination::TimeLimit => 2,
TrajectoryTermination::VelocityFloor => 3,
};
left.0
.total_cmp(&right.0)
.then_with(|| priority(left.1).cmp(&priority(right.1)))
})
.ok_or_else(|| {
BallisticsError::from(
"trajectory integration stopped without crossing a supported boundary",
)
})?;
let mut position = previous.position + (post_position - previous.position) * fraction;
match termination {
TrajectoryTermination::MaxRange => position.x = self.max_range,
TrajectoryTermination::GroundThreshold => {
position.y = self.inputs.ground_threshold;
}
TrajectoryTermination::TimeLimit | TrajectoryTermination::VelocityFloor => {}
}
let velocity_magnitude = previous.velocity_magnitude
+ (post_velocity.magnitude() - previous.velocity_magnitude) * fraction;
let mut time = previous.time + (post_time - previous.time) * fraction;
if termination == TrajectoryTermination::TimeLimit {
time = TRAJECTORY_TIME_LIMIT_S;
}
let kinetic_energy =
0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
if position.y > *max_height {
*max_height = position.y;
}
let terminal_point = TrajectoryPoint {
time,
position,
velocity_magnitude,
kinetic_energy,
drag_coefficient: None,
};
if terminal_point.position.x < previous.position.x {
return Err(BallisticsError::from(
"trajectory terminal state reversed downrange before the crossed boundary",
));
}
if terminal_point.position.x == previous.position.x {
let last = points.last_mut().ok_or_else(|| {
BallisticsError::from("trajectory points disappeared during terminal finalization")
})?;
*last = terminal_point;
} else {
self.push_trajectory_point(points, terminal_point)?;
}
Ok(termination)
}
fn gravity_acceleration(&self) -> Vector3<f64> {
let theta = self.inputs.shooting_angle;
Vector3::new(
-crate::constants::G_ACCEL_MPS2 * theta.sin(),
-crate::constants::G_ACCEL_MPS2 * theta.cos(),
0.0,
)
}
fn get_wind_at_altitude(&self, altitude_m: f64) -> Vector3<f64> {
let model = match self.inputs.wind_shear_model.as_str() {
"logarithmic" => WindShearModel::Logarithmic,
"power_law" | "powerlaw" | "exponential" => WindShearModel::PowerLaw,
"ekman_spiral" | "ekman" => WindShearModel::EkmanSpiral,
"custom_layers" | "custom" => WindShearModel::CustomLayers,
_ => WindShearModel::PowerLaw,
};
let speed_ratio = crate::wind_shear::boundary_layer_speed_ratio(altitude_m, model);
crate::wind::wind_vector(self.wind.speed, self.wind.direction, 0.0) * speed_ratio
+ Vector3::new(0.0, self.wind.vertical_speed, 0.0)
}
pub fn solve(&self) -> Result<TrajectoryResult, BallisticsError> {
self.validate_for_solve()?;
let mut result = if self.inputs.use_rk4 {
if self.inputs.use_adaptive_rk45 {
self.solve_rk45()?
} else {
self.solve_rk4()?
}
} else {
self.solve_euler()?
};
self.apply_spin_drift(&mut result);
self.validate_result_sanity(&result)?;
Ok(result)
}
fn apply_spin_drift(&self, result: &mut TrajectoryResult) {
if !self.inputs.use_enhanced_spin_drift {
return;
}
let d_in = self.inputs.bullet_diameter / 0.0254; let m_gr = self.inputs.bullet_mass / crate::constants::GRAINS_TO_KG; let twist_in = self.inputs.twist_rate; if d_in <= 0.0 || m_gr <= 0.0 || twist_in <= 0.0 {
return;
}
let sg = self.effective_spin_drift_sg();
for p in result.points.iter_mut() {
if p.time <= 0.0 {
continue;
}
p.position.z +=
crate::spin_drift::litz_drift_meters(sg, p.time, self.inputs.is_twist_right);
}
if let Some(samples) = result.sampled_points.as_mut() {
for s in samples.iter_mut() {
if s.time_s <= 0.0 {
continue;
}
s.wind_drift_m +=
crate::spin_drift::litz_drift_meters(sg, s.time_s, self.inputs.is_twist_right);
}
}
}
fn effective_spin_drift_sg(&self) -> f64 {
let (_, _, temp_c, press_hpa) = self.resolved_atmosphere();
crate::spin_drift::effective_sg_from_inputs(&self.inputs, temp_c, press_hpa)
}
fn initial_position(&self) -> Vector3<f64> {
if self.inputs.cant_angle == 0.0 && self.inputs.sight_offset_lateral_m == 0.0 {
return Vector3::new(0.0, self.inputs.muzzle_height, 0.0);
}
let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
let sh = self.inputs.sight_height;
let off = self.inputs.sight_offset_lateral_m;
Vector3::new(
0.0,
self.inputs.muzzle_height + sh * (1.0 - cos_c) + off * sin_c,
-sh * sin_c - off * cos_c,
)
}
fn build_sampled_points(
&self,
points: &[TrajectoryPoint],
max_height: f64,
transonic_distances: Vec<f64>,
mach_transitions: &MachTransitionTracker,
) -> Result<Option<Vec<TrajectorySample>>, BallisticsError> {
if !self.inputs.enable_trajectory_sampling {
return Ok(None);
}
let last_point = points.last().ok_or("No trajectory points generated")?;
let trajectory_data = TrajectoryData {
times: points.iter().map(|p| p.time).collect(),
positions: points.iter().map(|p| p.position).collect(),
velocities: points
.iter()
.map(|p| {
Vector3::new(0.0, 0.0, p.velocity_magnitude)
})
.collect(),
transonic_distances, mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
};
let sight_position_m = self.inputs.muzzle_height + self.inputs.sight_height;
let target_reference = self.inputs.drops_reference == DropsReference::Target;
let target_vertical_height_m = if target_reference && self.inputs.target_height != 0.0 {
self.inputs.target_height
} else {
sight_position_m
};
let outputs = TrajectoryOutputs {
target_distance_horiz_m: last_point.position.x, target_vertical_height_m,
time_of_flight_s: last_point.time,
max_ord_dist_horiz_m: max_height,
sight_height_m: sight_position_m,
};
let mut samples = sample_trajectory(
&trajectory_data,
&outputs,
self.inputs.sample_interval,
self.inputs.bullet_mass,
)?;
if target_reference {
let cos_theta = self.inputs.shooting_angle.cos();
for sample in &mut samples {
sample.drop_m /= cos_theta;
}
}
Ok(Some(samples))
}
fn solve_euler(&self) -> Result<TrajectoryResult, BallisticsError> {
let mut time = 0.0;
let mut position = self.initial_position();
let aj_components = self.aerodynamic_jump_components();
let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
let mut velocity = Vector3::new(
horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
let mut points = Vec::new();
let mut max_height = position.y;
let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
let mut mach_transitions = MachTransitionTracker::default();
let mut angular_state = if self.inputs.enable_precession_nutation {
Some(AngularState {
pitch_angle: 0.001, yaw_angle: 0.001,
pitch_rate: 0.0,
yaw_rate: 0.0,
precession_angle: 0.0,
nutation_phase: 0.0,
})
} else {
None
};
let mut max_yaw_angle = 0.0;
let mut max_precession_angle = 0.0;
let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
self.resolved_atmosphere();
let base_ratio = air_density / 1.225;
let wind_vector =
crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
self.inputs.bullet_model.as_deref().unwrap_or("default"),
);
while position.x < self.max_range
&& position.y > self.inputs.ground_threshold
&& time < TRAJECTORY_TIME_LIMIT_S
{
let velocity_magnitude = velocity.magnitude();
let kinetic_energy =
0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
self.push_trajectory_point(
&mut points,
TrajectoryPoint {
time,
position,
velocity_magnitude,
kinetic_energy,
drag_coefficient: None,
},
)?;
{
let mach_here = if speed_of_sound > 0.0 {
velocity_magnitude / speed_of_sound
} else {
0.0
};
mach_transitions.record_downward_crossings(
mach_here,
position.x,
&mut transonic_distances,
);
}
if position.y > max_height {
max_height = position.y;
}
if self.inputs.enable_pitch_damping {
let mach = velocity_magnitude / speed_of_sound;
if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
transonic_mach = Some(mach);
}
let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
if pitch_damping < min_pitch_damping {
min_pitch_damping = pitch_damping;
}
}
if self.inputs.enable_precession_nutation {
if let Some(ref mut state) = angular_state {
let velocity_magnitude = velocity.magnitude();
let params = self.precession_nutation_params(
velocity_magnitude,
air_density,
speed_of_sound,
);
*state = calculate_combined_angular_motion(
¶ms,
state,
time,
self.time_step,
0.001, );
if state.yaw_angle.abs() > max_yaw_angle {
max_yaw_angle = state.yaw_angle.abs();
}
if state.precession_angle.abs() > max_precession_angle {
max_precession_angle = state.precession_angle.abs();
}
}
}
let acceleration = self.calculate_acceleration(
&position,
&velocity,
&wind_vector,
(resolved_temp_c, resolved_press_hpa, base_ratio),
);
velocity += acceleration * self.time_step;
position += velocity * self.time_step;
time += self.time_step;
self.validate_integration_state(&position, &velocity, time)?;
}
let termination =
self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
self.annotate_drag_coefficients(&mut points, speed_of_sound);
let last_point = points.last().ok_or("No trajectory points generated")?;
let sampled_points = self.build_sampled_points(
&points,
max_height,
transonic_distances,
&mach_transitions,
)?;
Ok(TrajectoryResult {
max_range: last_point.position.x, max_height,
time_of_flight: last_point.time,
impact_velocity: last_point.velocity_magnitude,
impact_energy: last_point.kinetic_energy,
projectile_mass_kg: self.inputs.bullet_mass,
line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
station_speed_of_sound_mps: speed_of_sound,
termination,
points,
sampled_points,
min_pitch_damping: if self.inputs.enable_pitch_damping {
Some(min_pitch_damping)
} else {
None
},
transonic_mach,
angular_state,
max_yaw_angle: if self.inputs.enable_precession_nutation {
Some(max_yaw_angle)
} else {
None
},
max_precession_angle: if self.inputs.enable_precession_nutation {
Some(max_precession_angle)
} else {
None
},
aerodynamic_jump: aj_components,
mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
})
}
fn solve_rk4(&self) -> Result<TrajectoryResult, BallisticsError> {
let mut time = 0.0;
let mut position = self.initial_position();
let aj_components = self.aerodynamic_jump_components();
let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
let mut velocity = Vector3::new(
horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
let mut points = Vec::new();
let mut max_height = position.y;
let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
let mut mach_transitions = MachTransitionTracker::default();
let mut angular_state = if self.inputs.enable_precession_nutation {
Some(AngularState {
pitch_angle: 0.001, yaw_angle: 0.001,
pitch_rate: 0.0,
yaw_rate: 0.0,
precession_angle: 0.0,
nutation_phase: 0.0,
})
} else {
None
};
let mut max_yaw_angle = 0.0;
let mut max_precession_angle = 0.0;
let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
self.resolved_atmosphere();
let base_ratio = air_density / 1.225;
let wind_vector =
crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
self.inputs.bullet_model.as_deref().unwrap_or("default"),
);
while position.x < self.max_range
&& position.y > self.inputs.ground_threshold
&& time < TRAJECTORY_TIME_LIMIT_S
{
let velocity_magnitude = velocity.magnitude();
let kinetic_energy =
0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
self.push_trajectory_point(
&mut points,
TrajectoryPoint {
time,
position,
velocity_magnitude,
kinetic_energy,
drag_coefficient: None,
},
)?;
{
let mach_here = if speed_of_sound > 0.0 {
velocity_magnitude / speed_of_sound
} else {
0.0
};
mach_transitions.record_downward_crossings(
mach_here,
position.x,
&mut transonic_distances,
);
}
if position.y > max_height {
max_height = position.y;
}
if self.inputs.enable_pitch_damping {
let mach = velocity_magnitude / speed_of_sound;
if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
transonic_mach = Some(mach);
}
let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
if pitch_damping < min_pitch_damping {
min_pitch_damping = pitch_damping;
}
}
if self.inputs.enable_precession_nutation {
if let Some(ref mut state) = angular_state {
let velocity_magnitude = velocity.magnitude();
let params = self.precession_nutation_params(
velocity_magnitude,
air_density,
speed_of_sound,
);
*state = calculate_combined_angular_motion(
¶ms,
state,
time,
self.time_step,
0.001, );
if state.yaw_angle.abs() > max_yaw_angle {
max_yaw_angle = state.yaw_angle.abs();
}
if state.precession_angle.abs() > max_precession_angle {
max_precession_angle = state.precession_angle.abs();
}
}
}
let dt = self.time_step;
let acc1 = self.calculate_acceleration(
&position,
&velocity,
&wind_vector,
(resolved_temp_c, resolved_press_hpa, base_ratio),
);
let pos2 = position + velocity * (dt * 0.5);
let vel2 = velocity + acc1 * (dt * 0.5);
let acc2 = self.calculate_acceleration(
&pos2,
&vel2,
&wind_vector,
(resolved_temp_c, resolved_press_hpa, base_ratio),
);
let pos3 = position + vel2 * (dt * 0.5);
let vel3 = velocity + acc2 * (dt * 0.5);
let acc3 = self.calculate_acceleration(
&pos3,
&vel3,
&wind_vector,
(resolved_temp_c, resolved_press_hpa, base_ratio),
);
let pos4 = position + vel3 * dt;
let vel4 = velocity + acc3 * dt;
let acc4 = self.calculate_acceleration(
&pos4,
&vel4,
&wind_vector,
(resolved_temp_c, resolved_press_hpa, base_ratio),
);
position += (velocity + vel2 * 2.0 + vel3 * 2.0 + vel4) * (dt / 6.0);
velocity += (acc1 + acc2 * 2.0 + acc3 * 2.0 + acc4) * (dt / 6.0);
time += dt;
self.validate_integration_state(&position, &velocity, time)?;
}
let termination =
self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
self.annotate_drag_coefficients(&mut points, speed_of_sound);
let last_point = points.last().ok_or("No trajectory points generated")?;
let sampled_points = self.build_sampled_points(
&points,
max_height,
transonic_distances,
&mach_transitions,
)?;
Ok(TrajectoryResult {
max_range: last_point.position.x, max_height,
time_of_flight: last_point.time,
impact_velocity: last_point.velocity_magnitude,
impact_energy: last_point.kinetic_energy,
projectile_mass_kg: self.inputs.bullet_mass,
line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
station_speed_of_sound_mps: speed_of_sound,
termination,
points,
sampled_points,
min_pitch_damping: if self.inputs.enable_pitch_damping {
Some(min_pitch_damping)
} else {
None
},
transonic_mach,
angular_state,
max_yaw_angle: if self.inputs.enable_precession_nutation {
Some(max_yaw_angle)
} else {
None
},
max_precession_angle: if self.inputs.enable_precession_nutation {
Some(max_precession_angle)
} else {
None
},
aerodynamic_jump: aj_components,
mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
})
}
fn solve_rk45(&self) -> Result<TrajectoryResult, BallisticsError> {
let mut time = 0.0;
let mut position = self.initial_position();
let aj_components = self.aerodynamic_jump_components();
let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
let mut velocity = Vector3::new(
horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
let mut points = Vec::new();
let mut max_height = position.y;
let mut dt = 0.001;
let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
self.resolved_atmosphere();
let base_ratio = air_density / 1.225;
let wind_vector =
crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
let mut transonic_distances: Vec<f64> = Vec::new();
let mut mach_transitions = MachTransitionTracker::default();
let mut min_pitch_damping = f64::INFINITY;
let mut transonic_mach: Option<f64> = None;
let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
self.inputs.bullet_model.as_deref().unwrap_or("default"),
);
let mut angular_state = if self.inputs.enable_precession_nutation {
Some(AngularState {
pitch_angle: 0.001,
yaw_angle: 0.001,
pitch_rate: 0.0,
yaw_rate: 0.0,
precession_angle: 0.0,
nutation_phase: 0.0,
})
} else {
None
};
let mut max_yaw_angle = 0.0;
let mut max_precession_angle = 0.0;
while position.x < self.max_range
&& position.y > self.inputs.ground_threshold
&& time < TRAJECTORY_TIME_LIMIT_S
{
let velocity_magnitude = velocity.magnitude();
let kinetic_energy = 0.5 * self.inputs.bullet_mass * velocity_magnitude.powi(2);
self.push_trajectory_point(
&mut points,
TrajectoryPoint {
time,
position,
velocity_magnitude,
kinetic_energy,
drag_coefficient: None,
},
)?;
{
let mach_here = if speed_of_sound > 0.0 {
velocity_magnitude / speed_of_sound
} else {
0.0
};
mach_transitions.record_downward_crossings(
mach_here,
position.x,
&mut transonic_distances,
);
}
if position.y > max_height {
max_height = position.y;
}
if self.inputs.enable_pitch_damping {
let mach = velocity_magnitude / speed_of_sound;
if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
transonic_mach = Some(mach);
}
let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
if pitch_damping < min_pitch_damping {
min_pitch_damping = pitch_damping;
}
}
let accepted_step = self.adaptive_rk45_step(
&position,
&velocity,
dt,
&wind_vector,
(resolved_temp_c, resolved_press_hpa, base_ratio),
);
debug_assert!(
accepted_step.error <= RK45_TOLERANCE || accepted_step.used_dt <= RK45_MIN_DT
);
if self.inputs.enable_precession_nutation {
if let Some(ref mut state) = angular_state {
let params = self.precession_nutation_params(
velocity_magnitude,
air_density,
speed_of_sound,
);
*state = calculate_combined_angular_motion(
¶ms,
state,
time,
accepted_step.used_dt,
0.001,
);
if state.yaw_angle.abs() > max_yaw_angle {
max_yaw_angle = state.yaw_angle.abs();
}
if state.precession_angle.abs() > max_precession_angle {
max_precession_angle = state.precession_angle.abs();
}
}
}
position = accepted_step.position;
velocity = accepted_step.velocity;
time += accepted_step.used_dt;
self.validate_integration_state(&position, &velocity, time)?;
dt = accepted_step.next_dt;
}
if points.is_empty() {
return Err(BallisticsError::from("No trajectory points calculated"));
}
let termination =
self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
self.annotate_drag_coefficients(&mut points, speed_of_sound);
let last_point = points.last().unwrap();
let sampled_points = self.build_sampled_points(
&points,
max_height,
transonic_distances,
&mach_transitions,
)?;
Ok(TrajectoryResult {
max_range: last_point.position.x, max_height,
time_of_flight: last_point.time,
impact_velocity: last_point.velocity_magnitude,
impact_energy: last_point.kinetic_energy,
projectile_mass_kg: self.inputs.bullet_mass,
line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
station_speed_of_sound_mps: speed_of_sound,
termination,
points,
sampled_points,
min_pitch_damping: if self.inputs.enable_pitch_damping {
Some(min_pitch_damping)
} else {
None
},
transonic_mach,
angular_state,
max_yaw_angle: if self.inputs.enable_precession_nutation {
Some(max_yaw_angle)
} else {
None
},
max_precession_angle: if self.inputs.enable_precession_nutation {
Some(max_precession_angle)
} else {
None
},
aerodynamic_jump: aj_components,
mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
})
}
fn adaptive_rk45_step(
&self,
position: &Vector3<f64>,
velocity: &Vector3<f64>,
initial_dt: f64,
wind_vector: &Vector3<f64>,
resolved_atmo: (f64, f64, f64),
) -> Rk45AcceptedStep {
let mut trial_dt = initial_dt;
loop {
let trial = self.rk45_step(
position,
velocity,
trial_dt,
wind_vector,
RK45_TOLERANCE,
resolved_atmo,
);
let next_dt = if trial.suggested_dt.is_finite() {
(RK45_SAFETY_FACTOR * trial.suggested_dt).clamp(RK45_MIN_DT, RK45_MAX_DT)
} else {
RK45_MIN_DT
};
if trial.error <= RK45_TOLERANCE || trial_dt <= RK45_MIN_DT {
return Rk45AcceptedStep {
position: trial.position,
velocity: trial.velocity,
used_dt: trial_dt,
next_dt,
error: trial.error,
};
}
trial_dt = next_dt;
}
}
fn rk45_step(
&self,
position: &Vector3<f64>,
velocity: &Vector3<f64>,
dt: f64,
wind_vector: &Vector3<f64>,
tolerance: f64,
resolved_atmo: (f64, f64, f64), ) -> Rk45Trial {
const A21: f64 = 1.0 / 5.0;
const A31: f64 = 3.0 / 40.0;
const A32: f64 = 9.0 / 40.0;
const A41: f64 = 44.0 / 45.0;
const A42: f64 = -56.0 / 15.0;
const A43: f64 = 32.0 / 9.0;
const A51: f64 = 19372.0 / 6561.0;
const A52: f64 = -25360.0 / 2187.0;
const A53: f64 = 64448.0 / 6561.0;
const A54: f64 = -212.0 / 729.0;
const A61: f64 = 9017.0 / 3168.0;
const A62: f64 = -355.0 / 33.0;
const A63: f64 = 46732.0 / 5247.0;
const A64: f64 = 49.0 / 176.0;
const A65: f64 = -5103.0 / 18656.0;
const A71: f64 = 35.0 / 384.0;
const A73: f64 = 500.0 / 1113.0;
const A74: f64 = 125.0 / 192.0;
const A75: f64 = -2187.0 / 6784.0;
const A76: f64 = 11.0 / 84.0;
const B1: f64 = 35.0 / 384.0;
const B3: f64 = 500.0 / 1113.0;
const B4: f64 = 125.0 / 192.0;
const B5: f64 = -2187.0 / 6784.0;
const B6: f64 = 11.0 / 84.0;
const B1_ERR: f64 = 5179.0 / 57600.0;
const B3_ERR: f64 = 7571.0 / 16695.0;
const B4_ERR: f64 = 393.0 / 640.0;
const B5_ERR: f64 = -92097.0 / 339200.0;
const B6_ERR: f64 = 187.0 / 2100.0;
const B7_ERR: f64 = 1.0 / 40.0;
let k1_v = self.calculate_acceleration(position, velocity, wind_vector, resolved_atmo);
let k1_p = *velocity;
let p2 = position + dt * A21 * k1_p;
let v2 = velocity + dt * A21 * k1_v;
let k2_v = self.calculate_acceleration(&p2, &v2, wind_vector, resolved_atmo);
let k2_p = v2;
let p3 = position + dt * (A31 * k1_p + A32 * k2_p);
let v3 = velocity + dt * (A31 * k1_v + A32 * k2_v);
let k3_v = self.calculate_acceleration(&p3, &v3, wind_vector, resolved_atmo);
let k3_p = v3;
let p4 = position + dt * (A41 * k1_p + A42 * k2_p + A43 * k3_p);
let v4 = velocity + dt * (A41 * k1_v + A42 * k2_v + A43 * k3_v);
let k4_v = self.calculate_acceleration(&p4, &v4, wind_vector, resolved_atmo);
let k4_p = v4;
let p5 = position + dt * (A51 * k1_p + A52 * k2_p + A53 * k3_p + A54 * k4_p);
let v5 = velocity + dt * (A51 * k1_v + A52 * k2_v + A53 * k3_v + A54 * k4_v);
let k5_v = self.calculate_acceleration(&p5, &v5, wind_vector, resolved_atmo);
let k5_p = v5;
let p6 = position + dt * (A61 * k1_p + A62 * k2_p + A63 * k3_p + A64 * k4_p + A65 * k5_p);
let v6 = velocity + dt * (A61 * k1_v + A62 * k2_v + A63 * k3_v + A64 * k4_v + A65 * k5_v);
let k6_v = self.calculate_acceleration(&p6, &v6, wind_vector, resolved_atmo);
let k6_p = v6;
let p7 = position + dt * (A71 * k1_p + A73 * k3_p + A74 * k4_p + A75 * k5_p + A76 * k6_p);
let v7 = velocity + dt * (A71 * k1_v + A73 * k3_v + A74 * k4_v + A75 * k5_v + A76 * k6_v);
let k7_v = self.calculate_acceleration(&p7, &v7, wind_vector, resolved_atmo);
let k7_p = v7;
let new_pos = position + dt * (B1 * k1_p + B3 * k3_p + B4 * k4_p + B5 * k5_p + B6 * k6_p);
let new_vel = velocity + dt * (B1 * k1_v + B3 * k3_v + B4 * k4_v + B5 * k5_v + B6 * k6_v);
let pos_err = position
+ dt * (B1_ERR * k1_p
+ B3_ERR * k3_p
+ B4_ERR * k4_p
+ B5_ERR * k5_p
+ B6_ERR * k6_p
+ B7_ERR * k7_p);
let vel_err = velocity
+ dt * (B1_ERR * k1_v
+ B3_ERR * k3_v
+ B4_ERR * k4_v
+ B5_ERR * k5_v
+ B6_ERR * k6_v
+ B7_ERR * k7_v);
let error = cli_rk45_error_norm(position, velocity, &new_pos, &new_vel, &pos_err, &vel_err);
let dt_new = if error < tolerance {
dt * (tolerance / error).powf(0.2).min(2.0)
} else {
dt * (tolerance / error).powf(0.25).max(0.1)
};
Rk45Trial {
position: new_pos,
velocity: new_vel,
suggested_dt: dt_new,
error,
}
}
fn apply_cluster_bc_correction(&self, base_bc: f64, velocity_fps: f64) -> f64 {
if let Some(ref cluster_bc) = self.cluster_bc {
cluster_bc.apply_correction_for_drag_model(
base_bc,
self.inputs.caliber_inches,
self.inputs.weight_grains,
velocity_fps,
self.inputs.bc_type,
)
} else {
base_bc
}
}
fn calculate_acceleration(
&self,
position: &Vector3<f64>,
velocity: &Vector3<f64>,
wind_vector: &Vector3<f64>,
resolved_atmo: (f64, f64, f64), ) -> Vector3<f64> {
let actual_wind = if let Some(ref sock) = self.wind_sock {
sock.vector_for_range_stateless(position.x)
} else if self.inputs.enable_wind_shear {
self.get_wind_at_altitude(position.y)
} else {
*wind_vector
};
let actual_wind =
crate::derivatives::level_vector_to_shot_frame(actual_wind, self.inputs.shooting_angle);
let relative_velocity = velocity - actual_wind;
let velocity_magnitude = relative_velocity.magnitude();
if velocity_magnitude < 0.001 {
return self.gravity_acceleration();
}
let (base_temp_c, base_press_hpa, station_ratio) = resolved_atmo;
let (drag_base_temp_c, drag_base_press_hpa, drag_base_ratio, drag_humidity_percent) =
if let Some(ref sock) = self.atmo_sock {
let (zone_temp_c, zone_press_hpa, zone_humidity) = sock.atmo_for_range(position.x);
let zone_base_ratio = crate::atmosphere::calculate_air_density_cimp(
zone_temp_c,
zone_press_hpa,
zone_humidity,
) / 1.225;
(zone_temp_c, zone_press_hpa, zone_base_ratio, zone_humidity)
} else {
(
base_temp_c,
base_press_hpa,
station_ratio,
self.atmosphere.humidity,
)
};
let local_alt = crate::atmosphere::shot_frame_altitude(
self.atmosphere.altitude,
position.x,
position.y,
self.inputs.shooting_angle,
);
let (air_density, speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
local_alt,
self.atmosphere.altitude,
drag_base_temp_c,
drag_base_press_hpa,
drag_base_ratio,
drag_humidity_percent,
);
let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
let velocity_fps = velocity_magnitude * 3.28084;
let cd_to_retard = crate::constants::CD_TO_RETARD;
let standard_factor = cd * cd_to_retard;
let density_scale = air_density / 1.225;
let a_drag_ft_s2 =
(velocity_fps * velocity_fps) * standard_factor * density_scale / retard_denom;
let a_drag_m_s2 = a_drag_ft_s2 * 0.3048;
let drag_acceleration = -a_drag_m_s2 * (relative_velocity / velocity_magnitude);
let mut accel = drag_acceleration + self.gravity_acceleration();
if self.inputs.enable_coriolis {
if let Some(lat_deg) = self.inputs.latitude {
let omega_earth = 7.2921159e-5_f64; let lat = lat_deg.to_radians();
let az = self.inputs.shot_azimuth; let omega = Vector3::new(
omega_earth * lat.cos() * az.cos(), omega_earth * lat.sin(), -omega_earth * lat.cos() * az.sin(), );
let omega = crate::derivatives::level_vector_to_shot_frame(
omega,
self.inputs.shooting_angle,
);
accel += -2.0 * omega.cross(velocity);
}
}
if self.inputs.enable_magnus
&& !self.inputs.use_enhanced_spin_drift
&& self.inputs.bullet_diameter > 0.0
&& self.inputs.twist_rate > 0.0
{
let diameter_m = self.inputs.bullet_diameter;
let (spin_rad_s, spin_param) = crate::spin_drift::calculate_magnus_spin_state(
self.inputs.muzzle_velocity,
velocity_magnitude,
self.inputs.twist_rate,
diameter_m,
);
let mach = velocity_magnitude / speed_of_sound;
let d_in = self.inputs.bullet_diameter / 0.0254;
let m_gr = self.inputs.bullet_mass / crate::constants::GRAINS_TO_KG;
let l_in = if self.inputs.bullet_length > 0.0 {
self.inputs.bullet_length / 0.0254
} else {
let est_m = crate::stability::estimate_bullet_length_m(
self.inputs.bullet_diameter,
self.inputs.bullet_mass,
);
if est_m > 0.0 {
est_m / 0.0254
} else {
4.5 * d_in
}
};
let sg = crate::spin_drift::calculate_dynamic_stability(
m_gr,
velocity_magnitude,
spin_rad_s,
d_in,
l_in,
air_density,
);
let (yaw_rad, _) = crate::spin_drift::calculate_yaw_of_repose(
sg,
velocity_magnitude,
spin_rad_s,
0.0, 0.0, air_density,
d_in,
l_in,
m_gr,
mach,
"match",
false,
);
let c_np = crate::derivatives::calculate_magnus_moment_coefficient(mach);
let area = std::f64::consts::PI * (diameter_m / 2.0).powi(2);
let magnus_force = 0.5
* air_density
* velocity_magnitude.powi(2)
* area
* c_np
* spin_param
* yaw_rad.sin();
if magnus_force.abs() > 1e-12 {
if let Some(dir) = crate::derivatives::yaw_of_repose_magnus_direction(
relative_velocity,
self.gravity_acceleration(),
self.inputs.is_twist_right,
) {
accel += (magnus_force / self.inputs.bullet_mass) * dir;
}
}
}
accel
}
fn drag_terms(&self, velocity_magnitude: f64, speed_of_sound: f64) -> (f64, f64) {
let cd = self.calculate_drag_coefficient(velocity_magnitude, speed_of_sound);
let velocity_fps = velocity_magnitude * 3.28084;
let (base_bc, bc_from_segments) = if let Some(segments) = self
.inputs
.bc_segments_data
.as_ref()
.filter(|segments| self.inputs.use_bc_segments && !segments.is_empty())
{
(
crate::bc_estimation::velocity_segment_bc(
velocity_fps,
segments,
self.inputs.bc_value,
),
true,
)
} else if let Some(segments) = self
.inputs
.bc_segments
.as_ref()
.filter(|segments| !segments.is_empty())
{
(
crate::derivatives::interpolated_bc(
velocity_magnitude / speed_of_sound,
segments,
Some(&self.inputs),
),
true,
)
} else {
(self.inputs.bc_value, false)
};
let effective_bc = if bc_from_segments {
base_bc
} else {
self.apply_cluster_bc_correction(base_bc, velocity_fps)
};
let effective_bc = effective_bc.max(1e-6);
let retard_denom = if self.inputs.custom_drag_table.is_some() {
self.inputs.custom_drag_denominator(effective_bc)
} else {
effective_bc
};
(cd, retard_denom)
}
pub fn effective_drag_coefficient(
&self,
velocity_magnitude: f64,
speed_of_sound: f64,
) -> Option<f64> {
if !velocity_magnitude.is_finite() || speed_of_sound <= 1e-9 {
return None;
}
let sectional_density = self.inputs.sectional_density_lb_in2()?;
let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
if retard_denom <= 0.0 {
return None;
}
let effective = cd * sectional_density / retard_denom;
effective.is_finite().then_some(effective)
}
fn annotate_drag_coefficients(&self, points: &mut [TrajectoryPoint], speed_of_sound: f64) {
for point in points.iter_mut() {
point.drag_coefficient =
self.effective_drag_coefficient(point.velocity_magnitude, speed_of_sound);
}
}
fn calculate_drag_coefficient(&self, velocity: f64, speed_of_sound: f64) -> f64 {
let mach = velocity / speed_of_sound;
if let Some(ref table) = self.inputs.custom_drag_table {
return table.interpolate(mach) * self.inputs.cd_scale;
}
crate::drag::get_drag_coefficient(mach, &self.inputs.bc_type)
}
}
#[derive(Debug, Clone)]
pub struct MonteCarloParams {
pub num_simulations: usize,
pub velocity_std_dev: f64,
pub angle_std_dev: f64,
pub bc_std_dev: f64,
pub wind_speed_std_dev: f64,
pub target_distance: Option<f64>,
pub base_wind_speed: f64,
pub base_wind_direction: f64,
pub azimuth_std_dev: f64, }
impl Default for MonteCarloParams {
fn default() -> Self {
Self {
num_simulations: 1000,
velocity_std_dev: 1.0,
angle_std_dev: 0.001,
bc_std_dev: 0.01,
wind_speed_std_dev: 1.0,
target_distance: None,
base_wind_speed: 0.0,
base_wind_direction: 0.0,
azimuth_std_dev: 0.001, }
}
}
#[derive(Debug, Clone)]
pub struct MonteCarloResults {
pub ranges: Vec<f64>,
pub impact_velocities: Vec<f64>,
pub impact_positions: Vec<Vector3<f64>>,
}
pub const DEFAULT_HIT_RADIUS_M: f64 = 0.3;
pub const TARGET_NOT_REACHED_SENTINEL_M: f64 = -1.0e9;
impl MonteCarloResults {
pub fn position_reached_target(position: &Vector3<f64>) -> bool {
position.iter().all(|component| component.is_finite())
&& position.y != TARGET_NOT_REACHED_SENTINEL_M
}
pub fn target_arrival_count(&self) -> usize {
self.impact_positions
.iter()
.filter(|position| Self::position_reached_target(position))
.count()
}
pub fn target_shortfall_fraction(&self) -> f64 {
if self.impact_positions.is_empty() {
return 0.0;
}
(self.impact_positions.len() - self.target_arrival_count()) as f64
/ self.impact_positions.len() as f64
}
pub fn target_plane_cep(&self) -> Option<f64> {
let mut radial_misses: Vec<f64> = self
.impact_positions
.iter()
.filter(|position| Self::position_reached_target(position))
.map(Vector3::norm)
.filter(|miss| miss.is_finite())
.collect();
radial_misses.sort_by(f64::total_cmp);
if radial_misses.is_empty() {
None
} else {
Some(radial_misses[radial_misses.len() / 2])
}
}
pub fn hit_probability(&self, hit_radius_m: f64) -> f64 {
if self.impact_positions.is_empty() {
return 0.0;
}
let hits = self
.impact_positions
.iter()
.filter(|position| Self::position_is_hit(position, hit_radius_m))
.count();
hits as f64 / self.impact_positions.len() as f64
}
pub fn position_is_hit(position: &Vector3<f64>, hit_radius_m: f64) -> bool {
Self::position_reached_target(position) && position.norm() < hit_radius_m
}
pub fn hit_probability_wilson(
&self,
hit_radius_m: f64,
level: ConfidenceLevel,
) -> (f64, (f64, f64), u64) {
let trials = self.impact_positions.len() as u64;
let hits = self
.impact_positions
.iter()
.filter(|position| Self::position_is_hit(position, hit_radius_m))
.count() as u64;
(
self.hit_probability(hit_radius_m),
wilson_interval(hits, trials, level),
trials,
)
}
pub fn rect_hit_probability(&self, width_m: f64, height_m: f64) -> f64 {
let dimensions_invalid = width_m.is_nan()
|| width_m <= 0.0
|| height_m.is_nan()
|| height_m <= 0.0;
if self.impact_positions.is_empty() || dimensions_invalid {
return 0.0;
}
let half_width = width_m / 2.0;
let half_height = height_m / 2.0;
let hits = self
.impact_positions
.iter()
.filter(|position| {
Self::position_reached_target(position)
&& position.z.abs() <= half_width
&& position.y.abs() <= half_height
})
.count();
hits as f64 / self.impact_positions.len() as f64
}
}
fn wind_from_signed_speed_sample(
signed_speed: f64,
sampled_direction: f64,
vertical_speed: f64,
) -> WindConditions {
if signed_speed < 0.0 {
WindConditions {
speed: -signed_speed,
direction: sampled_direction + std::f64::consts::PI,
vertical_speed,
}
} else {
WindConditions {
speed: signed_speed,
direction: sampled_direction,
vertical_speed,
}
}
}
struct MonteCarloWindSampler {
speed: rand_distr::Normal<f64>,
direction: rand_distr::Normal<f64>,
vertical_speed: f64,
}
impl MonteCarloWindSampler {
fn new(
base_wind: &WindConditions,
wind_speed_std_dev: f64,
wind_direction_std_dev: f64,
) -> Result<Self, BallisticsError> {
use rand_distr::Normal;
if !wind_direction_std_dev.is_finite() || wind_direction_std_dev < 0.0 {
return Err("Wind direction standard deviation must be finite and non-negative".into());
}
let speed = Normal::new(base_wind.speed, wind_speed_std_dev)
.map_err(|e| format!("Invalid wind speed distribution: {e}"))?;
let direction = Normal::new(base_wind.direction, wind_direction_std_dev)
.map_err(|e| format!("Invalid wind direction distribution: {e}"))?;
Ok(Self { speed, direction, vertical_speed: base_wind.vertical_speed })
}
fn sample<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> WindConditions {
use rand_distr::Distribution;
wind_from_signed_speed_sample(
self.speed.sample(rng),
self.direction.sample(rng),
self.vertical_speed,
)
}
}
#[derive(Debug, Clone, Copy)]
struct TrialOutcome {
range: f64,
impact_velocity: f64,
impact_position: Vector3<f64>,
}
struct MonteCarloTrialSampler {
base_inputs: BallisticInputs,
atmosphere: AtmosphericConditions,
solver_max_range: f64,
target_distance: f64,
baseline_at_target: Vector3<f64>,
velocity_delta_dist: rand_distr::Normal<f64>,
angle_dist: rand_distr::Normal<f64>,
bc_dist: rand_distr::Normal<f64>,
wind_sampler: MonteCarloWindSampler,
azimuth_dist: rand_distr::Normal<f64>,
}
impl MonteCarloTrialSampler {
fn new(
base_inputs: BallisticInputs,
base_wind: &WindConditions,
params: &MonteCarloParams,
wind_direction_std_dev: f64,
) -> Result<Self, BallisticsError> {
use rand_distr::Normal;
let atmosphere = AtmosphericConditions {
temperature: base_inputs.temperature,
pressure: base_inputs.pressure,
humidity: base_inputs.humidity_percent(),
altitude: base_inputs.altitude,
};
let target_hint = params
.target_distance
.unwrap_or(base_inputs.target_distance);
let solver_max_range = target_hint.max(1000.0) * 2.0;
let mut baseline_solver =
TrajectorySolver::new(base_inputs.clone(), base_wind.clone(), atmosphere.clone());
baseline_solver.set_max_range(solver_max_range);
let baseline_result = baseline_solver.solve()?;
let target_distance = params.target_distance.unwrap_or(baseline_result.max_range);
let baseline_at_target = baseline_result
.position_at_range(target_distance)
.ok_or("Could not interpolate baseline at target distance")?;
let velocity_delta_dist = Normal::new(0.0, params.velocity_std_dev)
.map_err(|e| format!("Invalid velocity distribution: {}", e))?;
let angle_dist = Normal::new(base_inputs.muzzle_angle, params.angle_std_dev)
.map_err(|e| format!("Invalid angle distribution: {}", e))?;
let bc_dist = Normal::new(base_inputs.bc_value, params.bc_std_dev)
.map_err(|e| format!("Invalid BC distribution: {}", e))?;
let wind_sampler = MonteCarloWindSampler::new(
base_wind,
params.wind_speed_std_dev,
wind_direction_std_dev,
)?;
let azimuth_dist = Normal::new(base_inputs.azimuth_angle, params.azimuth_std_dev)
.map_err(|e| format!("Invalid azimuth distribution: {}", e))?;
Ok(Self {
base_inputs,
atmosphere,
solver_max_range,
target_distance,
baseline_at_target,
velocity_delta_dist,
angle_dist,
bc_dist,
wind_sampler,
azimuth_dist,
})
}
fn sample_one_trial<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> Option<TrialOutcome> {
use rand_distr::Distribution;
let mut inputs = self.base_inputs.clone();
let muzzle_velocity_delta = self.velocity_delta_dist.sample(&mut *rng);
inputs.muzzle_angle = self.angle_dist.sample(&mut *rng);
inputs.bc_value = self.bc_dist.sample(&mut *rng).max(0.01);
inputs.azimuth_angle = self.azimuth_dist.sample(&mut *rng);
let wind = self.wind_sampler.sample(&mut *rng);
let mut solver = TrajectorySolver::new(inputs, wind, self.atmosphere.clone());
solver.inputs.muzzle_velocity =
(solver.inputs.muzzle_velocity + muzzle_velocity_delta).max(0.0);
solver.set_max_range(self.solver_max_range);
let result = solver.solve().ok()?;
let impact_position = if result.max_range < self.target_distance {
Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)
} else {
let pos_at_target = result.position_at_range(self.target_distance)?;
Vector3::new(
0.0,
pos_at_target.y - self.baseline_at_target.y,
pos_at_target.z - self.baseline_at_target.z,
)
};
Some(TrialOutcome {
range: result.max_range,
impact_velocity: result.impact_velocity,
impact_position,
})
}
}
pub fn run_monte_carlo(
base_inputs: BallisticInputs,
params: MonteCarloParams,
) -> Result<MonteCarloResults, BallisticsError> {
run_monte_carlo_with_direction_std_dev(base_inputs, params, 0.0)
}
pub fn run_monte_carlo_with_direction_std_dev(
base_inputs: BallisticInputs,
params: MonteCarloParams,
wind_direction_std_dev: f64,
) -> Result<MonteCarloResults, BallisticsError> {
let base_wind = WindConditions {
speed: params.base_wind_speed,
direction: params.base_wind_direction,
vertical_speed: 0.0,
};
run_monte_carlo_with_wind_and_direction_std_dev(
base_inputs,
base_wind,
params,
wind_direction_std_dev,
)
}
pub fn run_monte_carlo_with_wind(
base_inputs: BallisticInputs,
base_wind: WindConditions,
params: MonteCarloParams,
) -> Result<MonteCarloResults, BallisticsError> {
run_monte_carlo_with_wind_and_direction_std_dev(base_inputs, base_wind, params, 0.0)
}
pub fn run_monte_carlo_with_wind_and_direction_std_dev(
base_inputs: BallisticInputs,
base_wind: WindConditions,
params: MonteCarloParams,
wind_direction_std_dev: f64,
) -> Result<MonteCarloResults, BallisticsError> {
let mut rng = rand::rng();
run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
base_inputs,
base_wind,
params,
wind_direction_std_dev,
&mut rng,
)
}
pub fn run_monte_carlo_with_wind_and_direction_std_dev_seeded(
base_inputs: BallisticInputs,
base_wind: WindConditions,
params: MonteCarloParams,
wind_direction_std_dev: f64,
seed: u64,
) -> Result<MonteCarloResults, BallisticsError> {
use rand::{rngs::StdRng, SeedableRng};
let mut rng = StdRng::seed_from_u64(seed);
run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
base_inputs,
base_wind,
params,
wind_direction_std_dev,
&mut rng,
)
}
fn run_monte_carlo_with_wind_and_direction_std_dev_using_rng<R: rand::Rng + ?Sized>(
base_inputs: BallisticInputs,
base_wind: WindConditions,
params: MonteCarloParams,
wind_direction_std_dev: f64,
rng: &mut R,
) -> Result<MonteCarloResults, BallisticsError> {
let mut ranges = Vec::new();
let mut impact_velocities = Vec::new();
let mut impact_positions = Vec::new();
let sampler = MonteCarloTrialSampler::new(
base_inputs,
&base_wind,
¶ms,
wind_direction_std_dev,
)?;
for _ in 0..params.num_simulations {
if let Some(outcome) = sampler.sample_one_trial(rng) {
ranges.push(outcome.range);
impact_velocities.push(outcome.impact_velocity);
impact_positions.push(outcome.impact_position);
}
}
if ranges.is_empty() {
return Err("No successful simulations".into());
}
Ok(MonteCarloResults {
ranges,
impact_velocities,
impact_positions,
})
}
pub const MC_ADAPTIVE_SCHEMA_VERSION_V1: u32 = 1;
pub const MC_ADAPTIVE_METHOD_V1: &str = "anytime_beta_binomial_mixture_cs_v1";
pub const MC_ADAPTIVE_ASSUMPTIONS_V1: [&str; 4] = [
"Sampling uncertainty only: intervals cover Monte Carlo sampling error, not model error in the trajectory solver or its inputs.",
"Anytime-valid stopping: the beta-binomial mixture confidence sequence keeps its coverage guarantee despite stopping the moment the target half-width is met.",
"Input dispersions are the independent normal distributions declared in MonteCarloParams; correlations between inputs are not modeled.",
"Continuous statistics are streaming Welford moments over trials that reached the target plane, reported with sample (n-1) standard deviations; hit probability's denominator includes all trials.",
];
#[derive(Debug, Clone)]
pub struct McConvergence {
pub level: ConfidenceLevel,
pub target_half_width: f64,
pub min_samples: u64,
pub max_samples: u64,
pub batch_size: u64,
}
impl Default for McConvergence {
fn default() -> Self {
Self {
level: ConfidenceLevel::P95,
target_half_width: 0.02,
min_samples: 1_000,
max_samples: 100_000,
batch_size: 500,
}
}
}
impl McConvergence {
pub fn validate(&self) -> Result<(), String> {
if !self.target_half_width.is_finite() || self.target_half_width <= 0.0 {
return Err(format!(
"McConvergence.target_half_width must be a finite value greater than zero (got {})",
self.target_half_width
));
}
if self.batch_size == 0 {
return Err("McConvergence.batch_size must be greater than zero".to_string());
}
if self.max_samples == 0 {
return Err("McConvergence.max_samples must be greater than zero".to_string());
}
if self.max_samples < self.min_samples {
return Err(format!(
"McConvergence.max_samples ({}) must be at least McConvergence.min_samples ({})",
self.max_samples, self.min_samples
));
}
Ok(())
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize)]
#[serde(rename_all = "snake_case")]
pub enum McStopReason {
TargetHalfWidthMet,
MaxSamplesReached,
}
#[derive(Debug, Clone, serde::Serialize)]
pub struct AdaptiveMcReportV1 {
pub schema_version: u32,
pub method: String,
pub assumptions: Vec<String>,
pub confidence_percent: u32,
pub hit_probability: f64,
pub ci_low: f64,
pub ci_high: f64,
pub samples: u64,
pub attempts: u64,
pub arrivals: u64,
pub stop_reason: McStopReason,
pub hit_radius_m: f64,
pub target_distance_m: f64,
pub mean_impact_velocity_mps: f64,
pub std_impact_velocity_mps: f64,
pub mean_drop_at_target_m: f64,
pub std_drop_at_target_m: f64,
pub mean_wind_drift_at_target_m: f64,
pub std_wind_drift_at_target_m: f64,
}
pub fn run_monte_carlo_adaptive_seeded(
base_inputs: &BallisticInputs,
base_wind: &WindConditions,
params: &MonteCarloParams,
convergence: &McConvergence,
hit_radius_m: f64,
seed: u64,
) -> Result<AdaptiveMcReportV1, String> {
use rand::{rngs::StdRng, SeedableRng};
convergence.validate()?;
let sampler = MonteCarloTrialSampler::new(base_inputs.clone(), base_wind, params, 0.0)
.map_err(|e| e.to_string())?;
let mut rng = StdRng::seed_from_u64(seed);
let mut hits_cs = BernoulliConfidenceSequence::new(convergence.level);
let mut impact_velocity = Welford::new();
let mut drop_at_target = Welford::new();
let mut drift_at_target = Welford::new();
let mut attempts: u64 = 0;
let mut stop_reason = McStopReason::MaxSamplesReached;
while attempts < convergence.max_samples {
let batch = convergence.batch_size.min(convergence.max_samples - attempts);
let mut batch_hits: u64 = 0;
let mut batch_trials: u64 = 0;
for _ in 0..batch {
attempts += 1;
let Some(outcome) = sampler.sample_one_trial(&mut rng) else {
continue; };
batch_trials += 1;
if MonteCarloResults::position_is_hit(&outcome.impact_position, hit_radius_m) {
batch_hits += 1;
}
if MonteCarloResults::position_reached_target(&outcome.impact_position) {
impact_velocity.push(outcome.impact_velocity);
drop_at_target.push(outcome.impact_position.y);
drift_at_target.push(outcome.impact_position.z);
}
}
hits_cs.update_batch(batch_hits, batch_trials);
if hits_cs.trials() >= convergence.min_samples
&& hits_cs.half_width() <= convergence.target_half_width
{
stop_reason = McStopReason::TargetHalfWidthMet;
break;
}
}
let samples = hits_cs.trials();
if samples == 0 {
return Err("No successful simulations".to_string());
}
let (ci_low, ci_high) = hits_cs.bounds();
Ok(AdaptiveMcReportV1 {
schema_version: MC_ADAPTIVE_SCHEMA_VERSION_V1,
method: MC_ADAPTIVE_METHOD_V1.to_string(),
assumptions: MC_ADAPTIVE_ASSUMPTIONS_V1
.iter()
.map(|s| s.to_string())
.collect(),
confidence_percent: convergence.level.as_percent(),
hit_probability: hits_cs.successes() as f64 / samples as f64,
ci_low,
ci_high,
samples,
attempts,
arrivals: drop_at_target.count(),
stop_reason,
hit_radius_m,
target_distance_m: sampler.target_distance,
mean_impact_velocity_mps: impact_velocity.mean(),
std_impact_velocity_mps: impact_velocity.sample_std(),
mean_drop_at_target_m: drop_at_target.mean(),
std_drop_at_target_m: drop_at_target.sample_std(),
mean_wind_drift_at_target_m: drift_at_target.mean(),
std_wind_drift_at_target_m: drift_at_target.sample_std(),
})
}
pub fn calculate_zero_angle(
inputs: BallisticInputs,
target_distance: f64,
target_height: f64,
) -> Result<f64, BallisticsError> {
calculate_zero_angle_with_conditions(
inputs,
target_distance,
target_height,
WindConditions::default(),
AtmosphericConditions::default(),
)
}
pub fn calculate_zero_angle_with_conditions(
inputs: BallisticInputs,
target_distance: f64,
target_height: f64,
wind: WindConditions,
atmosphere: AtmosphericConditions,
) -> Result<f64, BallisticsError> {
let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
}
pub fn calculate_zero_angle_with_resolved_conditions(
inputs: BallisticInputs,
target_distance: f64,
target_height: f64,
wind: WindConditions,
atmosphere: AtmosphericConditions,
) -> Result<f64, BallisticsError> {
let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
}
pub const ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M: f64 = 2000.0;
pub fn calculate_zero_range_from_angle_with_conditions(
inputs: BallisticInputs,
zero_angle_rad: f64,
target_height: f64,
wind: WindConditions,
atmosphere: AtmosphericConditions,
) -> Result<ZeroCrossings, BallisticsError> {
let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
}
pub fn calculate_zero_range_from_angle_with_resolved_conditions(
inputs: BallisticInputs,
zero_angle_rad: f64,
target_height: f64,
wind: WindConditions,
atmosphere: AtmosphericConditions,
) -> Result<ZeroCrossings, BallisticsError> {
let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum BcFitMode {
Drop,
Velocity,
}
#[derive(Debug, Clone, Copy)]
pub struct BcEstimate {
pub bc: f64,
pub rms_error: f64,
pub drag_model: DragModel,
pub mode: BcFitMode,
pub at_bound: bool,
}
fn fit_value_at(
points: &[TrajectoryPoint],
target_dist: f64,
mode: BcFitMode,
drop_offset: f64,
) -> Option<f64> {
let val = |p: &TrajectoryPoint| match mode {
BcFitMode::Drop => drop_offset - p.position.y,
BcFitMode::Velocity => p.velocity_magnitude,
};
for i in 0..points.len() {
if points[i].position.x >= target_dist {
if i == 0 {
return Some(val(&points[0]));
}
let p1 = &points[i - 1];
let p2 = &points[i];
let dx = p2.position.x - p1.position.x;
if dx.abs() < 1e-9 {
return Some(val(p2));
}
let t = (target_dist - p1.position.x) / dx;
return Some(val(p1) + t * (val(p2) - val(p1)));
}
}
None
}
fn fit_residual_sse(
trajectory: &[TrajectoryPoint],
observations: &[(f64, f64)],
mode: BcFitMode,
drop_offset: f64,
) -> Option<f64> {
if observations.is_empty() {
return None;
}
let mut total = 0.0;
for (target_dist, target_val) in observations {
let value = fit_value_at(trajectory, *target_dist, mode, drop_offset)?;
let error = value - target_val;
total += error * error;
}
Some(total)
}
#[allow(clippy::too_many_arguments)] pub fn estimate_bc_fit(
velocity: f64,
mass: f64,
diameter: f64,
points: &[(f64, f64)],
drag_model: DragModel,
mode: BcFitMode,
atmosphere: AtmosphericConditions,
zero_range: Option<f64>,
sight_height: f64,
) -> Result<BcEstimate, BallisticsError> {
if points.is_empty() {
return Err(BallisticsError::from(
"No data points provided for BC estimation.".to_string(),
));
}
let max_dist = points.iter().map(|(d, _)| *d).fold(0.0_f64, f64::max);
let drop_offset = if zero_range.is_some() { sight_height } else { 0.0 };
let sse = |bc_value: f64| -> Option<f64> {
let mut inputs = BallisticInputs {
muzzle_velocity: velocity,
bc_value,
bc_type: drag_model,
bullet_mass: mass,
bullet_diameter: diameter,
sight_height,
..Default::default()
};
if let Some(zr) = zero_range {
let za = calculate_zero_angle_with_conditions(
inputs.clone(),
zr,
sight_height,
WindConditions::default(),
atmosphere.clone(),
)
.ok()?;
inputs.muzzle_angle = za;
}
let mut solver =
TrajectorySolver::new(inputs, WindConditions::default(), atmosphere.clone());
solver.set_max_range(max_dist * 1.5);
let result = solver.solve().ok()?;
fit_residual_sse(&result.points, points, mode, drop_offset)
};
let (bc_min, bc_max) = match drag_model {
DragModel::G7 => (0.05, 0.70),
_ => (0.10, 1.20),
};
let mut best_bc = f64::NAN;
let mut best_sse = f64::MAX;
let mut bc = bc_min;
while bc <= bc_max + 1e-9 {
if let Some(s) = sse(bc) {
if s < best_sse {
best_sse = s;
best_bc = bc;
}
}
bc += 0.01;
}
if !best_bc.is_finite() {
return Err(BallisticsError::from(
"Unable to estimate BC from provided data. Check that the values and units are correct."
.to_string(),
));
}
let lo = (best_bc - 0.01).max(bc_min);
let hi = (best_bc + 0.01).min(bc_max);
let mut bc = lo;
while bc <= hi + 1e-9 {
if let Some(s) = sse(bc) {
if s < best_sse {
best_sse = s;
best_bc = bc;
}
}
bc += 0.001;
}
let at_bound = best_bc <= bc_min + 0.011 || best_bc >= bc_max - 0.011;
let rms_error = (best_sse / points.len() as f64).sqrt();
Ok(BcEstimate {
bc: best_bc,
rms_error,
drag_model,
mode,
at_bound,
})
}
pub fn estimate_bc_from_trajectory(
velocity: f64,
mass: f64,
diameter: f64,
points: &[(f64, f64)], ) -> Result<f64, BallisticsError> {
estimate_bc_fit(
velocity,
mass,
diameter,
points,
DragModel::G1,
BcFitMode::Drop,
AtmosphericConditions::default(),
None,
0.05,
)
.map(|e| e.bc)
}
use rand;
use rand_distr;
#[cfg(test)]
mod mba737_powder_resolution_tests {
use super::*;
#[test]
fn linear_model_cold_powder_subtracts() {
let v = resolve_powder_adjusted_velocity(823.0, 11.1, true, 0.5486, 21.1, None, None);
assert!((v - (823.0 + 0.5486 * (11.1 - 21.1))).abs() < 1e-12);
assert!(v < 823.0);
}
#[test]
fn linear_model_hot_powder_adds() {
let v = resolve_powder_adjusted_velocity(823.0, 31.1, true, 0.5486, 21.1, None, None);
assert!((v - (823.0 + 0.5486 * 10.0)).abs() < 1e-12);
}
#[test]
fn disabled_flag_is_passthrough() {
let v = resolve_powder_adjusted_velocity(823.0, 40.0, false, 0.5486, 21.1, None, None);
assert_eq!(v, 823.0);
}
#[test]
fn curve_overrides_linear_and_interpolates_at_powder_temp() {
let curve = [(4.4, 798.6), (21.1, 823.0), (37.8, 841.2)];
let v = resolve_powder_adjusted_velocity(823.0, 30.0, true, 99.0, 21.1, Some(&curve), Some(4.4));
assert!((v - 798.6).abs() < 1e-9);
}
#[test]
fn curve_falls_back_to_ambient_and_clamps() {
let curve = [(4.4, 798.6), (37.8, 841.2)];
let v = resolve_powder_adjusted_velocity(823.0, -40.0, true, 1.0, 21.1, Some(&curve), None);
assert!((v - 798.6).abs() < 1e-9);
let v_hot = resolve_powder_adjusted_velocity(823.0, 60.0, true, 1.0, 21.1, Some(&curve), None);
assert!((v_hot - 841.2).abs() < 1e-9);
}
#[test]
fn empty_curve_suppresses_linear_fallback() {
let v = resolve_powder_adjusted_velocity(823.0, 40.0, true, 0.5486, 21.1, Some(&[]), None);
assert_eq!(v, 823.0);
}
#[test]
fn sweep_huge_range_errors_instead_of_overflowing() {
assert!(parse_powder_sweep("0:1e20:1").is_err());
assert!(parse_powder_sweep("0:1e308:1e-3").is_err());
}
#[test]
fn sweep_fractional_step_keeps_end_row() {
let rows = parse_powder_sweep("0:0.3:0.1").unwrap();
assert_eq!(rows.len(), 4);
assert!((rows[3] - 0.3).abs() < 1e-9);
}
#[test]
fn solver_and_helper_agree_on_linear_model() {
let inputs = BallisticInputs {
use_powder_sensitivity: true,
powder_temp_sensitivity: 0.5486,
powder_temp: 21.1,
temperature: 4.4,
..Default::default()
};
let expected = resolve_powder_adjusted_velocity(
inputs.muzzle_velocity,
inputs.temperature,
true,
0.5486,
21.1,
None,
None,
);
let solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
assert!((solver.inputs.muzzle_velocity - expected).abs() < 1e-12);
}
}
#[cfg(test)]
mod mba1302_solver_seam_tests {
use super::*;
use crate::wind::WindSegment;
#[test]
fn authoritative_station_atmosphere_preserves_explicit_standard_values_at_altitude() {
let atmosphere = AtmosphericConditions {
temperature: 15.0,
pressure: 1013.25,
humidity: 50.0,
altitude: 2_000.0,
};
let legacy = TrajectorySolver::new(
BallisticInputs::default(),
WindConditions::default(),
atmosphere.clone(),
);
let authoritative = TrajectorySolver::new_with_resolved_station_atmosphere(
BallisticInputs::default(),
WindConditions::default(),
atmosphere,
);
let (legacy_density, _, legacy_temp_c, legacy_pressure_hpa) = legacy.resolved_atmosphere();
let (authoritative_density, _, authoritative_temp_c, authoritative_pressure_hpa) =
authoritative.resolved_atmosphere();
let (icao_temp_k, icao_pressure_pa) =
crate::atmosphere::calculate_icao_standard_atmosphere(2_000.0);
let (expected_authoritative_density, _) =
crate::atmosphere::calculate_atmosphere(2_000.0, Some(15.0), Some(1013.25), 50.0);
assert!((legacy_temp_c - (icao_temp_k - 273.15)).abs() < 1e-12);
assert!((legacy_pressure_hpa - icao_pressure_pa / 100.0).abs() < 1e-12);
assert_eq!(authoritative_temp_c.to_bits(), 15.0_f64.to_bits());
assert_eq!(authoritative_pressure_hpa.to_bits(), 1013.25_f64.to_bits());
assert_eq!(
authoritative_density.to_bits(),
expected_authoritative_density.to_bits()
);
assert!(
(authoritative_density - legacy_density).abs() > 0.1,
"explicit standard values at altitude must differ from ICAO-at-altitude: explicit={authoritative_density}, ICAO={legacy_density}"
);
}
#[test]
fn precomputed_absolute_resolution_via_authoritative_matches_legacy_new() {
for (temperature, pressure, altitude) in [
(15.0, 1013.25, 0.0), (15.0, 1013.25, 2000.0), (-5.0, 850.0, 2000.0), (22.0, 950.0, 500.0),
] {
let atmosphere = AtmosphericConditions {
temperature,
pressure,
humidity: 50.0,
altitude,
};
let legacy = TrajectorySolver::new(
BallisticInputs::default(),
WindConditions::default(),
atmosphere.clone(),
);
let (resolved_temp_c, resolved_pressure_hpa) =
crate::atmosphere::resolve_station_conditions_with_pressure_mode(
temperature,
pressure,
altitude,
crate::atmosphere::PressureReferenceMode::Absolute,
);
let precomputed_atmosphere = AtmosphericConditions {
temperature: resolved_temp_c,
pressure: resolved_pressure_hpa,
humidity: 50.0,
altitude,
};
let precomputed = TrajectorySolver::new_with_resolved_station_atmosphere(
BallisticInputs::default(),
WindConditions::default(),
precomputed_atmosphere,
);
let (legacy_density, legacy_sos, legacy_temp_c, legacy_pressure_hpa) =
legacy.resolved_atmosphere();
let (pre_density, pre_sos, pre_temp_c, pre_pressure_hpa) =
precomputed.resolved_atmosphere();
assert_eq!(
legacy_temp_c.to_bits(),
pre_temp_c.to_bits(),
"temperature=({temperature}, {pressure}, {altitude})"
);
assert_eq!(
legacy_pressure_hpa.to_bits(),
pre_pressure_hpa.to_bits(),
"pressure=({temperature}, {pressure}, {altitude})"
);
assert_eq!(legacy_density.to_bits(), pre_density.to_bits());
assert_eq!(legacy_sos.to_bits(), pre_sos.to_bits());
}
}
fn configured_euler_zero(vertical_wind_mps: f64, time_step_s: f64) -> TrajectorySolver {
let inputs = BallisticInputs {
muzzle_velocity: 800.0,
bc_value: 0.5,
bc_type: DragModel::G7,
bullet_mass: 0.0109,
bullet_diameter: 0.00782,
bullet_length: 0.0309,
sight_height: 0.05,
ground_threshold: -100.0,
use_rk4: false,
use_adaptive_rk45: false,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(300.0);
solver.set_time_step(time_step_s);
if vertical_wind_mps != 0.0 {
solver.set_wind_segments(vec![WindSegment {
speed_kmh: 0.0,
angle_deg: 0.0,
until_m: 400.0,
vertical_mps: vertical_wind_mps,
}]);
}
solver
}
#[test]
fn inclined_shot_zeroes_like_a_level_rifle() {
const ZERO_DISTANCE_M: f64 = 91.44; const SIGHT_HEIGHT_M: f64 = 0.0381;
let inputs = BallisticInputs {
bc_value: 0.5,
bullet_mass: 150.0 * 0.06479891 / 1000.0,
muzzle_velocity: 2700.0 * 0.3048,
sight_height: SIGHT_HEIGHT_M,
..Default::default()
};
let mut level = inputs.clone();
level.shooting_angle = 0.0;
let level_angle = TrajectorySolver::new(level, Default::default(), Default::default())
.find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
.expect("level zero must solve");
let mut inclined = inputs;
inclined.shooting_angle = 5.71_f64.to_radians();
let inclined_angle =
TrajectorySolver::new(inclined, Default::default(), Default::default())
.find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
.expect("MBA-1412: a 5.71 deg incline at a 100 yd zero must be solvable");
assert!(
(inclined_angle - level_angle).abs() < 1e-9,
"zeroing is level-rifle sight geometry; incline must not move the solved zero: \
level={level_angle}, inclined={inclined_angle}"
);
}
#[test]
fn configured_zero_keeps_segments_method_and_time_step_then_sets_base_angle() {
const TARGET_DISTANCE_M: f64 = 150.0;
const TARGET_HEIGHT_M: f64 = 0.05;
let mut segmented = configured_euler_zero(-10.0, 0.02);
let coarse_height = segmented
.zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
.expect("coarse configured trial")
.expect("coarse trial reaches target");
let mut fine = segmented.clone();
fine.set_time_step(0.001);
let fine_height = fine
.zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
.expect("fine configured trial")
.expect("fine trial reaches target");
assert!(
(coarse_height - fine_height).abs() > 1e-5,
"configured Euler step must affect zero trials: coarse={coarse_height}, fine={fine_height}"
);
let segmented_angle = segmented
.calculate_and_set_zero_angle(TARGET_DISTANCE_M, TARGET_HEIGHT_M, ZeroTargetFrame::SightLine)
.expect("segmented zero");
assert_eq!(
segmented.inputs.muzzle_angle.to_bits(),
segmented_angle.to_bits(),
"successful zero must install its angle on the configured solver"
);
assert_eq!(segmented.time_step.to_bits(), 0.02_f64.to_bits());
assert_eq!(segmented.max_range.to_bits(), 300.0_f64.to_bits());
assert!(segmented.wind_sock.is_some());
assert_eq!(
segmented.station_atmosphere_resolution,
StationAtmosphereResolution::Authoritative
);
let zero_height = segmented
.zero_trial_height_at(segmented_angle, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
.expect("verify segmented zero")
.expect("zeroed trial reaches target");
assert!(
(zero_height - TARGET_HEIGHT_M).abs() < 0.0001,
"configured zero missed target: height={zero_height}"
);
let mut calm = configured_euler_zero(0.0, 0.02);
let calm_angle = calm
.calculate_and_set_zero_angle(TARGET_DISTANCE_M, TARGET_HEIGHT_M, ZeroTargetFrame::SightLine)
.expect("calm zero");
assert!(
(segmented_angle - calm_angle).abs() > 1e-5,
"segmented vertical wind must participate in zero trials: segmented={segmented_angle}, calm={calm_angle}"
);
}
}
#[cfg(test)]
mod result_sanity_tests {
use super::*;
fn default_solver() -> TrajectorySolver {
TrajectorySolver::new(
BallisticInputs::default(),
WindConditions::default(),
AtmosphericConditions::default(),
)
}
fn minimal_result() -> TrajectoryResult {
TrajectoryResult {
max_range: 100.0,
max_height: 1.0,
time_of_flight: 0.5,
impact_velocity: 700.0,
impact_energy: 2450.0,
projectile_mass_kg: 0.01,
line_of_sight_height_m: 1.5,
station_speed_of_sound_mps: 340.0,
termination: TrajectoryTermination::MaxRange,
points: vec![],
sampled_points: None,
min_pitch_damping: None,
transonic_mach: None,
angular_state: None,
max_yaw_angle: None,
max_precession_angle: None,
aerodynamic_jump: None,
mach_1_2_distance_m: None,
mach_1_0_distance_m: None,
mach_0_9_distance_m: None,
}
}
#[test]
fn mba1293_negative_scalars_fail_the_result_postcondition() {
let solver = default_solver();
solver
.validate_result_sanity(&minimal_result())
.expect("a sane result must pass");
for (name, mutate) in [
("max_range", (|r| r.max_range = -50.588) as fn(&mut TrajectoryResult)),
("time_of_flight", |r| r.time_of_flight = -1.0),
("impact_velocity", |r| r.impact_velocity = -700.0),
("impact_energy", |r| r.impact_energy = -1.0),
] {
let mut result = minimal_result();
mutate(&mut result);
let error = solver
.validate_result_sanity(&result)
.expect_err("negative scalar must fail");
assert!(
error.to_string().contains(name),
"error for {name} did not name the field: {error}"
);
}
}
#[test]
fn mba1293_speed_budget_bounds_legitimate_states_and_rejects_divergence() {
let solver = default_solver();
let mv = solver.inputs.muzzle_velocity;
let position = Vector3::new(10.0, 0.0, 0.0);
solver
.validate_integration_state(&position, &Vector3::new(mv, 0.0, 0.0), 0.01)
.expect("muzzle-speed state must pass");
let error = solver
.validate_integration_state(&position, &Vector3::new(-13.0 * mv, 0.0, 0.0), 0.01)
.expect_err("13x muzzle speed must fail the budget");
assert!(error.to_string().contains("diverged"), "{error}");
let after_fall = mv + crate::constants::G_ACCEL_MPS2 * 60.0;
solver
.validate_integration_state(&position, &Vector3::new(0.0, -after_fall, 0.0), 60.0)
.expect("gravity-accelerated speed within g*t must pass");
}
}
#[cfg(test)]
mod trajectory_point_budget_tests {
use super::*;
use crate::MAX_TRAJECTORY_SAMPLES;
fn solver_with_budget(
use_rk4: bool,
use_adaptive_rk45: bool,
point_budget: usize,
max_range: f64,
) -> TrajectorySolver {
let inputs = BallisticInputs {
use_rk4,
use_adaptive_rk45,
ground_threshold: f64::NEG_INFINITY,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.max_trajectory_points = point_budget;
solver.set_max_range(max_range);
solver.set_time_step(0.001);
solver
}
#[test]
fn mba1283_every_solver_errors_instead_of_exceeding_point_budget() {
for (mode, use_rk4, use_adaptive_rk45) in [
("Euler", false, false),
("RK4", true, false),
("RK45", true, true),
] {
let error = solver_with_budget(use_rk4, use_adaptive_rk45, 3, 10.0)
.solve()
.expect_err("a solve requiring more than three points must fail");
assert!(
error.to_string().contains("point limit of 3"),
"unexpected {mode} point-budget error: {error}"
);
}
}
#[test]
fn mba1283_interpolated_endpoint_counts_toward_point_budget() {
for (mode, use_rk4, use_adaptive_rk45) in [
("Euler", false, false),
("RK4", true, false),
("RK45", true, true),
] {
let result = solver_with_budget(use_rk4, use_adaptive_rk45, 2, 0.1)
.solve()
.expect("the initial point plus exact endpoint fit a two-point budget");
assert_eq!(result.points.len(), 2, "unexpected {mode} point count");
let error = solver_with_budget(use_rk4, use_adaptive_rk45, 1, 0.1)
.solve()
.expect_err("the exact endpoint must not exceed a one-point budget");
assert!(
error.to_string().contains("point limit of 1"),
"unexpected {mode} endpoint-budget error: {error}"
);
}
}
#[test]
fn mba1299_every_solver_preflights_the_sample_budget() {
for (mode, use_rk4, use_adaptive_rk45) in [
("Euler", false, false),
("RK4", true, false),
("RK45", true, true),
] {
let inputs = BallisticInputs {
use_rk4,
use_adaptive_rk45,
enable_trajectory_sampling: true,
sample_interval: 1.0,
ground_threshold: f64::NEG_INFINITY,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(MAX_TRAJECTORY_SAMPLES as f64);
solver.max_trajectory_points = 0;
let error = solver
.solve()
.expect_err("an over-limit sample grid must fail before integration");
assert!(
error
.to_string()
.contains("trajectory sample limit of 250000 exceeded"),
"unexpected {mode} sample-budget error: {error}"
);
}
}
#[test]
fn mba1299_normal_sampling_does_not_change_solver_results() {
for (mode, use_rk4, use_adaptive_rk45) in [
("Euler", false, false),
("RK4", true, false),
("RK45", true, true),
] {
let solve = |enable_trajectory_sampling| {
let inputs = BallisticInputs {
use_rk4,
use_adaptive_rk45,
enable_trajectory_sampling,
sample_interval: 0.5,
ground_threshold: f64::NEG_INFINITY,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(2.0);
solver.solve().expect("normal short-range solve")
};
let baseline = solve(false);
let sampled = solve(true);
for (field, left, right) in [
("max_range", baseline.max_range, sampled.max_range),
("max_height", baseline.max_height, sampled.max_height),
(
"time_of_flight",
baseline.time_of_flight,
sampled.time_of_flight,
),
(
"impact_velocity",
baseline.impact_velocity,
sampled.impact_velocity,
),
(
"impact_energy",
baseline.impact_energy,
sampled.impact_energy,
),
] {
assert_eq!(
left.to_bits(),
right.to_bits(),
"{mode} sampling changed {field}"
);
}
assert_eq!(baseline.points.len(), sampled.points.len());
for (index, (left, right)) in baseline
.points
.iter()
.zip(&sampled.points)
.enumerate()
{
assert_eq!(left.time.to_bits(), right.time.to_bits(), "{mode} point {index}");
assert_eq!(
left.position.map(f64::to_bits),
right.position.map(f64::to_bits),
"{mode} point {index} position"
);
assert_eq!(
left.velocity_magnitude.to_bits(),
right.velocity_magnitude.to_bits(),
"{mode} point {index} velocity"
);
assert_eq!(
left.kinetic_energy.to_bits(),
right.kinetic_energy.to_bits(),
"{mode} point {index} energy"
);
}
assert!(baseline.sampled_points.is_none());
let samples = sampled
.sampled_points
.expect("sampling-enabled solve should return observations");
assert_eq!(
samples
.iter()
.map(|sample| sample.distance_m)
.collect::<Vec<_>>(),
vec![0.0, 0.5, 1.0, 1.5, 2.0],
"{mode} normal sampling grid changed"
);
}
}
}
#[cfg(test)]
mod monte_carlo_result_tests {
use super::*;
fn make_results(impact_positions: Vec<Vector3<f64>>) -> MonteCarloResults {
let count = impact_positions.len();
MonteCarloResults {
ranges: vec![500.0; count],
impact_velocities: vec![300.0; count],
impact_positions,
}
}
#[test]
fn target_plane_cep_excludes_shortfall_markers() {
let mut positions: Vec<Vector3<f64>> = (1..=5)
.map(|radius| Vector3::new(0.0, radius as f64, 0.0))
.collect();
positions.extend(
(0..5).map(|_| Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)),
);
let results = make_results(positions);
assert_eq!(results.target_arrival_count(), 5);
assert_eq!(results.target_shortfall_fraction(), 0.5);
assert_eq!(results.target_plane_cep(), Some(3.0));
let one_shortfall = make_results(vec![
Vector3::new(0.0, 1.0, 0.0),
Vector3::new(0.0, 2.0, 0.0),
Vector3::new(0.0, 3.0, 0.0),
Vector3::new(0.0, 4.0, 0.0),
Vector3::new(0.0, 5.0, 0.0),
Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
]);
assert_eq!(one_shortfall.target_plane_cep(), Some(3.0));
}
#[test]
fn all_shortfalls_have_no_cep_but_still_count_as_misses() {
let all_shortfalls = make_results(vec![
Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
]);
assert_eq!(all_shortfalls.target_arrival_count(), 0);
assert_eq!(all_shortfalls.target_shortfall_fraction(), 1.0);
assert_eq!(all_shortfalls.target_plane_cep(), None);
assert_eq!(all_shortfalls.hit_probability(0.3), 0.0);
let one_hit_one_shortfall = make_results(vec![
Vector3::new(0.0, 0.1, 0.0),
Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
]);
assert_eq!(one_hit_one_shortfall.hit_probability(0.3), 0.5);
}
#[test]
fn rect_hit_probability_checks_independent_axis_halves() {
let results = make_results(vec![
Vector3::new(0.0, 0.1, 0.1),
Vector3::new(0.0, 0.0, 0.2),
Vector3::new(0.0, 0.0, 0.201),
Vector3::new(0.0, 0.301, 0.0),
Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
]);
assert!((results.rect_hit_probability(0.4, 0.6) - 0.4).abs() < 1e-12);
}
#[test]
fn rect_hit_probability_matches_circular_hit_probability_for_a_centered_hit() {
let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
assert_eq!(results.rect_hit_probability(0.5, 0.5), 1.0);
assert_eq!(results.hit_probability(0.3), 1.0);
}
#[test]
fn rect_hit_probability_is_zero_for_empty_or_nonpositive_dimensions() {
let empty = make_results(vec![]);
assert_eq!(empty.rect_hit_probability(1.0, 1.0), 0.0);
let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
assert_eq!(results.rect_hit_probability(0.0, 1.0), 0.0);
assert_eq!(results.rect_hit_probability(1.0, 0.0), 0.0);
assert_eq!(results.rect_hit_probability(-1.0, 1.0), 0.0);
}
}
#[cfg(test)]
mod monte_carlo_seeded_tests {
use super::*;
fn seeded_test_fixture() -> (BallisticInputs, WindConditions) {
(
BallisticInputs {
muzzle_velocity: 800.0,
..BallisticInputs::default()
},
WindConditions::default(),
)
}
fn loose_params() -> MonteCarloParams {
MonteCarloParams {
num_simulations: 1, velocity_std_dev: 3.0,
angle_std_dev: 3.5e-4,
bc_std_dev: 0.01,
wind_speed_std_dev: 1.0,
target_distance: Some(300.0),
base_wind_speed: 0.0,
base_wind_direction: 0.0,
azimuth_std_dev: 3.5e-4,
}
}
fn mixed_arrival_params() -> MonteCarloParams {
MonteCarloParams {
num_simulations: 1,
target_distance: Some(1920.0),
..MonteCarloParams::default()
}
}
#[test]
fn legacy_seeded_estimates_are_pinned_bit_for_bit() {
let (inputs, wind) = seeded_test_fixture();
let params = MonteCarloParams {
num_simulations: 200,
target_distance: Some(500.0),
..MonteCarloParams::default()
};
let results = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
inputs,
wind,
params,
0.01,
0x1352_5EED,
)
.expect("seeded legacy run");
assert_eq!(results.ranges.len(), 200, "ranges length");
assert_eq!(results.impact_velocities.len(), 200, "impact_velocities length");
assert_eq!(results.impact_positions.len(), 200, "impact_positions length");
assert_eq!(
results.hit_probability(DEFAULT_HIT_RADIUS_M).to_bits(),
0.14_f64.to_bits(),
"hit_probability = {:?}",
results.hit_probability(DEFAULT_HIT_RADIUS_M)
);
let expected_ranges: [f64; 3] =
[1907.972891143359, 1936.408435469319, 1912.8150447617645];
let expected_velocities: [f64; 3] =
[238.6187151542299, 239.91923651600106, 243.14112455427164];
let expected_positions: [(f64, f64, f64); 3] = [
(0.0, -0.0643556039548101, 0.7344970252014579),
(0.0, 0.5769422971539162, 0.27227201756386726),
(0.0, -0.7440425792472842, 0.1541804446282822),
];
for (i, expected) in expected_ranges.iter().enumerate() {
assert_eq!(
results.ranges[i].to_bits(),
expected.to_bits(),
"ranges[{i}] = {:?}, pinned {expected:?}",
results.ranges[i]
);
}
for (i, expected) in expected_velocities.iter().enumerate() {
assert_eq!(
results.impact_velocities[i].to_bits(),
expected.to_bits(),
"impact_velocities[{i}] = {:?}, pinned {expected:?}",
results.impact_velocities[i]
);
}
for (i, (x, y, z)) in expected_positions.iter().enumerate() {
let actual = results.impact_positions[i];
assert_eq!(actual.x.to_bits(), x.to_bits(), "impact_positions[{i}].x = {:?}", actual.x);
assert_eq!(actual.y.to_bits(), y.to_bits(), "impact_positions[{i}].y = {:?}", actual.y);
assert_eq!(actual.z.to_bits(), z.to_bits(), "impact_positions[{i}].z = {:?}", actual.z);
}
}
#[test]
fn seeded_runs_are_deterministic_and_match_the_using_rng_path() {
let inputs = BallisticInputs {
muzzle_velocity: 800.0,
..BallisticInputs::default()
};
let params = MonteCarloParams {
num_simulations: 64,
target_distance: Some(200.0),
..MonteCarloParams::default()
};
let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
inputs.clone(),
WindConditions::default(),
params.clone(),
0.01,
42,
)
.expect("seeded run a");
let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
inputs,
WindConditions::default(),
params,
0.01,
42,
)
.expect("seeded run b");
assert_eq!(a.ranges.len(), b.ranges.len());
for (ra, rb) in a.ranges.iter().zip(b.ranges.iter()) {
assert_eq!(ra.to_bits(), rb.to_bits());
}
for (pa, pb) in a.impact_positions.iter().zip(b.impact_positions.iter()) {
assert_eq!(pa.x.to_bits(), pb.x.to_bits());
assert_eq!(pa.y.to_bits(), pb.y.to_bits());
assert_eq!(pa.z.to_bits(), pb.z.to_bits());
}
}
#[test]
fn different_seeds_generally_produce_different_draws() {
let inputs = BallisticInputs {
muzzle_velocity: 800.0,
..BallisticInputs::default()
};
let params = MonteCarloParams {
num_simulations: 32,
velocity_std_dev: 5.0,
target_distance: Some(200.0),
..MonteCarloParams::default()
};
let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
inputs.clone(),
WindConditions::default(),
params.clone(),
0.0,
1,
)
.expect("seeded run a");
let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
inputs,
WindConditions::default(),
params,
0.0,
2,
)
.expect("seeded run b");
assert_ne!(a.impact_velocities, b.impact_velocities);
}
#[test]
fn adaptive_stops_at_target_half_width_on_an_easy_case() {
let (inputs, wind) = seeded_test_fixture();
let conv = McConvergence {
target_half_width: 0.05,
..Default::default()
};
let r = run_monte_carlo_adaptive_seeded(
&inputs,
&wind,
&loose_params(),
&conv,
DEFAULT_HIT_RADIUS_M,
0x1352_ADA9,
)
.unwrap();
assert_eq!(r.stop_reason, McStopReason::TargetHalfWidthMet);
assert!(
(r.ci_high - r.ci_low) / 2.0 <= 0.05 + 1e-12,
"half-width {} exceeds the requested 0.05",
(r.ci_high - r.ci_low) / 2.0
);
assert!(r.samples >= conv.min_samples, "stopped below min_samples");
assert!(r.samples < conv.max_samples, "did not actually stop early");
assert!(r.samples.is_multiple_of(conv.batch_size) || r.samples == conv.min_samples);
assert!(r.ci_low <= r.hit_probability && r.hit_probability <= r.ci_high);
assert_eq!(r.hit_radius_m, DEFAULT_HIT_RADIUS_M);
assert_eq!(r.target_distance_m, 300.0);
assert_eq!(r.confidence_percent, 95);
assert!(
r.samples > 1,
"params.num_simulations must be ignored by the adaptive driver"
);
assert!(
r.mean_impact_velocity_mps > 0.0,
"no impact velocity accumulated"
);
assert!(
r.std_drop_at_target_m > 0.0 && r.std_wind_drift_at_target_m > 0.0,
"dispersion collapsed: drop sd {} drift sd {}",
r.std_drop_at_target_m,
r.std_wind_drift_at_target_m
);
}
#[test]
fn adaptive_caps_at_max_samples_on_an_impossible_target() {
let (inputs, wind) = seeded_test_fixture();
let conv = McConvergence {
target_half_width: 1e-6,
max_samples: 3_000,
batch_size: 500,
min_samples: 1_000,
level: ConfidenceLevel::P95,
};
let r = run_monte_carlo_adaptive_seeded(
&inputs,
&wind,
&loose_params(),
&conv,
DEFAULT_HIT_RADIUS_M,
7,
)
.unwrap();
assert_eq!(r.stop_reason, McStopReason::MaxSamplesReached);
assert_eq!(r.samples, 3_000);
}
#[test]
fn adaptive_stops_between_the_floor_and_the_ceiling() {
let (inputs, wind) = seeded_test_fixture();
let conv = McConvergence {
level: ConfidenceLevel::P95,
target_half_width: 0.03,
min_samples: 0,
max_samples: 10_000,
batch_size: 100,
};
let r = run_monte_carlo_adaptive_seeded(
&inputs,
&wind,
&loose_params(),
&conv,
DEFAULT_HIT_RADIUS_M,
0x1352_5A1D,
)
.unwrap();
assert_eq!(r.stop_reason, McStopReason::TargetHalfWidthMet);
assert!(r.samples > 0);
assert!(
r.samples.is_multiple_of(conv.batch_size),
"samples {} is not a whole number of batches",
r.samples
);
assert!(
r.samples > conv.min_samples,
"stopped on the floor, not on the data"
);
assert!(
r.samples < conv.max_samples,
"ran to the ceiling, so nothing adaptive was exercised"
);
assert!(
r.samples > conv.batch_size,
"stopped on the very first look ({} samples); the multi-batch path is untested",
r.samples
);
assert!((r.ci_high - r.ci_low) / 2.0 <= 0.03 + 1e-12);
assert!(r.ci_low <= r.hit_probability && r.hit_probability <= r.ci_high);
assert_eq!(r.attempts, r.samples, "no trial should have been dropped");
}
#[test]
fn adaptive_runs_a_truncated_final_batch_up_to_max_samples() {
let (inputs, wind) = seeded_test_fixture();
let conv = McConvergence {
level: ConfidenceLevel::P95,
target_half_width: 1e-6,
min_samples: 0,
max_samples: 750,
batch_size: 500,
};
let r = run_monte_carlo_adaptive_seeded(
&inputs,
&wind,
&loose_params(),
&conv,
DEFAULT_HIT_RADIUS_M,
0x1352_7B10,
)
.unwrap();
assert_eq!(r.stop_reason, McStopReason::MaxSamplesReached);
assert_eq!(
r.samples, 750,
"the 250-trial final batch did not run, or was not truncated"
);
assert_eq!(r.attempts, 750);
assert!(!r.samples.is_multiple_of(conv.batch_size));
}
#[test]
fn adaptive_is_deterministic_for_a_seed() {
let (inputs, wind) = seeded_test_fixture();
let conv = McConvergence::default();
let a = run_monte_carlo_adaptive_seeded(
&inputs,
&wind,
&loose_params(),
&conv,
DEFAULT_HIT_RADIUS_M,
99,
)
.unwrap();
let b = run_monte_carlo_adaptive_seeded(
&inputs,
&wind,
&loose_params(),
&conv,
DEFAULT_HIT_RADIUS_M,
99,
)
.unwrap();
assert_eq!(a.hit_probability.to_bits(), b.hit_probability.to_bits());
assert_eq!(a.samples, b.samples);
assert_eq!(a.ci_low.to_bits(), b.ci_low.to_bits());
assert_eq!(a.ci_high.to_bits(), b.ci_high.to_bits());
assert_eq!(
a.mean_impact_velocity_mps.to_bits(),
b.mean_impact_velocity_mps.to_bits()
);
assert_eq!(
a.std_drop_at_target_m.to_bits(),
b.std_drop_at_target_m.to_bits()
);
}
#[test]
fn adaptive_report_carries_schema_method_and_all_four_assumptions() {
let (inputs, wind) = seeded_test_fixture();
let conv = McConvergence {
min_samples: 0,
max_samples: 50,
batch_size: 50,
target_half_width: 1.0,
level: ConfidenceLevel::P90,
};
let r = run_monte_carlo_adaptive_seeded(
&inputs,
&wind,
&mixed_arrival_params(),
&conv,
DEFAULT_HIT_RADIUS_M,
0x1352_D0C5,
)
.unwrap();
assert_eq!(r.schema_version, MC_ADAPTIVE_SCHEMA_VERSION_V1);
assert_eq!(r.schema_version, 1);
assert_eq!(r.method, "anytime_beta_binomial_mixture_cs_v1");
assert_eq!(r.confidence_percent, 90);
assert_eq!(r.attempts, 50, "one full batch was drawn");
assert_eq!(r.samples, 50, "no trial was dropped by the solver");
assert!(
r.arrivals > 0 && r.arrivals < r.samples,
"fixture must split the run: arrivals {} of samples {}",
r.arrivals,
r.samples
);
assert!(r.arrivals >= 2, "arrivals {} too few for a sample sd", r.arrivals);
assert!(r.std_drop_at_target_m > 0.0 && r.std_impact_velocity_mps > 0.0);
assert!(r.attempts >= r.samples && r.samples >= r.arrivals);
assert_eq!(r.assumptions.len(), 4, "exactly four assumptions expected");
assert_eq!(
r.assumptions[0],
"Sampling uncertainty only: intervals cover Monte Carlo sampling error, not model error in the trajectory solver or its inputs."
);
assert_eq!(
r.assumptions[1],
"Anytime-valid stopping: the beta-binomial mixture confidence sequence keeps its coverage guarantee despite stopping the moment the target half-width is met."
);
assert_eq!(
r.assumptions[2],
"Input dispersions are the independent normal distributions declared in MonteCarloParams; correlations between inputs are not modeled."
);
assert_eq!(
r.assumptions[3],
"Continuous statistics are streaming Welford moments over trials that reached the target plane, reported with sample (n-1) standard deviations; hit probability's denominator includes all trials."
);
assert_eq!(
serde_json::to_string(&McStopReason::TargetHalfWidthMet).unwrap(),
"\"target_half_width_met\""
);
assert_eq!(
serde_json::to_string(&McStopReason::MaxSamplesReached).unwrap(),
"\"max_samples_reached\""
);
}
#[test]
fn adaptive_rejects_nonsense_convergence() {
let (inputs, wind) = seeded_test_fixture();
let run = |conv: McConvergence| {
run_monte_carlo_adaptive_seeded(
&inputs,
&wind,
&loose_params(),
&conv,
DEFAULT_HIT_RADIUS_M,
1,
)
.unwrap_err()
};
for bad_width in [0.0, -0.01, f64::NAN] {
let err = run(McConvergence {
target_half_width: bad_width,
..Default::default()
});
assert!(
err.contains("target_half_width"),
"error must name the field, got: {err}"
);
}
let err = run(McConvergence {
batch_size: 0,
..Default::default()
});
assert!(err.contains("batch_size"), "got: {err}");
let err = run(McConvergence {
min_samples: 5_000,
max_samples: 1_000,
..Default::default()
});
assert!(err.contains("max_samples"), "got: {err}");
assert!(err.contains("min_samples"), "got: {err}");
let err = run(McConvergence {
max_samples: 0,
min_samples: 0,
..Default::default()
});
assert!(err.contains("max_samples"), "got: {err}");
let err = McConvergence {
batch_size: 0,
..Default::default()
}
.validate()
.unwrap_err();
assert!(err.contains("batch_size"));
}
#[test]
fn wilson_companion_matches_hit_probability_and_wilson_interval() {
let (inputs, wind) = seeded_test_fixture();
let params = MonteCarloParams {
num_simulations: 128,
target_distance: Some(500.0),
..MonteCarloParams::default()
};
let results = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
inputs,
wind,
params,
0.01,
0x1352_C0DE,
)
.expect("seeded legacy run");
for level in [
ConfidenceLevel::P90,
ConfidenceLevel::P95,
ConfidenceLevel::P99,
] {
let (p_hat, (lo, hi), n) =
results.hit_probability_wilson(DEFAULT_HIT_RADIUS_M, level);
assert_eq!(
p_hat.to_bits(),
results.hit_probability(DEFAULT_HIT_RADIUS_M).to_bits()
);
assert_eq!(n, results.impact_positions.len() as u64);
let hits = results
.impact_positions
.iter()
.filter(|p| MonteCarloResults::position_is_hit(p, DEFAULT_HIT_RADIUS_M))
.count() as u64;
let (want_lo, want_hi) = wilson_interval(hits, n, level);
assert_eq!(lo.to_bits(), want_lo.to_bits());
assert_eq!(hi.to_bits(), want_hi.to_bits());
assert!(lo <= p_hat && p_hat <= hi, "interval excludes p_hat");
}
let empty = MonteCarloResults {
ranges: Vec::new(),
impact_velocities: Vec::new(),
impact_positions: Vec::new(),
};
assert_eq!(
empty.hit_probability_wilson(DEFAULT_HIT_RADIUS_M, ConfidenceLevel::P95),
(0.0, (0.0, 1.0), 0)
);
}
}
#[cfg(test)]
mod monte_carlo_powder_curve_tests {
use super::*;
use rand::{rngs::StdRng, SeedableRng};
#[test]
fn powder_curve_preserves_sampled_muzzle_velocity_dispersion() {
let inputs = BallisticInputs {
muzzle_velocity: 700.0,
powder_temp_curve: Some(vec![(15.0, 800.0)]),
powder_curve_temp_c: Some(15.0),
..BallisticInputs::default()
};
let params = MonteCarloParams {
num_simulations: 16,
velocity_std_dev: 20.0,
angle_std_dev: 1e-12,
bc_std_dev: 1e-12,
wind_speed_std_dev: 1e-12,
target_distance: Some(100.0),
azimuth_std_dev: 1e-12,
..MonteCarloParams::default()
};
let mut rng = StdRng::seed_from_u64(0x5EED_1176);
let results = run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
inputs,
WindConditions::default(),
params,
0.0,
&mut rng,
)
.expect("Monte Carlo solve");
let min_velocity = results
.impact_velocities
.iter()
.copied()
.fold(f64::INFINITY, f64::min);
let max_velocity = results
.impact_velocities
.iter()
.copied()
.fold(f64::NEG_INFINITY, f64::max);
assert!(
max_velocity - min_velocity > 1.0,
"20 m/s muzzle spread collapsed after curve resolution: impact-velocity span={} m/s",
max_velocity - min_velocity
);
}
}
#[cfg(test)]
mod monte_carlo_wind_sampling_tests {
use super::*;
use rand::{rngs::StdRng, SeedableRng};
#[test]
fn wind_speed_sigma_does_not_change_seeded_direction_draws() {
let base_wind = WindConditions {
speed: 100.0,
direction: 0.37,
vertical_speed: 0.0,
};
let narrow_speed = MonteCarloWindSampler::new(&base_wind, 0.5, 0.2).unwrap();
let wide_speed = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
let mut speed_changed = false;
for _ in 0..32 {
let narrow = narrow_speed.sample(&mut narrow_rng);
let wide = wide_speed.sample(&mut wide_rng);
assert!(narrow.speed > 0.0 && wide.speed > 0.0);
assert_eq!(narrow.direction.to_bits(), wide.direction.to_bits());
speed_changed |= narrow.speed.to_bits() != wide.speed.to_bits();
}
assert!(
speed_changed,
"different speed sigmas must still vary speed draws"
);
}
#[test]
fn zero_direction_sigma_has_no_angular_jitter() {
let base_wind = WindConditions {
speed: 100.0,
direction: 0.37,
vertical_speed: 0.0,
};
let sampler = MonteCarloWindSampler::new(&base_wind, 4.0, 0.0).unwrap();
let mut rng = StdRng::seed_from_u64(0x5EED_1223);
let mut speed_changed = false;
for _ in 0..32 {
let wind = sampler.sample(&mut rng);
speed_changed |= wind.speed.to_bits() != base_wind.speed.to_bits();
assert_eq!(wind.direction.to_bits(), base_wind.direction.to_bits());
}
assert!(speed_changed, "speed uncertainty should remain active");
}
#[test]
fn direction_sigma_controls_seeded_angular_spread_in_radians() {
let base_wind = WindConditions {
speed: 100.0,
direction: 0.37,
vertical_speed: 0.0,
};
let narrow = MonteCarloWindSampler::new(&base_wind, 4.0, 0.1).unwrap();
let wide = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
let mut nonzero_direction_draw = false;
for _ in 0..32 {
let narrow_wind = narrow.sample(&mut narrow_rng);
let wide_wind = wide.sample(&mut wide_rng);
assert_eq!(narrow_wind.speed.to_bits(), wide_wind.speed.to_bits());
let narrow_delta = narrow_wind.direction - base_wind.direction;
let wide_delta = wide_wind.direction - base_wind.direction;
assert!((wide_delta - 2.0 * narrow_delta).abs() < 1e-12);
nonzero_direction_draw |= narrow_delta.abs() > 1e-6;
}
assert!(
nonzero_direction_draw,
"positive radians sigma must vary direction"
);
}
#[test]
fn direction_sigma_rejects_negative_or_nonfinite_values() {
let base_wind = WindConditions::default();
for sigma in [-0.1, f64::NAN, f64::INFINITY] {
assert!(MonteCarloWindSampler::new(&base_wind, 1.0, sigma).is_err());
}
}
#[test]
fn base_vertical_wind_rides_into_every_mc_sample() {
use rand::SeedableRng;
let base_wind = WindConditions { vertical_speed: 4.2, ..Default::default() };
let sampler = MonteCarloWindSampler::new(&base_wind, 1.0, 0.2).unwrap();
let mut rng = rand::rngs::StdRng::seed_from_u64(7);
for _ in 0..32 {
let w = sampler.sample(&mut rng);
assert_eq!(w.vertical_speed, 4.2);
}
}
#[test]
fn negative_speed_sample_reverses_wind_direction() {
let direction = 0.25;
let signed_speed = -2.5;
let wind = wind_from_signed_speed_sample(signed_speed, direction, 0.0);
let positive_wind = wind_from_signed_speed_sample(2.5, direction, 0.0);
assert_eq!(wind.speed, 2.5);
assert!(
(wind.direction - (direction + std::f64::consts::PI)).abs() < f64::EPSILON,
"negative speed must reverse direction by pi: got {}",
wind.direction
);
assert_eq!(positive_wind.speed, 2.5);
assert_eq!(positive_wind.direction, direction);
let normalized_x = -wind.speed * wind.direction.cos();
let normalized_z = -wind.speed * wind.direction.sin();
let signed_x = -signed_speed * direction.cos();
let signed_z = -signed_speed * direction.sin();
assert!((normalized_x - signed_x).abs() < 1e-12);
assert!((normalized_z - signed_z).abs() < 1e-12);
}
}
#[cfg(test)]
mod bc_fit_objective_tests {
use super::*;
fn velocity_point(range_m: f64, velocity_mps: f64) -> TrajectoryPoint {
TrajectoryPoint {
time: 0.0,
position: Vector3::new(range_m, 0.0, 0.0),
velocity_magnitude: velocity_mps,
kinetic_energy: 0.0,
drag_coefficient: None,
}
}
#[test]
fn candidate_that_misses_an_observation_has_no_score() {
let trajectory = vec![velocity_point(0.0, 800.0), velocity_point(100.0, 700.0)];
let observations = vec![(50.0, 750.0), (150.0, 600.0)];
assert!(
fit_residual_sse(&trajectory, &observations, BcFitMode::Velocity, 0.0).is_none(),
"a candidate that reaches only one of two observations must not compete on partial SSE"
);
let complete_observations = vec![(50.0, 740.0), (100.0, 680.0)];
assert_eq!(
fit_residual_sse(
&trajectory,
&complete_observations,
BcFitMode::Velocity,
0.0,
),
Some(500.0)
);
}
}
#[cfg(test)]
mod cluster_bc_reference_space_tests {
use super::*;
fn acceleration_at_1100_fps(inputs: BallisticInputs) -> Vector3<f64> {
let solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
let position = Vector3::zeros();
let velocity = Vector3::new(1100.0 / 3.28084, 0.0, 0.0);
let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
solver.calculate_acceleration(
&position,
&velocity,
&Vector3::zeros(),
(temp_c, pressure_hpa, density / 1.225),
)
}
#[test]
fn solver_passes_g7_reference_model_to_cluster_classifier() {
let inputs = BallisticInputs {
bc_value: 0.190,
bc_type: DragModel::G7,
bullet_mass: 77.0 * crate::constants::GRAINS_TO_KG,
bullet_diameter: 0.224 * 0.0254,
use_cluster_bc: true,
..BallisticInputs::default()
};
let solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
let corrected = solver.apply_cluster_bc_correction(0.190, 2800.0);
assert!(
(corrected / 0.190 - 1.004).abs() < 1e-12,
"solver selected the wrong G7 cluster multiplier: {}",
corrected / 0.190
);
}
#[test]
fn velocity_bc_segments_are_not_cluster_corrected_twice() {
let segmented_clustered = BallisticInputs {
bc_value: 0.5,
bc_type: DragModel::G7,
use_bc_segments: true,
bc_segments_data: Some(vec![
crate::BCSegmentData {
velocity_min: 0.0,
velocity_max: 1_600.0,
bc_value: 0.4,
},
crate::BCSegmentData {
velocity_min: 1_600.0,
velocity_max: 5_000.0,
bc_value: 0.45,
},
]),
use_cluster_bc: true,
..BallisticInputs::default()
};
let mut segmented_only = segmented_clustered.clone();
segmented_only.use_cluster_bc = false;
let mut constant_clustered = segmented_clustered.clone();
constant_clustered.bc_value = 0.4;
constant_clustered.bc_segments_data = None;
let stacked = acceleration_at_1100_fps(segmented_clustered);
let segment_only = acceleration_at_1100_fps(segmented_only);
let cluster_only = acceleration_at_1100_fps(constant_clustered);
assert!(
(stacked.x - segment_only.x).abs() < 1e-12,
"segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
stacked.x,
segment_only.x
);
assert!(
(cluster_only.x - segment_only.x).abs() > 1e-6,
"cluster correction must remain active for a constant BC"
);
}
#[test]
fn mach_bc_segments_are_not_cluster_corrected_twice() {
let mach_segmented_clustered = BallisticInputs {
bc_value: 0.5,
bc_type: DragModel::G7,
use_bc_segments: false,
bc_segments: Some(vec![(0.5, 0.3), (1.5, 0.5)]),
use_cluster_bc: true,
..BallisticInputs::default()
};
let mut mach_segmented_only = mach_segmented_clustered.clone();
mach_segmented_only.use_cluster_bc = false;
let stacked = acceleration_at_1100_fps(mach_segmented_clustered);
let segment_only = acceleration_at_1100_fps(mach_segmented_only);
assert!(
(stacked.x - segment_only.x).abs() < 1e-12,
"Mach segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
stacked.x,
segment_only.x
);
}
}
#[cfg(test)]
mod velocity_bc_flag_tests {
use super::*;
fn acceleration_at_600_mps(inputs: BallisticInputs) -> Vector3<f64> {
let solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
solver.calculate_acceleration(
&Vector3::zeros(),
&Vector3::new(600.0, 0.0, 0.0),
&Vector3::zeros(),
(temp_c, pressure_hpa, density / 1.225),
)
}
#[test]
fn velocity_bc_data_requires_opt_in_in_trajectory_solver() {
let scalar_inputs = BallisticInputs {
bc_value: 0.5,
bc_type: DragModel::G7,
..BallisticInputs::default()
};
let mut disabled_inputs = scalar_inputs.clone();
disabled_inputs.bc_segments_data = Some(vec![crate::BCSegmentData {
velocity_min: 0.0,
velocity_max: 4_000.0,
bc_value: 0.46,
}]);
disabled_inputs.use_bc_segments = false;
let mut enabled_inputs = disabled_inputs.clone();
enabled_inputs.use_bc_segments = true;
let mut mach_only_inputs = scalar_inputs.clone();
mach_only_inputs.bc_segments = Some(vec![(0.0, 0.4), (3.0, 0.4)]);
let mut disabled_with_both = mach_only_inputs.clone();
disabled_with_both.bc_segments_data = disabled_inputs.bc_segments_data.clone();
let scalar = acceleration_at_600_mps(scalar_inputs);
let disabled = acceleration_at_600_mps(disabled_inputs);
let enabled = acceleration_at_600_mps(enabled_inputs);
let mach_only = acceleration_at_600_mps(mach_only_inputs);
let disabled_with_both = acceleration_at_600_mps(disabled_with_both);
assert_eq!(
disabled.x.to_bits(),
scalar.x.to_bits(),
"a populated velocity table must not change drag while use_bc_segments is false"
);
assert!(
enabled.x < disabled.x - 1.0,
"enabling the lower BC table must increase drag: disabled ax={} enabled ax={}",
disabled.x,
enabled.x
);
assert_eq!(
disabled_with_both.x.to_bits(),
mach_only.x.to_bits(),
"disabling velocity data must fall through to an explicit Mach table"
);
}
}
#[cfg(test)]
mod mach_bc_segment_tests {
use super::*;
#[test]
fn trajectory_solver_interpolates_explicit_mach_bc_segments() {
let segmented_inputs = BallisticInputs {
bc_value: 0.8,
use_bc_segments: false,
bc_segments: Some(vec![(1.0, 0.2), (2.0, 0.4)]),
bc_segments_data: None,
..BallisticInputs::default()
};
let mut expected_inputs = segmented_inputs.clone();
expected_inputs.bc_value = 0.3;
expected_inputs.bc_segments = None;
let atmosphere = AtmosphericConditions::default();
let segmented_solver = TrajectorySolver::new(
segmented_inputs,
WindConditions::default(),
atmosphere.clone(),
);
let expected_solver = TrajectorySolver::new(
expected_inputs,
WindConditions::default(),
atmosphere,
);
let position = Vector3::zeros();
let (density, _, temp_c, pressure_hpa) = segmented_solver.resolved_atmosphere();
let (_, local_speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
segmented_solver.atmosphere.altitude,
segmented_solver.atmosphere.altitude,
temp_c,
pressure_hpa,
density / 1.225,
segmented_solver.atmosphere.humidity,
);
let velocity = Vector3::new(1.5 * local_speed_of_sound, 0.0, 0.0);
let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
let segmented_acceleration = segmented_solver.calculate_acceleration(
&position,
&velocity,
&Vector3::zeros(),
resolved_atmo,
);
let expected_acceleration = expected_solver.calculate_acceleration(
&position,
&velocity,
&Vector3::zeros(),
resolved_atmo,
);
assert!(
(segmented_acceleration.x - expected_acceleration.x).abs() < 1e-12,
"Mach 1.5 must interpolate BC 0.3: segmented ax={} expected ax={}",
segmented_acceleration.x,
expected_acceleration.x
);
}
}
#[cfg(test)]
mod custom_drag_table_validation_tests {
use super::*;
#[test]
fn solve_accepts_zero_bc_when_custom_table_present() {
let inputs = BallisticInputs {
bc_value: 0.0, bullet_mass: 0.0106,
bullet_diameter: 0.00782,
muzzle_velocity: 850.0,
custom_drag_table: Some(crate::drag::DragTable::new(
vec![0.5, 1.0, 2.0, 3.0],
vec![0.23, 0.40, 0.30, 0.26],
)),
..BallisticInputs::default()
};
let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
assert!(solver.solve().is_ok());
}
#[test]
fn solve_still_requires_bc_without_table() {
let inputs = BallisticInputs {
bc_value: 0.0,
bullet_mass: 0.0106,
bullet_diameter: 0.00782,
muzzle_velocity: 850.0,
..BallisticInputs::default()
};
let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
assert!(solver.solve().is_err());
}
}
#[cfg(test)]
mod cd_scale_tests {
use super::*;
fn deck() -> crate::drag::DragTable {
crate::drag::DragTable::new(vec![0.5, 1.0, 2.0, 3.0], vec![0.23, 0.40, 0.30, 0.26])
}
fn deck_inputs(cd_scale: f64) -> BallisticInputs {
BallisticInputs {
bullet_mass: 0.0106,
bullet_diameter: 0.00782,
muzzle_velocity: 850.0,
custom_drag_table: Some(deck()),
cd_scale,
..BallisticInputs::default()
}
}
#[test]
fn default_cd_scale_is_one() {
assert_eq!(BallisticInputs::default().cd_scale, 1.0);
}
#[test]
fn cd_scale_absent_is_byte_identical_to_explicit_one() {
let omitted = BallisticInputs {
bullet_mass: 0.0106,
bullet_diameter: 0.00782,
muzzle_velocity: 850.0,
custom_drag_table: Some(deck()),
..BallisticInputs::default()
};
let explicit = BallisticInputs {
cd_scale: 1.0,
..omitted.clone()
};
let solver_omitted =
TrajectorySolver::new(omitted, WindConditions::default(), AtmosphericConditions::default());
let solver_explicit =
TrajectorySolver::new(explicit, WindConditions::default(), AtmosphericConditions::default());
let cd_omitted = solver_omitted.calculate_drag_coefficient(700.0, 340.0);
let cd_explicit = solver_explicit.calculate_drag_coefficient(700.0, 340.0);
assert_eq!(
cd_omitted.to_bits(),
cd_explicit.to_bits(),
"default cd_scale must be bit-identical to an explicit 1.0"
);
let result = solver_omitted.solve();
assert!(result.is_ok(), "existing custom-deck solves must pass unchanged");
}
#[test]
fn cd_scale_multiplies_the_interpolated_cd_exactly() {
let velocity = 700.0;
let speed_of_sound = 340.0;
let mach = velocity / speed_of_sound;
let expected_unscaled = deck().interpolate(mach);
for &scale in &[0.90, 1.0, 1.10, 1.5] {
let solver = TrajectorySolver::new(
deck_inputs(scale),
WindConditions::default(),
AtmosphericConditions::default(),
);
let cd = solver.calculate_drag_coefficient(velocity, speed_of_sound);
assert!(
(cd - expected_unscaled * scale).abs() < 1e-12,
"scale={scale}: cd={cd} expected={}",
expected_unscaled * scale
);
}
}
#[test]
fn cd_scale_direction_on_cli_api_solver() {
let solve = |scale: f64| {
TrajectorySolver::new(
deck_inputs(scale),
WindConditions::default(),
AtmosphericConditions::default(),
)
.solve()
.expect("custom-deck solve should succeed")
};
let baseline = solve(1.0);
let scaled_up = solve(1.10);
let scaled_down = solve(0.90);
assert!(
scaled_up.impact_velocity < baseline.impact_velocity,
"cd_scale=1.10 must increase drag -> lower impact velocity: base={} up={}",
baseline.impact_velocity,
scaled_up.impact_velocity
);
assert!(
scaled_down.impact_velocity > baseline.impact_velocity,
"cd_scale=0.90 must decrease drag -> higher impact velocity: base={} down={}",
baseline.impact_velocity,
scaled_down.impact_velocity
);
}
#[test]
fn validate_for_solve_rejects_invalid_cd_scale() {
for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
let solver = TrajectorySolver::new(
deck_inputs(bad),
WindConditions::default(),
AtmosphericConditions::default(),
);
assert!(
solver.solve().is_err(),
"cd_scale={bad} must be rejected by validate_for_solve"
);
}
}
#[test]
fn validate_for_solve_rejects_invalid_cd_scale_without_a_custom_drag_table() {
for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
let inputs = BallisticInputs {
bc_value: 0.5,
bc_type: crate::DragModel::G1,
bullet_mass: 0.0106,
bullet_diameter: 0.00782,
muzzle_velocity: 850.0,
cd_scale: bad,
..BallisticInputs::default()
};
assert!(inputs.custom_drag_table.is_none(), "precondition: no custom deck");
let solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
assert!(
solver.solve().is_err(),
"cd_scale={bad} must be rejected by validate_for_solve even without a custom \
drag table"
);
}
}
#[test]
fn cd_scale_is_inert_without_a_custom_drag_table() {
let make = |cd_scale: f64| BallisticInputs {
bc_value: 0.5,
bc_type: crate::DragModel::G1,
bullet_mass: 0.0106,
bullet_diameter: 0.00782,
muzzle_velocity: 850.0,
cd_scale,
..BallisticInputs::default()
};
let solver_neutral = TrajectorySolver::new(
make(1.0),
WindConditions::default(),
AtmosphericConditions::default(),
);
let solver_far = TrajectorySolver::new(
make(1.5),
WindConditions::default(),
AtmosphericConditions::default(),
);
let cd_neutral = solver_neutral.calculate_drag_coefficient(700.0, 340.0);
let cd_far = solver_far.calculate_drag_coefficient(700.0, 340.0);
assert_eq!(
cd_neutral.to_bits(),
cd_far.to_bits(),
"cd_scale must not affect the G-model/BC drag path"
);
}
#[test]
fn cd_scale_shifts_all_three_solver_paths_in_the_same_direction() {
let cli_solve = |scale: f64| {
TrajectorySolver::new(
deck_inputs(scale),
WindConditions::default(),
AtmosphericConditions::default(),
)
.solve()
.expect("cli_api custom-deck solve should succeed")
};
let cli_baseline = cli_solve(1.0);
let cli_scaled = cli_solve(1.10);
assert!(
cli_scaled.impact_velocity < cli_baseline.impact_velocity,
"cli_api: cd_scale=1.10 must lower impact velocity"
);
let derivatives_accel_x = |scale: f64| {
let inputs = deck_inputs(scale);
crate::derivatives::compute_derivatives(
nalgebra::Vector3::zeros(),
nalgebra::Vector3::new(700.0, 0.0, 0.0),
&inputs,
nalgebra::Vector3::zeros(),
(1.225, 340.0, 0.0, 0.0),
inputs.bc_value,
None,
0.0,
None,
)[3]
};
let deriv_baseline = derivatives_accel_x(1.0);
let deriv_scaled = derivatives_accel_x(1.10);
assert!(
deriv_scaled < deriv_baseline,
"derivatives: cd_scale=1.10 must make x-acceleration more negative (more drag): \
base={deriv_baseline} scaled={deriv_scaled}"
);
let fast_final_speed = |scale: f64| {
let inputs = deck_inputs(scale);
let wind_sock = crate::wind::WindSock::new(vec![]);
let params = crate::fast_trajectory::FastIntegrationParams {
horiz: 500.0,
vert: 0.0,
initial_state: [0.0, 0.0, 0.0, 850.0, 0.0, 0.0],
t_span: (0.0, 5.0),
atmo_params: (0.0, 15.0, 1013.25, 1.0),
atmo_sock: None,
};
let solution = crate::fast_trajectory::fast_integrate(&inputs, &wind_sock, params);
assert!(solution.success, "fast_integrate must succeed for scale={scale}");
let last = solution.t.len() - 1;
let (vx, vy, vz) = (
solution.y[3][last],
solution.y[4][last],
solution.y[5][last],
);
(vx * vx + vy * vy + vz * vz).sqrt()
};
let fast_baseline = fast_final_speed(1.0);
let fast_scaled = fast_final_speed(1.10);
assert!(
fast_scaled < fast_baseline,
"fast_trajectory: cd_scale=1.10 must lower final speed: base={fast_baseline} scaled={fast_scaled}"
);
}
}
#[cfg(test)]
mod humid_local_mach_tests {
use super::*;
fn solver_with_station_humidity(humidity_percent: f64) -> TrajectorySolver {
let inputs = BallisticInputs {
custom_drag_table: Some(crate::drag::DragTable::new(vec![0.5, 1.5], vec![0.1, 1.1])),
..BallisticInputs::default()
};
TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions {
temperature: 30.0,
pressure: 1013.25,
humidity: humidity_percent,
altitude: 0.0,
},
)
}
fn acceleration(solver: &TrajectorySolver, base_ratio: f64) -> Vector3<f64> {
solver.calculate_acceleration(
&Vector3::zeros(),
&Vector3::new(350.0, 0.0, 0.0),
&Vector3::zeros(),
(30.0, 1013.25, base_ratio),
)
}
#[test]
fn local_mach_uses_station_humidity_when_density_is_held_constant() {
let dry = acceleration(&solver_with_station_humidity(0.0), 1.0);
let humid = acceleration(&solver_with_station_humidity(100.0), 1.0);
assert!(
humid.x > dry.x,
"humid sound speed should lower Mach and drag on the rising test curve: dry ax={} humid ax={}",
dry.x,
humid.x
);
}
#[test]
fn active_atmosphere_zone_uses_zone_humidity_instead_of_station_humidity() {
let zone_humidity = 80.0;
let zone_ratio =
crate::atmosphere::calculate_air_density_cimp(30.0, 1013.25, zone_humidity) / 1.225;
let station_solver = solver_with_station_humidity(zone_humidity);
let mut zoned_solver = solver_with_station_humidity(0.0);
zoned_solver.set_atmo_segments(vec![(30.0, 1013.25, zone_humidity, 1_000.0)]);
let station = acceleration(&station_solver, zone_ratio);
let zoned = acceleration(&zoned_solver, zone_ratio);
assert!(
(zoned - station).norm() < 1e-12,
"active zone T/P/RH should override the station atmosphere: station={station:?} zoned={zoned:?}"
);
}
}
#[cfg(test)]
mod inclined_atmosphere_frame_tests {
use super::*;
fn expected_shot_frame_vector(level: Vector3<f64>, angle: f64) -> Vector3<f64> {
let (sin_angle, cos_angle) = angle.sin_cos();
Vector3::new(
level.x * cos_angle + level.y * sin_angle,
-level.x * sin_angle + level.y * cos_angle,
level.z,
)
}
#[test]
fn inclined_positions_at_same_world_altitude_have_same_solver_acceleration() {
let angle = std::f64::consts::FRAC_PI_6;
let inputs = BallisticInputs {
shooting_angle: angle,
..BallisticInputs::default()
};
let atmosphere = AtmosphericConditions {
altitude: 100.0,
..AtmosphericConditions::default()
};
let solver = TrajectorySolver::new(inputs, WindConditions::default(), atmosphere);
let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
let velocity = Vector3::new(600.0, 0.0, 0.0);
let along_slant = Vector3::new(1_000.0, 0.0, 0.0);
let across_slant = Vector3::new(0.0, 500.0 / angle.cos(), 0.0);
let a = solver.calculate_acceleration(
&along_slant,
&velocity,
&Vector3::zeros(),
resolved_atmo,
);
let b = solver.calculate_acceleration(
&across_slant,
&velocity,
&Vector3::zeros(),
resolved_atmo,
);
assert!(
(a - b).norm() < 1e-10,
"solver acceleration differs at equal world altitude: {a:?} vs {b:?}"
);
}
#[test]
fn inclined_headwind_is_rotated_into_solver_frame() {
let angle = std::f64::consts::FRAC_PI_6;
let inputs = BallisticInputs {
shooting_angle: angle,
..BallisticInputs::default()
};
let solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
let level_headwind = Vector3::new(-100.0, 0.0, 0.0);
let velocity = expected_shot_frame_vector(level_headwind, angle);
let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
let actual = solver.calculate_acceleration(
&Vector3::zeros(),
&velocity,
&level_headwind,
(temp_c, pressure_hpa, density / 1.225),
);
assert!(
(actual - solver.gravity_acceleration()).norm() < 1e-12,
"co-moving horizontal wind must leave only shot-frame gravity: {actual:?}"
);
}
#[test]
fn inclined_coriolis_is_rotated_into_solver_frame() {
let angle = std::f64::consts::FRAC_PI_6;
let latitude_deg = 45.0_f64;
let shot_azimuth = 0.4_f64;
let velocity = Vector3::new(600.0, 20.0, 5.0);
let base_inputs = BallisticInputs {
shooting_angle: angle,
latitude: Some(latitude_deg),
shot_azimuth,
..BallisticInputs::default()
};
let acceleration = |enable_coriolis| {
let mut inputs = base_inputs.clone();
inputs.enable_coriolis = enable_coriolis;
let solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
solver.calculate_acceleration(
&Vector3::zeros(),
&velocity,
&Vector3::zeros(),
(temp_c, pressure_hpa, density / 1.225),
)
};
let omega_earth = 7.2921159e-5_f64;
let latitude = latitude_deg.to_radians();
let level_omega = Vector3::new(
omega_earth * latitude.cos() * shot_azimuth.cos(),
omega_earth * latitude.sin(),
-omega_earth * latitude.cos() * shot_azimuth.sin(),
);
let expected = -2.0 * expected_shot_frame_vector(level_omega, angle).cross(&velocity);
let actual = acceleration(true) - acceleration(false);
assert!(
(actual - expected).norm() < 1e-12,
"inclined Coriolis mismatch: actual={actual:?}, expected={expected:?}"
);
}
}
#[cfg(test)]
mod terminal_range_interpolation_tests {
use super::*;
#[test]
fn terminal_finalizer_selects_the_earliest_crossed_boundary() {
let inputs = BallisticInputs {
ground_threshold: 0.0,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(120.0);
let previous_speed = 700.0;
let mut points = vec![TrajectoryPoint {
time: 99.0,
position: Vector3::new(90.0, 1.0, -1.0),
velocity_magnitude: previous_speed,
kinetic_energy: 0.5 * solver.inputs.bullet_mass * previous_speed.powi(2),
drag_coefficient: None,
}];
let mut max_height = 1.0;
let termination = solver
.append_terminal_endpoint(
&mut points,
Vector3::new(130.0, -3.0, 3.0),
Vector3::new(600.0, 0.0, 0.0),
101.0,
&mut max_height,
)
.expect("the final step brackets supported boundaries");
assert_eq!(termination, TrajectoryTermination::GroundThreshold);
assert_eq!(points.len(), 2);
let terminal = points.last().expect("terminal point");
assert_eq!(terminal.time, 99.5);
assert_eq!(terminal.position, Vector3::new(100.0, 0.0, 0.0));
assert_eq!(terminal.velocity_magnitude, 675.0);
assert_eq!(
terminal.kinetic_energy,
0.5 * solver.inputs.bullet_mass * 675.0_f64.powi(2)
);
solver.set_max_range(100.0);
let mut tied_points = vec![points[0].clone()];
assert_eq!(
solver
.append_terminal_endpoint(
&mut tied_points,
Vector3::new(130.0, -3.0, 3.0),
Vector3::new(600.0, 0.0, 0.0),
101.0,
&mut max_height,
)
.expect("tied boundaries remain a valid terminal"),
TrajectoryTermination::GroundThreshold
);
}
#[test]
fn sub_ulp_terminal_crossing_replaces_instead_of_duplicating_range() {
let ground_threshold = f64::from_bits(1.0_f64.to_bits() - 1);
let inputs = BallisticInputs {
ground_threshold,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(1_000.0);
let speed = 700.0;
let mut points = vec![TrajectoryPoint {
time: 0.0,
position: Vector3::new(100.0, 1.0, 0.0),
velocity_magnitude: speed,
kinetic_energy: 0.5 * solver.inputs.bullet_mass * speed.powi(2),
drag_coefficient: None,
}];
let mut max_height = 1.0;
let termination = solver
.append_terminal_endpoint(
&mut points,
Vector3::new(101.0, 0.0, 0.0),
Vector3::new(699.0, 0.0, 0.0),
1.0,
&mut max_height,
)
.expect("sub-ULP ground crossing remains representable as one terminal state");
assert_eq!(termination, TrajectoryTermination::GroundThreshold);
assert_eq!(points.len(), 1);
assert_eq!(points[0].position.x, 100.0);
assert_eq!(points[0].position.y.to_bits(), ground_threshold.to_bits());
assert!(points[0].time > 0.0);
}
#[test]
fn every_solver_appends_an_exact_max_range_endpoint() {
let target_range = 0.1;
let modes = [
("Euler", false, false),
("RK4", true, false),
("RK45", true, true),
];
for (name, use_rk4, use_adaptive_rk45) in modes {
let inputs = BallisticInputs {
use_rk4,
use_adaptive_rk45,
ground_threshold: f64::NEG_INFINITY,
enable_trajectory_sampling: true,
sample_interval: target_range,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(target_range);
let result = solver.solve().expect("short-range solve should succeed");
let terminal = result.points.last().expect("terminal point is missing");
let muzzle = result.points.first().expect("muzzle point is missing");
assert_eq!(result.termination, TrajectoryTermination::MaxRange);
assert_eq!(
terminal.position.x.to_bits(),
target_range.to_bits(),
"{name} did not terminate exactly at max_range"
);
assert_eq!(result.max_range.to_bits(), target_range.to_bits());
assert!(
result.time_of_flight > 0.0 && result.time_of_flight < solver.time_step,
"{name} terminal time was not interpolated within the crossing step: {}",
result.time_of_flight
);
assert_eq!(result.time_of_flight.to_bits(), terminal.time.to_bits());
assert_eq!(
result.impact_velocity.to_bits(),
terminal.velocity_magnitude.to_bits()
);
assert_eq!(
result.impact_energy.to_bits(),
terminal.kinetic_energy.to_bits()
);
let expected_energy = 0.5 * solver.inputs.bullet_mass * result.impact_velocity.powi(2);
assert!((result.impact_energy - expected_energy).abs() < 1e-12);
assert!(terminal.velocity_magnitude < muzzle.velocity_magnitude);
assert!(terminal.kinetic_energy < muzzle.kinetic_energy);
let terminal_sample = result
.sampled_points
.as_ref()
.and_then(|samples| samples.last())
.expect("terminal trajectory sample is missing");
assert_eq!(
terminal_sample.distance_m.to_bits(),
target_range.to_bits(),
"{name} sampling did not include max_range"
);
assert_eq!(
terminal_sample.time_s.to_bits(),
result.time_of_flight.to_bits()
);
assert_eq!(
terminal_sample.velocity_mps.to_bits(),
result.impact_velocity.to_bits()
);
assert!((terminal_sample.energy_j - result.impact_energy).abs() < 1e-12);
}
}
}
#[cfg(test)]
mod precession_inertia_wiring_tests {
use super::*;
#[test]
fn solver_uses_projectile_specific_moments_of_inertia() {
let mass_kg = 55.0 * crate::constants::GRAINS_TO_KG;
let caliber_m = 0.224 * 0.0254;
let length_m = 0.75 * 0.0254;
let inputs = BallisticInputs {
bullet_mass: mass_kg,
bullet_diameter: caliber_m,
bullet_length: length_m,
muzzle_velocity: 800.0,
twist_rate: 7.0,
enable_precession_nutation: true,
use_rk4: false,
use_adaptive_rk45: false,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(0.1);
let (air_density, speed_of_sound, _, _) = solver.resolved_atmosphere();
let velocity_mps = solver.inputs.muzzle_velocity;
let velocity_fps = velocity_mps * 3.28084;
let twist_rate_ft = solver.inputs.twist_rate / 12.0;
let spin_rate_rad_s = (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI;
let initial_state = AngularState {
pitch_angle: 0.001,
yaw_angle: 0.001,
pitch_rate: 0.0,
yaw_rate: 0.0,
precession_angle: 0.0,
nutation_phase: 0.0,
};
let params = PrecessionNutationParams {
mass_kg,
caliber_m,
length_m,
spin_rate_rad_s,
spin_inertia: crate::spin_decay::calculate_moment_of_inertia(
mass_kg, caliber_m, length_m, "ogive",
),
transverse_inertia: crate::pitch_damping::calculate_transverse_moment_of_inertia(
mass_kg, caliber_m, length_m, "ogive",
),
velocity_mps,
air_density_kg_m3: air_density,
mach: velocity_mps / speed_of_sound,
pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
nutation_damping_factor: 0.05,
};
let expected = calculate_combined_angular_motion(
¶ms,
&initial_state,
0.0,
solver.time_step,
0.001,
);
let actual = solver
.solve()
.expect("one-step solve should succeed")
.angular_state
.expect("precession/nutation was enabled");
assert!(
(actual.precession_angle - expected.precession_angle).abs() < 1e-15,
"precession phase used the wrong inertia: actual={}, expected={}",
actual.precession_angle,
expected.precession_angle
);
assert!(
(actual.nutation_phase - expected.nutation_phase).abs() < 1e-15,
"nutation phase used the wrong inertia: actual={}, expected={}",
actual.nutation_phase,
expected.nutation_phase
);
}
}
#[cfg(test)]
mod form_factor_drag_tests {
use super::*;
fn acceleration_with_form_factor_flag(enabled: bool) -> Vector3<f64> {
let inputs = BallisticInputs {
bc_value: 0.462,
bc_type: DragModel::G1,
bullet_model: Some("168gr SMK Match".to_string()),
use_form_factor: enabled,
..BallisticInputs::default()
};
let solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
solver.calculate_acceleration(
&Vector3::zeros(),
&Vector3::new(600.0, 0.0, 0.0),
&Vector3::zeros(),
(temp_c, pressure_hpa, density / 1.225),
)
}
#[test]
fn measured_bc_drag_does_not_apply_name_based_form_factor_again() {
let baseline = acceleration_with_form_factor_flag(false);
let flagged = acceleration_with_form_factor_flag(true);
assert!(
(flagged - baseline).norm() < 1e-12,
"published BC already encodes form factor: baseline={baseline:?} flagged={flagged:?}"
);
}
}
#[cfg(test)]
mod rk45_adaptivity_tests {
use super::*;
#[test]
fn cli_rk45_error_norm_scales_components_independently() {
let position = Vector3::new(1.0e9, 0.0, 0.0);
let velocity = Vector3::new(800.0, 0.0, 0.0);
let fifth_position = position;
let fifth_velocity = velocity;
let fourth_position = position;
let fourth_velocity = Vector3::new(800.0, 1.0e-3, 0.0);
let error = cli_rk45_error_norm(
&position,
&velocity,
&fifth_position,
&fifth_velocity,
&fourth_position,
&fourth_velocity,
);
let expected = 1.0e-3 / 6.0_f64.sqrt();
assert!(
(error - expected).abs() <= 1e-15,
"large downrange position masked a velocity-component error: {error}"
);
}
fn discontinuous_wind_solver() -> TrajectorySolver {
let inputs = BallisticInputs::default();
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_wind_segments(vec![
crate::wind::WindSegment::new(0.0, 90.0, 4.0),
crate::wind::WindSegment::new(1_000.0, 90.0, 10_000.0),
]);
solver
}
#[test]
fn rk45_retries_discontinuous_trial_before_advancing() {
let solver = discontinuous_wind_solver();
let position = Vector3::new(0.0, solver.inputs.muzzle_height, 0.0);
let velocity = Vector3::new(solver.inputs.muzzle_velocity, 0.0, 0.0);
let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
let dt = 0.01;
let rejected_trial = solver.rk45_step(
&position,
&velocity,
dt,
&Vector3::zeros(),
RK45_TOLERANCE,
resolved_atmo,
);
assert!(
rejected_trial.error > RK45_TOLERANCE,
"discontinuous full step must exceed tolerance, got {}",
rejected_trial.error
);
let accepted = solver.adaptive_rk45_step(
&position,
&velocity,
dt,
&Vector3::zeros(),
resolved_atmo,
);
assert!(accepted.used_dt < dt, "oversized trial was not retried");
assert!(
accepted.error <= RK45_TOLERANCE || accepted.used_dt <= RK45_MIN_DT,
"accepted error {} exceeds tolerance at dt {}",
accepted.error,
accepted.used_dt
);
let accepted_trial = solver.rk45_step(
&position,
&velocity,
accepted.used_dt,
&Vector3::zeros(),
RK45_TOLERANCE,
resolved_atmo,
);
assert_eq!(accepted.position, accepted_trial.position);
assert_eq!(accepted.velocity, accepted_trial.velocity);
assert!((RK45_MIN_DT..=RK45_MAX_DT).contains(&accepted.next_dt));
}
}
#[cfg(test)]
mod ground_termination_tests {
use super::*;
use crate::trajectory_observation::TrajectoryObservationFlag;
#[test]
fn every_solver_reports_one_exact_early_ground_endpoint() {
for (name, use_rk4, use_adaptive_rk45) in [
("Euler", false, false),
("RK4", true, false),
("RK45", true, true),
] {
let inputs = BallisticInputs {
muzzle_height: 1.0,
muzzle_angle: -0.2,
ground_threshold: 0.0,
use_rk4,
use_adaptive_rk45,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(1_000.0);
let result = solver.solve().expect("early-ground solve should succeed");
let terminal = result.points.last().expect("terminal point is missing");
assert_eq!(result.termination, TrajectoryTermination::GroundThreshold);
assert_eq!(terminal.position.y.to_bits(), 0.0_f64.to_bits());
assert!(
terminal.position.x < 1_000.0,
"{name} incorrectly reached max range"
);
assert_eq!(result.max_range.to_bits(), terminal.position.x.to_bits());
assert_eq!(
result
.points
.iter()
.filter(|point| point.position.y == 0.0)
.count(),
1,
"{name} did not retain exactly one ground endpoint"
);
let observations = result
.sample_observations(1.0, 100)
.expect("checked early-ground sampling should succeed");
assert!(observations[..observations.len() - 1]
.iter()
.all(|observation| observation.distance_m < terminal.position.x));
let terminal_observation = observations.last().expect("terminal observation");
assert_eq!(
terminal_observation.distance_m.to_bits(),
terminal.position.x.to_bits()
);
assert!(terminal_observation
.flags
.contains(&TrajectoryObservationFlag::Terminal));
assert!(terminal_observation
.flags
.contains(&TrajectoryObservationFlag::GroundThreshold));
assert_eq!(
observations
.iter()
.filter(|observation| observation
.flags
.contains(&TrajectoryObservationFlag::Terminal))
.count(),
1,
"{name} repeated the terminal observation"
);
}
}
#[test]
fn rk4_and_rk45_descend_to_ground_threshold() {
for adaptive in [false, true] {
let inputs = BallisticInputs {
muzzle_angle: 0.1, use_rk4: true,
use_adaptive_rk45: adaptive,
..BallisticInputs::default()
};
assert_eq!(
inputs.ground_threshold, -100.0,
"default ground_threshold is -100 m"
);
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(1.0e7);
let result = solver.solve().expect("solve should succeed");
let final_y = result
.points
.last()
.expect("trajectory has points")
.position
.y;
assert!(
final_y < -1.0,
"adaptive_rk45={adaptive}: final y = {final_y} m; a lofted shot should descend \
past launch level toward the ground_threshold floor, not stop at y = 0"
);
}
}
}
#[cfg(test)]
mod magnus_stability_tests {
use super::*;
#[test]
fn yaw_of_repose_magnus_force_is_vertical_and_twist_signed() {
let acceleration = |enable_magnus, is_twist_right| {
let inputs = BallisticInputs {
muzzle_velocity: 822.96,
bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
bullet_diameter: 0.308 * 0.0254,
bullet_length: 1.215 * 0.0254,
twist_rate: 10.0,
is_twist_right,
enable_magnus,
..BallisticInputs::default()
};
let solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
solver.calculate_acceleration(
&Vector3::zeros(),
&Vector3::new(822.96, 0.0, 0.0),
&Vector3::zeros(),
(temp_c, pressure_hpa, density / 1.225),
)
};
let baseline = acceleration(false, true);
let right_twist = acceleration(true, true) - baseline;
let left_twist = acceleration(true, false) - baseline;
assert!(
right_twist.y < 0.0,
"right-hand Magnus must point down, got {right_twist:?}"
);
assert!(
left_twist.y > 0.0,
"left-hand Magnus must point up, got {left_twist:?}"
);
assert!((right_twist.y + left_twist.y).abs() < 1e-12);
assert!(right_twist.x.abs() < 1e-12 && right_twist.z.abs() < 1e-12);
assert!(left_twist.x.abs() < 1e-12 && left_twist.z.abs() < 1e-12);
}
#[test]
fn magnus_uses_velocity_corrected_muzzle_stability_gate() {
let muzzle_velocity = 1_400.0 / 3.28084;
let inputs = BallisticInputs {
muzzle_velocity,
bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
bullet_diameter: 0.308 * 0.0254,
bullet_length: 1.215 * 0.0254,
twist_rate: 15.0,
enable_magnus: true,
..BallisticInputs::default()
};
let solver = TrajectorySolver::new(
inputs.clone(),
WindConditions::default(),
AtmosphericConditions::default(),
);
let bare_sg = crate::spin_drift::miller_stability(0.308, 168.0, 15.0, 1.215);
let canonical_sg = solver.effective_spin_drift_sg();
assert!(bare_sg > 1.0, "test requires bare Sg above the Magnus gate");
assert!(
canonical_sg < 1.0,
"velocity-corrected Sg must be below the gate, got {canonical_sg}"
);
let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
let acceleration = solver.calculate_acceleration(
&Vector3::zeros(),
&Vector3::new(muzzle_velocity, 0.0, 0.0),
&Vector3::zeros(),
(temp_c, pressure_hpa, density / 1.225),
);
let mut baseline_inputs = inputs;
baseline_inputs.enable_magnus = false;
let baseline_solver = TrajectorySolver::new(
baseline_inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
let baseline = baseline_solver.calculate_acceleration(
&Vector3::zeros(),
&Vector3::new(muzzle_velocity, 0.0, 0.0),
&Vector3::zeros(),
(temp_c, pressure_hpa, density / 1.225),
);
assert_eq!(
acceleration, baseline,
"canonical Sg below 1 must suppress every Magnus acceleration component"
);
}
#[test]
fn magnus_force_grows_as_fixed_spin_projectile_slows() {
let inputs = BallisticInputs {
muzzle_velocity: 800.0,
bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
bullet_diameter: 0.308 * 0.0254,
bullet_length: 1.215 * 0.0254,
twist_rate: 12.0,
enable_magnus: true,
..BallisticInputs::default()
};
let magnus_acceleration = |speed_mps| {
let evaluate = |enable_magnus| {
let mut run_inputs = inputs.clone();
run_inputs.enable_magnus = enable_magnus;
let solver = TrajectorySolver::new(
run_inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
solver
.calculate_acceleration(
&Vector3::zeros(),
&Vector3::new(speed_mps, 0.0, 0.0),
&Vector3::zeros(),
(temp_c, pressure_hpa, density / 1.225),
)
.y
};
(evaluate(true) - evaluate(false)).abs()
};
let fast = magnus_acceleration(200.0);
let slow = magnus_acceleration(100.0);
let ratio = slow / fast;
let expected_ratio = 2.0_f64.powf(5.0 / 3.0);
assert!(fast > 0.0 && slow > 0.0, "fast={fast}, slow={slow}");
assert!(
(ratio - expected_ratio).abs() < 1e-3,
"fixed-spin Magnus acceleration must grow downrange; slow/fast={ratio}, \
expected={expected_ratio}"
);
}
}
#[cfg(test)]
mod coriolis_direction_tests {
use super::*;
use std::f64::consts::FRAC_PI_2;
#[test]
fn supersonic_crossing_flags_a_positive_range_sample() {
use crate::trajectory_sampling::TrajectoryFlag;
for (solver_name, use_rk4, use_adaptive_rk45) in [
("Euler", false, false),
("RK4", true, false),
("RK45", true, true),
] {
let inputs = BallisticInputs {
muzzle_velocity: 850.0,
bc_value: 0.2,
bc_type: DragModel::G7,
muzzle_angle: 0.03,
enable_trajectory_sampling: true,
sample_interval: 50.0,
use_rk4,
use_adaptive_rk45,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(2000.0);
let samples = solver
.solve()
.expect("supersonic solve should succeed")
.sampled_points
.expect("sampling was enabled");
let flagged_distances: Vec<_> = samples
.iter()
.filter(|sample| sample.flags.contains(&TrajectoryFlag::MachTransition))
.map(|sample| sample.distance_m)
.collect();
assert!(
!flagged_distances.is_empty()
&& flagged_distances.iter().all(|distance| *distance > 0.0),
"{solver_name} must flag genuine crossings only at positive range: {flagged_distances:?}"
);
}
}
#[test]
fn subsonic_launch_does_not_flag_a_muzzle_transition() {
use crate::trajectory_sampling::TrajectoryFlag;
for (solver_name, use_rk4, use_adaptive_rk45) in [
("Euler", false, false),
("RK4", true, false),
("RK45", true, true),
] {
let inputs = BallisticInputs {
muzzle_velocity: 250.0,
muzzle_angle: 0.02,
enable_trajectory_sampling: true,
sample_interval: 25.0,
use_rk4,
use_adaptive_rk45,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(300.0);
let samples = solver
.solve()
.expect("subsonic solve should succeed")
.sampled_points
.expect("sampling was enabled");
assert!(
samples
.iter()
.all(|sample| !sample.flags.contains(&TrajectoryFlag::MachTransition)),
"{solver_name} marked a Mach transition for a launch already below Mach 1"
);
}
}
#[test]
fn mach_transition_tracker_requires_a_downward_crossing() {
fn record(mach_values: &[f64]) -> Vec<f64> {
let mut tracker = MachTransitionTracker::default();
let mut distances = Vec::new();
for (index, mach) in mach_values.iter().copied().enumerate() {
tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
}
distances
}
assert!(record(&[0.9, 0.8, 0.7]).is_empty());
assert_eq!(record(&[1.1, 1.05, 0.99]), vec![20.0]);
assert_eq!(record(&[1.2, 1.19, 1.0, 0.99]), vec![10.0, 30.0]);
assert_eq!(record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]), vec![20.0, 30.0]);
assert!(record(&[1.3, f64::NAN, 1.1]).is_empty());
}
#[test]
fn mach_transition_tracker_labels_0_9_without_touching_the_flat_vec() {
fn record(mach_values: &[f64]) -> (Vec<f64>, MachTransitionTracker) {
let mut tracker = MachTransitionTracker::default();
let mut distances = Vec::new();
for (index, mach) in mach_values.iter().copied().enumerate() {
tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
}
(distances, tracker)
}
let (distances, tracker) = record(&[0.9, 0.8, 0.7]);
assert!(distances.is_empty()); assert_eq!(tracker.mach_1_2_distance_m, None);
assert_eq!(tracker.mach_1_0_distance_m, None);
assert_eq!(tracker.mach_0_9_distance_m, Some(10.0));
let (distances, tracker) = record(&[1.1, 1.05, 0.99]);
assert_eq!(distances, vec![20.0]);
assert_eq!(tracker.mach_1_2_distance_m, None);
assert_eq!(tracker.mach_1_0_distance_m, Some(20.0));
assert_eq!(tracker.mach_0_9_distance_m, None);
let (distances, tracker) = record(&[1.2, 1.19, 1.0, 0.99]);
assert_eq!(distances, vec![10.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(10.0));
assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
assert_eq!(tracker.mach_0_9_distance_m, None);
let (distances, tracker) = record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]);
assert_eq!(distances, vec![20.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(20.0));
assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
assert_eq!(tracker.mach_0_9_distance_m, Some(50.0));
assert!(
tracker.mach_1_2_distance_m < tracker.mach_1_0_distance_m
&& tracker.mach_1_0_distance_m < tracker.mach_0_9_distance_m,
"labeled crossings must be strictly increasing downrange"
);
let (distances, tracker) = record(&[1.3, f64::NAN, 1.1]);
assert!(distances.is_empty());
assert_eq!(tracker.mach_1_2_distance_m, None);
assert_eq!(tracker.mach_1_0_distance_m, None);
assert_eq!(tracker.mach_0_9_distance_m, None);
}
#[test]
fn humidity_percent_converts_and_clamps() {
let mut i = BallisticInputs {
humidity: 0.5,
..BallisticInputs::default()
};
assert!((i.humidity_percent() - 50.0).abs() < 1e-9, "0.5 -> 50%");
i.humidity = 0.0;
assert_eq!(i.humidity_percent(), 0.0);
i.humidity = 1.0;
assert_eq!(i.humidity_percent(), 100.0);
i.humidity = 1.5; assert_eq!(i.humidity_percent(), 100.0);
}
fn vertical_at(shot_azimuth: f64, range_m: f64) -> f64 {
let inputs = BallisticInputs {
muzzle_velocity: 800.0,
bc_value: 0.5,
bc_type: DragModel::G7,
muzzle_angle: 0.02, enable_coriolis: true,
latitude: Some(45.0),
shot_azimuth,
ground_threshold: f64::NEG_INFINITY, ..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(range_m + 50.0);
let r = solver.solve().expect("solve");
let pts = &r.points;
for i in 1..pts.len() {
if pts[i].position.x >= range_m {
let p1 = &pts[i - 1];
let p2 = &pts[i];
let t = (range_m - p1.position.x) / (p2.position.x - p1.position.x);
return p1.position.y + t * (p2.position.y - p1.position.y);
}
}
panic!("range {range_m} not reached");
}
#[test]
fn eotvos_east_higher_than_west() {
let range = 600.0;
let east = vertical_at(FRAC_PI_2, range); let west = vertical_at(3.0 * FRAC_PI_2, range); let north = vertical_at(0.0, range); assert!(
east > west,
"east ({east:.5}) must be higher than west ({west:.5}) at {range} m (Eötvös)"
);
assert!(
east > north && north > west,
"north ({north:.5}) must lie between east ({east:.5}) and west ({west:.5})"
);
assert!(
(east - west) > 1e-3,
"E-W vertical separation ({:.6} m) should be physically meaningful, not FP noise",
east - west
);
}
#[test]
fn labeled_mach_crossings_match_pinned_pre_change_flat_vec_across_solvers() {
let cases = [
("Euler", false, false, 670.9878683238721_f64, 805.5274119916264_f64),
("RK4", true, false, 671.7257336844475_f64, 805.933409072171_f64),
("RK45", true, true, 672.4905711917901_f64, 806.5709746782849_f64),
];
for (solver_name, use_rk4, use_adaptive_rk45, expected_1_2, expected_1_0) in cases {
let inputs = BallisticInputs {
muzzle_velocity: 850.0,
bc_value: 0.2,
bc_type: DragModel::G7,
muzzle_angle: 0.03,
use_rk4,
use_adaptive_rk45,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(2000.0);
let result = solver.solve().expect("solve should succeed");
assert_eq!(
result.mach_1_2_distance_m,
Some(expected_1_2),
"{solver_name}: mach_1_2_distance_m must match the pinned pre-change flat-Vec value"
);
assert_eq!(
result.mach_1_0_distance_m,
Some(expected_1_0),
"{solver_name}: mach_1_0_distance_m must match the pinned pre-change flat-Vec value"
);
let mach_1_2 = result.mach_1_2_distance_m.expect("crosses 1.2");
let mach_1_0 = result.mach_1_0_distance_m.expect("crosses 1.0");
let mach_0_9 = result
.mach_0_9_distance_m
.expect("this trajectory also goes past 0.9 within 2000 m");
assert!(
mach_1_2 < mach_1_0 && mach_1_0 < mach_0_9,
"{solver_name}: labeled crossings must be strictly increasing downrange \
(1.2={mach_1_2}, 1.0={mach_1_0}, 0.9={mach_0_9})"
);
}
}
#[test]
fn labeled_mach_crossings_are_none_for_a_fully_supersonic_trajectory() {
for (solver_name, use_rk4, use_adaptive_rk45) in [
("Euler", false, false),
("RK4", true, false),
("RK45", true, true),
] {
let inputs = BallisticInputs {
muzzle_velocity: 850.0,
bc_value: 0.2,
bc_type: DragModel::G7,
muzzle_angle: 0.03,
use_rk4,
use_adaptive_rk45,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
solver.set_max_range(200.0);
let result = solver.solve().expect("solve should succeed");
assert_eq!(
result.mach_1_2_distance_m, None,
"{solver_name}: must not report a 1.2 crossing that never happens"
);
assert_eq!(
result.mach_1_0_distance_m, None,
"{solver_name}: must not report a 1.0 crossing that never happens"
);
assert_eq!(
result.mach_0_9_distance_m, None,
"{solver_name}: must not report a 0.9 crossing that never happens"
);
}
}
}
#[cfg(test)]
mod cant_tests {
use super::*;
fn base_inputs() -> BallisticInputs {
BallisticInputs {
muzzle_velocity: 800.0,
bc_value: 0.5,
bc_type: DragModel::G7,
bullet_mass: 0.0109,
bullet_diameter: 0.00782,
bullet_length: 0.0309,
sight_height: 0.05,
twist_rate: 10.0,
use_rk4: true,
..BallisticInputs::default()
}
}
fn solve_with(inputs: BallisticInputs, max_range: f64) -> TrajectoryResult {
let mut s = TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
);
s.set_max_range(max_range);
s.solve().expect("solve")
}
fn yz_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
let pts = &result.points;
for i in 1..pts.len() {
if pts[i].position.x >= x {
let (p1, p2) = (&pts[i - 1], &pts[i]);
let dx = p2.position.x - p1.position.x;
let t = if dx.abs() < 1e-12 { 0.0 } else { (x - p1.position.x) / dx };
return (
p1.position.y + t * (p2.position.y - p1.position.y),
p1.position.z + t * (p2.position.z - p1.position.z),
);
}
}
panic!("trajectory never reached {x} m");
}
#[test]
fn cant_sign_clockwise_up_offset_goes_right_and_low() {
let mut level = base_inputs();
level.muzzle_angle = 0.003; let mut canted = level.clone();
canted.cant_angle = 10f64.to_radians();
let (y0, z0) = yz_at(&solve_with(level, 400.0), 300.0);
let (y1, z1) = yz_at(&solve_with(canted, 400.0), 300.0);
assert!(z1 > z0 + 0.01, "clockwise cant must move POI right: z0={z0} z1={z1}");
assert!(y1 < y0 - 0.001, "clockwise cant must move POI low: y0={y0} y1={y1}");
}
#[test]
fn pure_cant_shows_bore_offset_near_range() {
let mut i = base_inputs();
i.muzzle_angle = 0.0;
i.cant_angle = 10f64.to_radians();
let sh = i.sight_height;
let r = solve_with(i, 60.0);
let first = &r.points[1]; let expected = -sh * 10f64.to_radians().sin();
assert!(
(first.position.z - expected).abs() < 0.005,
"near-muzzle lateral {} should be ~bore offset {expected}",
first.position.z
);
}
#[test]
fn zero_angle_is_independent_of_cant() {
let a = base_inputs();
let mut b = base_inputs();
b.cant_angle = 15f64.to_radians();
let za = calculate_zero_angle(a.clone(), 100.0, 0.0).expect("zero a");
let zb = calculate_zero_angle(b.clone(), 100.0, 0.0).expect("zero b");
assert_eq!(za.to_bits(), zb.to_bits(), "zeroing must ignore cant: {za} vs {zb}");
let _ = (a.cant_angle, b.cant_angle);
}
#[test]
fn nonfinite_cant_is_rejected() {
let mut i = base_inputs();
i.cant_angle = f64::NAN;
let s = TrajectorySolver::new(i, WindConditions::default(), AtmosphericConditions::default());
assert!(s.solve().is_err());
}
#[test]
fn incline_and_cant_compose_without_breaking() {
let mut flat = base_inputs();
flat.muzzle_angle = 0.003;
flat.shooting_angle = 15f64.to_radians();
let mut canted = flat.clone();
canted.cant_angle = 10f64.to_radians();
let (_, z_flat) = yz_at(&solve_with(flat, 400.0), 300.0);
let (_, z_cant) = yz_at(&solve_with(canted, 400.0), 300.0);
assert!(z_cant > z_flat, "cant must still deflect right on an incline");
}
}
#[cfg(test)]
mod vertical_wind_tests {
use super::*;
fn base_inputs() -> BallisticInputs {
BallisticInputs {
muzzle_velocity: 800.0,
bc_value: 0.5,
bc_type: DragModel::G7,
bullet_mass: 0.0109,
bullet_diameter: 0.00782,
bullet_length: 0.0309,
sight_height: 0.05,
twist_rate: 10.0,
use_rk4: true,
..BallisticInputs::default()
}
}
fn y_at(result: &TrajectoryResult, x: f64) -> f64 {
let pts = &result.points;
for i in 1..pts.len() {
if pts[i].position.x >= x {
let (p1, p2) = (&pts[i - 1], &pts[i]);
let dx = p2.position.x - p1.position.x;
let t = if dx.abs() < 1e-12 { 0.0 } else { (x - p1.position.x) / dx };
return p1.position.y + t * (p2.position.y - p1.position.y);
}
}
panic!("trajectory never reached {x} m");
}
fn solve_with(inputs: BallisticInputs, wind: WindConditions, max_range: f64) -> TrajectoryResult {
let mut s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
s.set_max_range(max_range);
s.solve().expect("solve")
}
#[test]
fn updraft_raises_poi_downrange() {
let calm_inputs = base_inputs();
let calm_wind = WindConditions::default();
let updraft = WindConditions {
vertical_speed: 5.0,
..Default::default()
};
let calm = solve_with(calm_inputs.clone(), calm_wind, 500.0);
let updraft_result = solve_with(calm_inputs, updraft, 500.0);
let y_calm = y_at(&calm, 400.0);
let y_updraft = y_at(&updraft_result, 400.0);
assert!(
y_updraft > y_calm,
"5 m/s updraft must raise POI at 400m: calm={y_calm}, updraft={y_updraft}"
);
}
#[test]
fn zero_vertical_is_default_and_finite_required() {
assert_eq!(WindConditions::default().vertical_speed, 0.0);
let inputs = base_inputs();
let wind = WindConditions {
vertical_speed: f64::NAN,
..Default::default()
};
let s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
assert!(
s.solve().is_err(),
"NaN wind.vertical_speed must be rejected by validate_for_solve"
);
}
}
#[cfg(test)]
mod bc_reference_standard_tests {
use super::*;
fn base_inputs() -> BallisticInputs {
BallisticInputs {
muzzle_velocity: 800.0,
bc_value: 0.5,
bc_type: DragModel::G7,
bullet_mass: 0.0109,
bullet_diameter: 0.00782,
bullet_length: 0.0309,
sight_height: 0.05,
twist_rate: 10.0,
use_rk4: true,
..BallisticInputs::default()
}
}
fn y_and_speed_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
let pts = &result.points;
for i in 1..pts.len() {
if pts[i].position.x >= x {
let (p1, p2) = (&pts[i - 1], &pts[i]);
let dx = p2.position.x - p1.position.x;
let t = if dx.abs() < 1e-12 {
0.0
} else {
(x - p1.position.x) / dx
};
return (
p1.position.y + t * (p2.position.y - p1.position.y),
p1.velocity_magnitude + t * (p2.velocity_magnitude - p1.velocity_magnitude),
);
}
}
panic!("trajectory never reached {x} m");
}
#[test]
fn asm_to_icao_ratio_matches_documented_value() {
assert!(
(crate::constants::ASM_TO_ICAO_BC - 0.98237).abs() < 1e-5,
"ASM_TO_ICAO_BC = {} must equal 0.98237 to 5 decimal places",
crate::constants::ASM_TO_ICAO_BC
);
assert_eq!(
crate::constants::ASM_TO_ICAO_BC,
crate::constants::ASM_DENSITY_LB_FT3 / crate::constants::ICAO_DENSITY_LB_FT3
);
}
#[test]
fn default_bc_reference_standard_is_icao() {
assert_eq!(
BallisticInputs::default().bc_reference_standard,
BcReferenceStandard::Icao
);
}
#[test]
fn icao_reference_leaves_bc_value_bit_identical() {
let raw_bc: f64 = 0.4372911; let inputs = BallisticInputs {
bc_value: raw_bc,
bc_reference_standard: BcReferenceStandard::Icao,
..base_inputs()
};
let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
assert_eq!(solver.inputs.bc_value.to_bits(), raw_bc.to_bits());
}
#[test]
fn default_inputs_solve_is_unaffected_by_the_new_field_existing() {
let a = TrajectorySolver::new(base_inputs(), WindConditions::default(), AtmosphericConditions::default())
.solve()
.expect("solve a");
let b = TrajectorySolver::new(
BallisticInputs { ..base_inputs() },
WindConditions::default(),
AtmosphericConditions::default(),
)
.solve()
.expect("solve b");
assert_eq!(a.impact_velocity.to_bits(), b.impact_velocity.to_bits());
assert_eq!(a.max_range.to_bits(), b.max_range.to_bits());
}
#[test]
fn army_standard_metro_scales_bc_value_by_exactly_the_derived_ratio() {
let raw_bc = 0.5;
let inputs = BallisticInputs {
bc_value: raw_bc,
bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
..base_inputs()
};
let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
assert_eq!(
solver.inputs.bc_value,
raw_bc * crate::constants::ASM_TO_ICAO_BC
);
}
#[test]
fn army_standard_metro_scales_mach_keyed_bc_segments() {
let inputs = BallisticInputs {
bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
bc_segments: Some(vec![(0.5, 0.40), (1.5, 0.30)]),
..base_inputs()
};
let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
let segments = solver.inputs.bc_segments.as_ref().expect("segments");
assert_eq!(segments[0], (0.5, 0.40 * crate::constants::ASM_TO_ICAO_BC));
assert_eq!(segments[1], (1.5, 0.30 * crate::constants::ASM_TO_ICAO_BC));
}
#[test]
fn army_standard_metro_scales_velocity_keyed_bc_segments_data() {
let inputs = BallisticInputs {
bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
bc_segments_data: Some(vec![
crate::BCSegmentData {
velocity_min: 0.0,
velocity_max: 500.0,
bc_value: 0.40,
},
crate::BCSegmentData {
velocity_min: 500.0,
velocity_max: 900.0,
bc_value: 0.45,
},
]),
..base_inputs()
};
let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
let segments = solver.inputs.bc_segments_data.as_ref().expect("segments");
assert_eq!(segments[0].bc_value, 0.40 * crate::constants::ASM_TO_ICAO_BC);
assert_eq!(segments[1].bc_value, 0.45 * crate::constants::ASM_TO_ICAO_BC);
assert_eq!(segments[0].velocity_min, 0.0);
assert_eq!(segments[1].velocity_max, 900.0);
}
#[test]
fn army_standard_metro_moves_impact_in_the_more_drag_direction() {
let solve_at = |standard: BcReferenceStandard| {
let inputs = BallisticInputs {
bc_value: 0.475,
bc_reference_standard: standard,
..base_inputs()
};
let mut solver =
TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
solver.set_max_range(500.0);
solver.solve().expect("solve")
};
let icao = solve_at(BcReferenceStandard::Icao);
let asm = solve_at(BcReferenceStandard::ArmyStandardMetro);
let (y_icao, v_icao) = y_and_speed_at(&icao, 400.0);
let (y_asm, v_asm) = y_and_speed_at(&asm, 400.0);
assert!(
y_asm < y_icao,
"ArmyStandardMetro must drop MORE (lower y) at 400m than Icao for the same raw \
bc_value: icao_y={y_icao}, asm_y={y_asm}"
);
assert!(
v_asm < v_icao,
"ArmyStandardMetro must retain LESS velocity at 400m than Icao for the same raw \
bc_value: icao_v={v_icao}, asm_v={v_asm}"
);
}
#[test]
fn monte_carlo_inherits_the_normalized_bc_reference() {
let base_inputs_asm = BallisticInputs {
bc_value: 0.475,
bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
..base_inputs()
};
let wind = WindConditions::default();
let mut direct_solver =
TrajectorySolver::new(base_inputs_asm.clone(), wind.clone(), AtmosphericConditions::default());
direct_solver.set_max_range(base_inputs_asm.target_distance.max(1000.0) * 2.0);
let direct = direct_solver.solve().expect("direct solve");
let mc_params = MonteCarloParams {
num_simulations: 1,
velocity_std_dev: 0.0,
angle_std_dev: 0.0,
bc_std_dev: 0.0,
wind_speed_std_dev: 0.0,
target_distance: None,
base_wind_speed: 0.0,
base_wind_direction: 0.0,
azimuth_std_dev: 0.0,
};
let mc = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
base_inputs_asm,
wind,
mc_params,
0.0,
42,
)
.expect("monte carlo");
assert_eq!(mc.ranges.len(), 1);
assert_eq!(
mc.ranges[0].to_bits(),
direct.max_range.to_bits(),
"a zero-dispersion single MC sample must match a plain solve of the same \
ASM-referenced inputs bit-for-bit"
);
assert_eq!(
mc.impact_velocities[0].to_bits(),
direct.impact_velocity.to_bits()
);
}
#[test]
fn estimate_bc_fit_recovers_an_icao_referenced_bc() {
let known_bc = 0.475;
let velocity = 800.0;
let mass = 0.0109;
let diameter = 0.00782;
let atmosphere = AtmosphericConditions::default();
let synth_inputs = BallisticInputs {
muzzle_velocity: velocity,
bc_value: known_bc,
bc_type: DragModel::G7,
bullet_mass: mass,
bullet_diameter: diameter,
bullet_length: 0.0309,
sight_height: 0.05,
twist_rate: 10.0,
use_rk4: true,
bc_reference_standard: BcReferenceStandard::Icao,
..BallisticInputs::default()
};
let mut solver = TrajectorySolver::new(synth_inputs, WindConditions::default(), atmosphere.clone());
solver.set_max_range(500.0);
let trajectory = solver.solve().expect("synthetic solve");
let points: Vec<(f64, f64)> = [100.0, 200.0, 300.0, 400.0]
.iter()
.map(|&d| {
let (y, _) = {
let pts = &trajectory.points;
let mut found = None;
for i in 1..pts.len() {
if pts[i].position.x >= d {
let (p1, p2) = (&pts[i - 1], &pts[i]);
let dx = p2.position.x - p1.position.x;
let t = if dx.abs() < 1e-12 {
0.0
} else {
(d - p1.position.x) / dx
};
found = Some((
p1.position.y + t * (p2.position.y - p1.position.y),
0.0,
));
break;
}
}
found.expect("trajectory reached observation distance")
};
(d, -y) })
.collect();
let estimate = estimate_bc_fit(
velocity,
mass,
diameter,
&points,
DragModel::G7,
BcFitMode::Drop,
atmosphere,
None,
0.05,
)
.expect("fit should converge");
assert!(
(estimate.bc - known_bc).abs() < 0.02,
"fit should recover the known ICAO-referenced bc={known_bc}, got {}",
estimate.bc
);
}
#[test]
fn custom_drag_table_makes_bc_reference_standard_numerically_inert() {
let table = crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0, 3.0], vec![0.3, 0.4, 0.3, 0.2])
.expect("valid table");
let solve_with = |standard: BcReferenceStandard| {
let inputs = BallisticInputs {
bc_value: 0.5, bc_reference_standard: standard,
custom_drag_table: Some(table.clone()),
..base_inputs()
};
let mut solver =
TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
solver.set_max_range(500.0);
solver.solve().expect("solve")
};
let icao = solve_with(BcReferenceStandard::Icao);
let asm = solve_with(BcReferenceStandard::ArmyStandardMetro);
assert_eq!(
icao.impact_velocity.to_bits(),
asm.impact_velocity.to_bits(),
"a custom drag table must make bc_reference_standard fully inert"
);
assert_eq!(icao.max_range.to_bits(), asm.max_range.to_bits());
}
#[test]
fn custom_drag_table_inert_warning_fires_only_for_army_standard_metro_with_a_table() {
let table = crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0], vec![0.3, 0.4, 0.3])
.expect("valid table");
let no_table_icao = base_inputs();
assert!(no_table_icao.bc_reference_standard_inert_warning().is_none());
let no_table_asm = BallisticInputs {
bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
..base_inputs()
};
assert!(no_table_asm.bc_reference_standard_inert_warning().is_none());
let table_icao = BallisticInputs {
custom_drag_table: Some(table.clone()),
..base_inputs()
};
assert!(table_icao.bc_reference_standard_inert_warning().is_none());
let table_asm = BallisticInputs {
custom_drag_table: Some(table),
bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
..base_inputs()
};
let warning = table_asm
.bc_reference_standard_inert_warning()
.expect("must warn");
assert!(warning.contains("--bc-reference"));
assert!(warning.contains("--drag-table"));
}
}
#[cfg(test)]
mod effective_drag_coefficient_tests {
use super::*;
fn inputs_175gr_g7() -> BallisticInputs {
let mut inputs = BallisticInputs {
bc_value: 0.243,
bc_type: DragModel::G7,
muzzle_velocity: 823.0,
..Default::default()
};
inputs.bullet_mass = 175.0 * crate::constants::GRAINS_TO_KG;
inputs.bullet_diameter = 0.308 * 0.0254;
inputs.weight_grains = 175.0;
inputs.caliber_inches = 0.308;
inputs
}
fn solver(inputs: BallisticInputs) -> TrajectorySolver {
TrajectorySolver::new(
inputs,
WindConditions::default(),
AtmosphericConditions::default(),
)
}
#[test]
fn reports_the_projectiles_own_cd_not_the_reference_tables() {
let inputs = inputs_175gr_g7();
let sd = inputs.sectional_density_lb_in2().expect("SD");
let solver = solver(inputs);
let sos = 340.0;
let velocity = 800.0;
let mach = velocity / sos;
let reference = crate::drag::get_drag_coefficient(mach, &DragModel::G7);
let reported = solver
.effective_drag_coefficient(velocity, sos)
.expect("mass and diameter are set");
let expected = reference * sd / 0.243;
assert!(
(reported - expected).abs() < 1e-12,
"reported {reported} != Cd_ref * SD / BC {expected}"
);
assert!(
(reported - reference).abs() > 1e-6,
"form factor collapsed to 1; this fixture no longer distinguishes the two values"
);
}
#[test]
fn a_custom_drag_table_passes_through_unscaled() {
let mut inputs = inputs_175gr_g7();
inputs.custom_drag_table = Some(crate::drag::DragTable::new(
vec![0.5, 3.0],
vec![0.15, 0.40],
));
let solver = solver(inputs);
let sos = 340.0;
let velocity = 0.9 * sos;
let table_value = solver
.inputs
.custom_drag_table
.as_ref()
.expect("table")
.interpolate(0.9);
let reported = solver
.effective_drag_coefficient(velocity, sos)
.expect("mass and diameter are set");
assert!(
(reported - table_value).abs() < 1e-12,
"custom table Cd {table_value} was rescaled to {reported}"
);
}
#[test]
fn a_velocity_segmented_bc_steps_the_reported_cd() {
let mut inputs = inputs_175gr_g7();
inputs.use_bc_segments = true;
inputs.bc_segments_data = Some(vec![
crate::BCSegmentData { velocity_min: 2400.0, velocity_max: 4000.0, bc_value: 0.243 },
crate::BCSegmentData { velocity_min: 0.0, velocity_max: 2400.0, bc_value: 0.200 },
]);
let solver = solver(inputs);
let sos = 340.0;
let above = solver.effective_drag_coefficient(2500.0 / 3.28084, sos).expect("cd");
let below = solver.effective_drag_coefficient(2300.0 / 3.28084, sos).expect("cd");
assert!(
below > above,
"expected the 0.200 band to report a higher Cd than the 0.243 band; got {below} vs {above}"
);
}
#[test]
fn is_absent_when_sectional_density_is_unknown() {
let mut inputs = inputs_175gr_g7();
inputs.weight_grains = 0.0;
inputs.bullet_mass = 0.0;
let solver = solver(inputs);
assert!(solver.effective_drag_coefficient(800.0, 340.0).is_none());
}
#[test]
fn the_json_emit_rule_is_flag_gated_and_absent_when_cd_is_unknown() {
let mut point = TrajectoryPoint {
time: 0.0,
position: nalgebra::Vector3::new(0.0, 0.0, 0.0),
velocity_magnitude: 800.0,
kinetic_energy: 3000.0,
drag_coefficient: Some(0.31),
};
assert_eq!(point.drag_coefficient_json_value(true), Some(0.31));
assert_eq!(
point.drag_coefficient_json_value(false),
None,
"without the flag the key must not exist, so default JSON stays byte-identical"
);
point.drag_coefficient = None;
assert_eq!(
point.drag_coefficient_json_value(true),
None,
"unknown sectional density must yield an ABSENT key, not null"
);
}
#[test]
fn every_point_of_a_solved_trajectory_carries_the_value() {
let mut solver = solver(inputs_175gr_g7());
solver.set_max_range(300.0);
let result = solver.solve().expect("solve");
assert!(!result.points.is_empty());
assert!(
result.points.iter().all(|p| p.drag_coefficient.is_some()),
"the post-integration pass missed at least one point"
);
}
}