1use crate::cluster_bc::ClusterBCDegradation;
3use crate::mc_stats::{wilson_interval, BernoulliConfidenceSequence, ConfidenceLevel, Welford};
4use crate::pitch_damping::{calculate_pitch_damping_coefficient, PitchDampingCoefficients};
5use crate::precession_nutation::{
6 calculate_combined_angular_motion, projectile_moments_of_inertia, AngularState,
7 PrecessionNutationParams,
8};
9use crate::trajectory_observation::{bracket_param, Bracket, TrajectoryTermination};
10use crate::trajectory_sampling::{
11 projected_sample_count, sample_trajectory, TrajectoryData, TrajectoryOutputs, TrajectorySample,
12};
13use crate::wind_shear::WindShearModel;
14use crate::DragModel;
15use nalgebra::{Vector3, Vector6};
16use std::error::Error;
17use std::fmt;
18
19#[derive(Debug, Clone, Copy, PartialEq)]
29#[cfg_attr(feature = "cli", derive(clap::ValueEnum))]
30pub enum UnitSystem {
31 Metric,
33 Imperial,
35}
36
37#[derive(Debug, Clone, Copy, PartialEq)]
39pub enum OutputFormat {
40 Table,
41 Json,
42 Csv,
43}
44
45#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
64#[cfg_attr(feature = "cli", derive(clap::ValueEnum))]
65pub enum BcReferenceStandard {
66 #[default]
69 Icao,
70 ArmyStandardMetro,
73}
74
75pub const BC_REFERENCE_STANDARD_INERT_WARNING: &str =
79 "warning: --bc-reference army-standard-metro has no effect together with a custom drag \
80 table (--drag-table): the deck's Cd is divided by sectional density, not a BC value, so \
81 no BC-reference conversion applies";
82
83#[derive(Debug)]
85pub struct BallisticsError {
86 message: String,
87}
88
89impl fmt::Display for BallisticsError {
90 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
91 write!(f, "{}", self.message)
92 }
93}
94
95impl Error for BallisticsError {}
96
97impl From<String> for BallisticsError {
98 fn from(msg: String) -> Self {
99 BallisticsError { message: msg }
100 }
101}
102
103impl From<&str> for BallisticsError {
104 fn from(msg: &str) -> Self {
105 BallisticsError {
106 message: msg.to_string(),
107 }
108 }
109}
110
111#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
119pub enum DropsReference {
120 #[default]
122 Los,
123 Target,
126}
127
128#[derive(Debug, Clone)]
132pub struct BallisticInputs {
133 pub bc_value: f64, pub bc_type: DragModel, pub bc_reference_standard: BcReferenceStandard,
141 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,
155 pub shooting_angle: f64, pub cant_angle: f64,
165 pub sight_height: f64, pub sight_offset_lateral_m: f64,
177 pub muzzle_height: f64, pub target_height: f64, pub zero_poi_vertical_m: f64,
188 pub zero_poi_horizontal_m: f64,
194 pub ground_threshold: f64, pub altitude: f64, pub temperature: f64, pub pressure: f64, pub humidity: f64,
205 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,
227 pub enable_magnus: bool, pub enable_coriolis: bool, pub use_powder_sensitivity: bool,
230 pub powder_temp_sensitivity: f64, pub powder_temp: f64, pub powder_temp_curve: Option<Vec<(f64, f64)>>,
239 pub powder_curve_temp_c: Option<f64>,
243 pub tipoff_yaw: f64, pub tipoff_decay_distance: f64, pub cd_delta2: f64,
251 pub use_bc_segments: bool,
261 pub bc_segments: Option<Vec<(f64, f64)>>, pub bc_segments_data: Option<Vec<crate::BCSegmentData>>, pub use_enhanced_spin_drift: bool,
264 pub use_form_factor: bool,
267 pub enable_wind_shear: bool,
268 pub wind_shear_model: String,
269 pub enable_trajectory_sampling: bool,
270 pub sample_interval: f64, pub drops_reference: DropsReference,
280 pub enable_pitch_damping: bool,
281 pub enable_precession_nutation: bool,
282 pub enable_aerodynamic_jump: bool,
285 pub use_cluster_bc: bool, pub custom_drag_table: Option<crate::drag::DragTable>,
289 pub cd_scale: f64,
297
298 pub bc_type_str: Option<String>,
300}
301
302impl BallisticInputs {
303 pub fn humidity_percent(&self) -> f64 {
308 (self.humidity * 100.0).clamp(0.0, 100.0)
309 }
310
311 pub fn windage_zero_bias_rad(&self, zero_distance_m: f64) -> f64 {
326 if zero_distance_m > 0.0 {
327 (self.zero_poi_horizontal_m + self.sight_offset_lateral_m) / zero_distance_m
328 } else {
329 0.0
330 }
331 }
332
333 pub fn sectional_density_lb_in2(&self) -> Option<f64> {
339 let weight_gr = if self.weight_grains > 0.0 {
340 self.weight_grains
341 } else {
342 self.bullet_mass / crate::constants::GRAINS_TO_KG };
344 let diameter_in = if self.caliber_inches > 0.0 {
345 self.caliber_inches
346 } else {
347 self.bullet_diameter / 0.0254 };
349 if weight_gr > 0.0 && diameter_in > 0.0 {
350 Some(weight_gr / 7000.0 / (diameter_in * diameter_in))
351 } else {
352 None
353 }
354 }
355
356 pub fn custom_drag_denominator(&self, fallback_bc: f64) -> f64 {
368 match self.sectional_density_lb_in2() {
369 Some(sd) => sd,
370 None => {
371 static WARN_ONCE: std::sync::Once = std::sync::Once::new();
372 WARN_ONCE.call_once(|| {
373 eprintln!(
374 "Warning: custom drag table active but bullet mass/diameter are \
375 unavailable; falling back to bc_value for the retardation denominator"
376 );
377 });
378 fallback_bc
379 }
380 }
381 }
382
383 pub fn bc_reference_standard_inert_warning(&self) -> Option<&'static str> {
394 if self.custom_drag_table.is_some()
395 && matches!(
396 self.bc_reference_standard,
397 BcReferenceStandard::ArmyStandardMetro
398 )
399 {
400 Some(BC_REFERENCE_STANDARD_INERT_WARNING)
401 } else {
402 None
403 }
404 }
405
406 pub fn normalize_for_solve(&mut self) {
425 if matches!(
436 self.bc_reference_standard,
437 BcReferenceStandard::ArmyStandardMetro
438 ) {
439 self.bc_value *= crate::constants::ASM_TO_ICAO_BC;
440 if let Some(segments) = self.bc_segments.as_mut() {
441 for (_mach, bc) in segments.iter_mut() {
442 *bc *= crate::constants::ASM_TO_ICAO_BC;
443 }
444 }
445 if let Some(segments) = self.bc_segments_data.as_mut() {
446 for segment in segments.iter_mut() {
447 segment.bc_value *= crate::constants::ASM_TO_ICAO_BC;
448 }
449 }
450 self.bc_reference_standard = BcReferenceStandard::Icao;
453 }
454
455 self.caliber_inches = self.bullet_diameter / 0.0254;
460 self.weight_grains = self.bullet_mass / crate::constants::GRAINS_TO_KG;
461
462 self.muzzle_velocity = resolve_powder_adjusted_velocity(
473 self.muzzle_velocity,
474 self.temperature,
475 self.use_powder_sensitivity,
476 self.powder_temp_sensitivity,
477 self.powder_temp,
478 self.powder_temp_curve.as_deref(),
479 self.powder_curve_temp_c,
480 );
481 }
482}
483
484impl Default for BallisticInputs {
485 fn default() -> Self {
486 let mass_kg = 0.01;
487 let diameter_m = 0.00762;
488 let bc = 0.5;
489 let muzzle_angle_rad = 0.0;
490 let bc_type = DragModel::G1;
491
492 Self {
493 bc_value: bc,
495 bc_type,
496 bc_reference_standard: BcReferenceStandard::Icao,
497 bullet_mass: mass_kg,
498 muzzle_velocity: 800.0,
499 bullet_diameter: diameter_m,
500 bullet_length: crate::stability::estimate_bullet_length_m(diameter_m, mass_kg),
504
505 muzzle_angle: muzzle_angle_rad,
507 target_distance: 100.0,
508 azimuth_angle: 0.0,
509 shot_azimuth: 0.0,
510 shooting_angle: 0.0,
511 cant_angle: 0.0,
512 sight_height: 0.05,
513 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,
518 ground_threshold: -100.0, altitude: 0.0,
522 temperature: 15.0,
523 pressure: 1013.25, humidity: 0.5, latitude: None,
526
527 wind_speed: 0.0,
529 wind_angle: 0.0,
530
531 twist_rate: 12.0, is_twist_right: true,
534 caliber_inches: diameter_m / 0.0254, weight_grains: mass_kg / crate::constants::GRAINS_TO_KG, manufacturer: None,
537 bullet_model: None,
538 bullet_id: None,
539 bullet_cluster: None,
540
541 use_rk4: true, use_adaptive_rk45: true, enable_advanced_effects: false,
547 enable_magnus: false,
548 enable_coriolis: false,
549 use_powder_sensitivity: false,
550 powder_temp_sensitivity: 0.0,
551 powder_temp: 15.0,
552 powder_temp_curve: None,
553 powder_curve_temp_c: None,
554 tipoff_yaw: 0.0,
555 tipoff_decay_distance: 50.0,
556 cd_delta2: 7.5,
557 use_bc_segments: false,
558 bc_segments: None,
559 bc_segments_data: None,
560 use_enhanced_spin_drift: false,
561 use_form_factor: false,
562 enable_wind_shear: false,
563 wind_shear_model: "none".to_string(),
564 enable_trajectory_sampling: false,
565 sample_interval: 10.0, drops_reference: DropsReference::Los, enable_pitch_damping: false,
568 enable_precession_nutation: false,
569 enable_aerodynamic_jump: false,
570 use_cluster_bc: false, custom_drag_table: None,
574 cd_scale: 1.0,
575
576 bc_type_str: None,
578 }
579 }
580}
581
582pub fn interpolate_powder_temp_curve(curve: &[(f64, f64)], temp_c: f64) -> f64 {
588 debug_assert!(!curve.is_empty());
589 if curve.is_empty() {
590 return 0.0;
591 }
592 let mut sorted;
595 let pts: &[(f64, f64)] = if curve.windows(2).all(|w| w[0].0 <= w[1].0) {
596 curve
597 } else {
598 sorted = curve.to_vec();
599 sorted.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
600 &sorted
601 };
602 let n = pts.len();
603 if temp_c <= pts[0].0 {
604 return pts[0].1; }
606 if temp_c >= pts[n - 1].0 {
607 return pts[n - 1].1; }
609 for i in 1..n {
610 let (t0, v0) = pts[i - 1];
611 let (t1, v1) = pts[i];
612 if temp_c <= t1 {
613 let span = t1 - t0;
614 if span.abs() < f64::EPSILON {
615 return v1; }
617 let f = (temp_c - t0) / span;
618 return v0 + f * (v1 - v0);
619 }
620 }
621 pts[n - 1].1
622}
623
624pub fn parse_powder_sweep(s: &str) -> Result<Vec<f64>, String> {
630 const MAX_SWEEP_ROWS: usize = 500;
631 let parts: Vec<&str> = s.split(':').collect();
632 if parts.len() != 3 {
633 return Err(format!(
634 "Invalid --sweep '{}': expected START:END:STEP (e.g. \"20:110:10\")",
635 s
636 ));
637 }
638 let parse = |p: &str, name: &str| -> Result<f64, String> {
639 p.trim()
640 .parse::<f64>()
641 .map_err(|_| format!("Invalid --sweep {}: '{}' is not a number", name, p.trim()))
642 };
643 let start = parse(parts[0], "START")?;
644 let end = parse(parts[1], "END")?;
645 let step = parse(parts[2], "STEP")?;
646 if !step.is_finite() || step <= 0.0 {
647 return Err(format!("Invalid --sweep STEP {}: must be positive", step));
648 }
649 if !start.is_finite() || !end.is_finite() || end < start {
650 return Err(format!(
651 "Invalid --sweep range {}:{}: END must be >= START",
652 start, end
653 ));
654 }
655 let n_f = ((end - start) / step + 1e-9).floor();
661 if !n_f.is_finite() || n_f + 1.0 > MAX_SWEEP_ROWS as f64 {
662 return Err(format!(
663 "--sweep would produce more than {} rows; use a larger STEP",
664 MAX_SWEEP_ROWS
665 ));
666 }
667 let n = n_f as usize + 1;
668 Ok((0..n).map(|i| start + step * i as f64).collect())
670}
671
672pub fn resolve_powder_adjusted_velocity(
681 nominal_velocity_mps: f64,
682 ambient_temperature_c: f64,
683 use_powder_sensitivity: bool,
684 powder_temp_sensitivity_mps_per_c: f64,
685 powder_reference_temp_c: f64,
686 powder_temp_curve: Option<&[(f64, f64)]>,
687 powder_curve_temp_c: Option<f64>,
688) -> f64 {
689 if let Some(curve) = powder_temp_curve {
690 if !curve.is_empty() {
691 let lookup_c = powder_curve_temp_c.unwrap_or(ambient_temperature_c);
692 return interpolate_powder_temp_curve(curve, lookup_c);
693 }
694 return nominal_velocity_mps;
697 }
698 if use_powder_sensitivity {
699 let temp_delta_c = ambient_temperature_c - powder_reference_temp_c;
700 return nominal_velocity_mps + powder_temp_sensitivity_mps_per_c * temp_delta_c;
701 }
702 nominal_velocity_mps
703}
704
705#[derive(Debug, Clone)]
707pub struct WindConditions {
708 pub speed: f64, pub direction: f64,
712 pub vertical_speed: f64,
720}
721
722impl Default for WindConditions {
723 fn default() -> Self {
724 Self {
725 speed: 0.0,
726 direction: 0.0,
727 vertical_speed: 0.0,
728 }
729 }
730}
731
732#[derive(Debug, Clone)]
734pub struct AtmosphericConditions {
735 pub temperature: f64, pub pressure: f64, pub humidity: f64,
741 pub altitude: f64, }
743
744impl Default for AtmosphericConditions {
745 fn default() -> Self {
746 Self {
747 temperature: 15.0,
748 pressure: 1013.25,
749 humidity: 50.0,
750 altitude: 0.0,
751 }
752 }
753}
754
755#[derive(Debug, Clone)]
757pub struct TrajectoryPoint {
758 pub time: f64,
759 pub position: Vector3<f64>,
760 pub velocity_magnitude: f64,
761 pub kinetic_energy: f64,
762 pub drag_coefficient: Option<f64>,
769}
770
771impl TrajectoryPoint {
772 pub fn drag_coefficient_json_value(&self, with_drag_coefficient: bool) -> Option<f64> {
780 if with_drag_coefficient {
781 self.drag_coefficient
782 } else {
783 None
784 }
785 }
786}
787
788#[derive(Debug, Clone)]
790pub struct TrajectoryResult {
791 pub max_range: f64,
792 pub max_height: f64,
793 pub time_of_flight: f64,
794 pub impact_velocity: f64,
795 pub impact_energy: f64,
796 pub projectile_mass_kg: f64,
798 pub line_of_sight_height_m: f64,
800 pub station_speed_of_sound_mps: f64,
802 pub termination: TrajectoryTermination,
804 pub points: Vec<TrajectoryPoint>,
805 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>,
814 pub mach_1_2_distance_m: Option<f64>,
819 pub mach_1_0_distance_m: Option<f64>,
823 pub mach_0_9_distance_m: Option<f64>,
831}
832
833const RK45_TOLERANCE: f64 = 1e-6;
834const RK45_SAFETY_FACTOR: f64 = 0.9;
835const RK45_MAX_DT: f64 = 0.01;
836const RK45_MIN_DT: f64 = 1e-6;
837const TRAJECTORY_TIME_LIMIT_S: f64 = 100.0;
838
839pub const MAX_TRAJECTORY_POINTS: usize = 250_000;
845
846fn cli_rk45_error_norm(
848 position: &Vector3<f64>,
849 velocity: &Vector3<f64>,
850 fifth_position: &Vector3<f64>,
851 fifth_velocity: &Vector3<f64>,
852 fourth_position: &Vector3<f64>,
853 fourth_velocity: &Vector3<f64>,
854) -> f64 {
855 let pack_state = |position: &Vector3<f64>, velocity: &Vector3<f64>| {
856 Vector6::new(
857 position.x, position.y, position.z, velocity.x, velocity.y, velocity.z,
858 )
859 };
860 let state = pack_state(position, velocity);
861 let fifth_order = pack_state(fifth_position, fifth_velocity);
862 let fourth_order = pack_state(fourth_position, fourth_velocity);
863
864 crate::trajectory_integration::rk45_error_norm(&state, &fifth_order, &fourth_order)
865}
866
867struct Rk45Trial {
868 position: Vector3<f64>,
869 velocity: Vector3<f64>,
870 suggested_dt: f64,
871 error: f64,
872}
873
874struct Rk45AcceptedStep {
875 position: Vector3<f64>,
876 velocity: Vector3<f64>,
877 used_dt: f64,
878 next_dt: f64,
879 error: f64,
880}
881
882#[derive(Default)]
896struct MachTransitionTracker {
897 previous_mach: Option<f64>,
898 crossed_transonic: bool,
899 crossed_subsonic: bool,
900 crossed_narrow: bool,
901 mach_1_2_distance_m: Option<f64>,
904 mach_1_0_distance_m: Option<f64>,
907 mach_0_9_distance_m: Option<f64>,
910}
911
912impl MachTransitionTracker {
913 fn record_downward_crossings(&mut self, mach: f64, downrange_m: f64, distances: &mut Vec<f64>) {
914 if !mach.is_finite() {
915 self.previous_mach = None;
916 return;
917 }
918
919 if let Some(previous_mach) = self.previous_mach {
920 if !self.crossed_transonic && previous_mach >= 1.2 && mach < 1.2 {
921 self.crossed_transonic = true;
922 distances.push(downrange_m);
923 self.mach_1_2_distance_m = Some(downrange_m);
924 }
925 if !self.crossed_subsonic && previous_mach >= 1.0 && mach < 1.0 {
926 self.crossed_subsonic = true;
927 distances.push(downrange_m);
928 self.mach_1_0_distance_m = Some(downrange_m);
929 }
930 if !self.crossed_narrow && previous_mach >= 0.9 && mach < 0.9 {
931 self.crossed_narrow = true;
932 self.mach_0_9_distance_m = Some(downrange_m);
934 }
935 }
936 self.previous_mach = Some(mach);
937 }
938}
939
940impl TrajectoryResult {
941 pub fn position_at_range(&self, target_range: f64) -> Option<Vector3<f64>> {
945 if self.points.is_empty() {
946 return None;
947 }
948
949 for i in 0..self.points.len() - 1 {
951 let p1 = &self.points[i];
952 let p2 = &self.points[i + 1];
953
954 if p1.position.x <= target_range && p2.position.x >= target_range {
956 let dx = p2.position.x - p1.position.x;
958 if dx.abs() < 1e-10 {
959 return Some(p1.position);
960 }
961 let t = (target_range - p1.position.x) / dx;
962
963 return Some(Vector3::new(
965 target_range,
966 p1.position.y + t * (p2.position.y - p1.position.y),
967 p1.position.z + t * (p2.position.z - p1.position.z),
968 ));
969 }
970 }
971
972 self.points.last().map(|p| p.position)
974 }
975}
976
977#[derive(Debug, Clone, Copy, PartialEq, Eq)]
979enum StationAtmosphereResolution {
980 LegacyDefaultSentinels,
983 Authoritative,
986}
987
988#[derive(Clone)]
989pub struct TrajectorySolver {
990 inputs: BallisticInputs,
991 wind: WindConditions,
992 atmosphere: AtmosphericConditions,
993 station_atmosphere_resolution: StationAtmosphereResolution,
994 max_range: f64,
995 time_step: f64,
996 max_trajectory_points: usize,
997 cluster_bc: Option<ClusterBCDegradation>,
998 precession_nutation_inertias: (f64, f64),
1000 wind_sock: Option<crate::wind::WindSock>,
1005 atmo_sock: Option<crate::atmosphere::AtmoSock>,
1012}
1013
1014#[derive(Debug, Clone, Copy, PartialEq, Eq)]
1026pub(crate) enum ZeroTargetFrame {
1027 SightLine,
1028 WorldVertical,
1029}
1030
1031#[derive(Debug, Clone, Copy, PartialEq)]
1044pub struct ZeroCrossings {
1045 pub near_m: Option<f64>,
1048 pub far_m: Option<f64>,
1051}
1052
1053impl TrajectorySolver {
1054 pub fn new(
1055 inputs: BallisticInputs,
1056 wind: WindConditions,
1057 atmosphere: AtmosphericConditions,
1058 ) -> Self {
1059 Self::new_with_station_atmosphere_resolution(
1060 inputs,
1061 wind,
1062 atmosphere,
1063 StationAtmosphereResolution::LegacyDefaultSentinels,
1064 )
1065 }
1066
1067 pub fn new_with_resolved_station_atmosphere(
1078 inputs: BallisticInputs,
1079 wind: WindConditions,
1080 atmosphere: AtmosphericConditions,
1081 ) -> Self {
1082 Self::new_with_station_atmosphere_resolution(
1083 inputs,
1084 wind,
1085 atmosphere,
1086 StationAtmosphereResolution::Authoritative,
1087 )
1088 }
1089
1090 fn new_with_station_atmosphere_resolution(
1091 mut inputs: BallisticInputs,
1092 wind: WindConditions,
1093 atmosphere: AtmosphericConditions,
1094 station_atmosphere_resolution: StationAtmosphereResolution,
1095 ) -> Self {
1096 inputs.normalize_for_solve();
1101
1102 let cluster_bc = if inputs.use_cluster_bc {
1104 Some(ClusterBCDegradation::new())
1105 } else {
1106 None
1107 };
1108 let precession_nutation_inertias = projectile_moments_of_inertia(
1109 inputs.bullet_mass,
1110 inputs.bullet_diameter,
1111 inputs.bullet_length,
1112 );
1113
1114 Self {
1115 inputs,
1116 wind,
1117 atmosphere,
1118 station_atmosphere_resolution,
1119 max_range: 1000.0,
1120 time_step: 0.001,
1121 max_trajectory_points: MAX_TRAJECTORY_POINTS,
1122 cluster_bc,
1123 precession_nutation_inertias,
1124 wind_sock: None,
1125 atmo_sock: None,
1126 }
1127 }
1128
1129 pub fn set_max_range(&mut self, range: f64) {
1130 self.max_range = range;
1131 }
1132
1133 pub fn set_time_step(&mut self, step: f64) {
1134 self.time_step = step;
1135 }
1136
1137 pub(crate) fn calculate_and_set_zero_angle(
1141 &mut self,
1142 target_distance_m: f64,
1143 target_height_m: f64,
1144 frame: ZeroTargetFrame,
1145 ) -> Result<f64, BallisticsError> {
1146 let angle = self.find_zero_angle(target_distance_m, target_height_m, frame)?;
1147 let angle = if target_distance_m > 0.0 {
1156 angle + self.inputs.zero_poi_vertical_m / target_distance_m
1157 } else {
1158 angle
1159 };
1160 self.inputs.muzzle_angle = angle;
1161 self.apply_windage_zero_bias(target_distance_m);
1162 Ok(angle)
1163 }
1164
1165 pub(crate) fn apply_windage_zero_bias(&mut self, target_distance_m: f64) {
1179 self.inputs.azimuth_angle += self.inputs.windage_zero_bias_rad(target_distance_m);
1180 }
1181
1182 fn find_zero_angle(
1183 &self,
1184 target_distance_m: f64,
1185 target_height_m: f64,
1186 frame: ZeroTargetFrame,
1187 ) -> Result<f64, BallisticsError> {
1188 let mut low_angle = 0.0;
1191 let mut high_angle = 0.2; let tolerance = 1e-7;
1193 let max_iterations = 60;
1194
1195 let low_height = self.zero_trial_height_at(low_angle, target_distance_m, frame)?;
1197 let high_height = self.zero_trial_height_at(high_angle, target_distance_m, frame)?;
1198
1199 match (low_height, high_height) {
1200 (Some(low_height), Some(high_height)) => {
1201 let low_error = low_height - target_height_m;
1202 let high_error = high_height - target_height_m;
1203
1204 if low_error > 0.0 && high_error > 0.0 {
1205 } else if low_error < 0.0 && high_error < 0.0 {
1208 let mut expanded = false;
1210 for multiplier in [2.0, 3.0, 4.0] {
1211 let new_high = (high_angle * multiplier).min(0.785);
1212 if let Ok(Some(height)) =
1213 self.zero_trial_height_at(new_high, target_distance_m, frame)
1214 {
1215 if height - target_height_m > 0.0 {
1216 high_angle = new_high;
1217 expanded = true;
1218 break;
1219 }
1220 }
1221 if new_high >= 0.785 {
1222 break;
1223 }
1224 }
1225 if !expanded {
1226 return Err("Cannot find zero angle: target beyond effective range even at maximum angle".into());
1227 }
1228 }
1229 }
1230 (None, Some(_)) => {
1231 }
1234 (Some(_), None) => {
1235 return Err(
1236 "Cannot find zero angle: high angle trajectory doesn't reach target distance"
1237 .into(),
1238 );
1239 }
1240 (None, None) => {
1241 return Err(
1242 "Cannot find zero angle: trajectory cannot reach target distance at any angle"
1243 .into(),
1244 );
1245 }
1246 }
1247
1248 for _ in 0..max_iterations {
1249 let mid_angle = (low_angle + high_angle) / 2.0;
1250 match self.zero_trial_height_at(mid_angle, target_distance_m, frame)? {
1251 Some(height) => {
1252 let error = height - target_height_m;
1253 if error.abs() < 0.0001 {
1256 return Ok(mid_angle);
1257 }
1258
1259 if (high_angle - low_angle).abs() < tolerance {
1262 if error.abs() < 0.01 {
1263 return Ok(mid_angle);
1264 }
1265 return Err("Zero angle did not converge: residual height error too large (target not reachable / not bracketed)".into());
1266 }
1267
1268 if error > 0.0 {
1269 high_angle = mid_angle;
1270 } else {
1271 low_angle = mid_angle;
1272 }
1273 }
1274 None => {
1275 low_angle = mid_angle;
1276 if (high_angle - low_angle).abs() < tolerance {
1277 return Err("Trajectory cannot reach target distance - angle converged without valid solution".into());
1278 }
1279 }
1280 }
1281 }
1282
1283 Err("Failed to find zero angle".into())
1284 }
1285
1286 fn zero_trial_height_at(
1289 &self,
1290 angle_rad: f64,
1291 target_distance_m: f64,
1292 frame: ZeroTargetFrame,
1293 ) -> Result<Option<f64>, BallisticsError> {
1294 let mut trial = self.clone();
1295 trial.inputs.muzzle_angle = angle_rad;
1296 trial.inputs.enable_aerodynamic_jump = false;
1299 trial.inputs.cant_angle = 0.0;
1302 if frame == ZeroTargetFrame::SightLine {
1309 trial.inputs.shooting_angle = 0.0;
1310 }
1311 trial.set_max_range(target_distance_m * 2.0);
1312 let result = trial.solve()?;
1313
1314 for (index, point) in result.points.iter().enumerate() {
1315 if point.position.x >= target_distance_m {
1316 let shot_y_m = if index == 0 {
1317 point.position.y
1318 } else {
1319 let previous = &result.points[index - 1];
1320 let span = point.position.x - previous.position.x;
1321 let fraction = (target_distance_m - previous.position.x) / span;
1322 previous.position.y + fraction * (point.position.y - previous.position.y)
1323 };
1324 return Ok(Some(crate::atmosphere::shot_frame_altitude(
1325 0.0,
1326 target_distance_m,
1327 shot_y_m,
1328 trial.inputs.shooting_angle,
1329 )));
1330 }
1331 }
1332 Ok(None)
1333 }
1334
1335 fn find_zero_range(
1362 &self,
1363 angle_rad: f64,
1364 target_height_m: f64,
1365 frame: ZeroTargetFrame,
1366 ) -> Result<ZeroCrossings, BallisticsError> {
1367 let mut trial = self.clone();
1368 trial.inputs.muzzle_angle = angle_rad;
1369 trial.inputs.enable_aerodynamic_jump = false;
1372 trial.inputs.cant_angle = 0.0;
1373 if frame == ZeroTargetFrame::SightLine {
1374 trial.inputs.shooting_angle = 0.0;
1375 }
1376 let result = trial.solve()?;
1377
1378 let mut near_crossing: Option<f64> = None;
1386 let mut far_crossing: Option<f64> = None;
1387 let mut previous: Option<(f64, f64)> = None; for point in &result.points {
1389 let height = crate::atmosphere::shot_frame_altitude(
1390 0.0,
1391 point.position.x,
1392 point.position.y,
1393 trial.inputs.shooting_angle,
1394 );
1395 let error = height - target_height_m;
1396 if let Some((prev_x, prev_error)) = previous {
1397 if prev_error == 0.0 {
1398 if near_crossing.is_none() {
1401 near_crossing = Some(prev_x);
1402 } else {
1403 far_crossing = Some(prev_x);
1404 }
1405 }
1406 if prev_error * error < 0.0 {
1407 let fraction = prev_error / (prev_error - error);
1408 let crossing = prev_x + fraction * (point.position.x - prev_x);
1409 if prev_error < 0.0 && error > 0.0 {
1410 if near_crossing.is_none() {
1412 near_crossing = Some(crossing);
1413 }
1414 } else {
1415 far_crossing = Some(crossing);
1417 }
1418 }
1419 }
1420 previous = Some((point.position.x, error));
1421 }
1422 if let Some((last_x, last_error)) = previous {
1425 if last_error == 0.0 {
1426 if near_crossing.is_none() {
1427 near_crossing = Some(last_x);
1428 } else {
1429 far_crossing = Some(last_x);
1430 }
1431 }
1432 }
1433
1434 if near_crossing.is_none() && far_crossing.is_none() {
1435 return Err(BallisticsError::from(
1436 "Cannot find zero range: this angle never crosses the target height within the \
1437 solved range (angle too shallow to reach it, or both crossings lie beyond \
1438 the solver's max range)."
1439 .to_string(),
1440 ));
1441 }
1442
1443 Ok(ZeroCrossings {
1444 near_m: near_crossing,
1445 far_m: far_crossing,
1446 })
1447 }
1448
1449 pub fn equivalent_horizontal_range(
1473 &self,
1474 target_range_m: f64,
1475 zero_distance_m: f64,
1476 ) -> Option<f64> {
1477 if !target_range_m.is_finite() || !zero_distance_m.is_finite() {
1478 return None;
1479 }
1480 if target_range_m <= zero_distance_m || target_range_m <= 0.0 {
1481 return None;
1482 }
1483
1484 fn path_y_at(points: &[TrajectoryPoint], distance_m: f64) -> Option<f64> {
1489 match bracket_param(points.len(), |i| points[i].position.x, distance_m) {
1490 Bracket::Below => Some(points[0].position.y),
1495 Bracket::Above | Bracket::Degenerate => None,
1496 Bracket::Inside { lo, t } => {
1497 let hi = lo + 1;
1498 Some(
1499 points[lo].position.y + t * (points[hi].position.y - points[lo].position.y),
1500 )
1501 }
1502 }
1503 }
1504
1505 let mut inclined = self.clone();
1509 inclined.inputs.enable_trajectory_sampling = false;
1510 let inclined_result = inclined.solve().ok()?;
1511 let los_height = inclined_result.line_of_sight_height_m;
1512 let inclined_drop = los_height - path_y_at(&inclined_result.points, target_range_m)?;
1513 let correction = inclined_drop / target_range_m;
1514 if correction <= 0.0 {
1515 return None;
1516 }
1517
1518 let mut flat = self.clone();
1520 flat.inputs.enable_trajectory_sampling = false;
1521 flat.inputs.shooting_angle = 0.0;
1522 let flat_result = flat.solve().ok()?;
1523 let flat_correction_at = |range_m: f64| -> Option<f64> {
1524 Some((los_height - path_y_at(&flat_result.points, range_m)?) / range_m)
1525 };
1526
1527 let flat_terminal_m = flat_result.points.last().map(|p| p.position.x)?;
1535 let mut low = zero_distance_m.max(1.0);
1536 let mut high = target_range_m.min(flat_terminal_m);
1537 if high <= low {
1538 return None;
1539 }
1540 if flat_correction_at(low)? - correction > 0.0 {
1541 return None; }
1543 if flat_correction_at(high)? - correction < 0.0 {
1544 return None; }
1546 for _ in 0..60 {
1547 let mid = 0.5 * (low + high);
1548 let error = flat_correction_at(mid)? - correction;
1549 if error.abs() == 0.0 {
1550 return Some(mid);
1551 }
1552 if error < 0.0 {
1553 low = mid;
1554 } else {
1555 high = mid;
1556 }
1557 if high - low < 0.01 {
1558 break;
1559 }
1560 }
1561 Some(0.5 * (low + high))
1562 }
1563
1564 fn validate_for_solve(&self) -> Result<(), BallisticsError> {
1570 let require_finite = |name: &str, value: f64| {
1571 if value.is_finite() {
1572 Ok(())
1573 } else {
1574 Err(BallisticsError::from(format!("{name} must be finite")))
1575 }
1576 };
1577 let require_positive = |name: &str, value: f64| {
1578 if value.is_finite() && value > 0.0 {
1579 Ok(())
1580 } else {
1581 Err(BallisticsError::from(format!(
1582 "{name} must be finite and greater than zero"
1583 )))
1584 }
1585 };
1586
1587 if self.inputs.custom_drag_table.is_none() {
1593 require_positive("bc_value", self.inputs.bc_value)?;
1594 }
1595 require_positive("bullet_mass", self.inputs.bullet_mass)?;
1596 require_positive("bullet_diameter", self.inputs.bullet_diameter)?;
1597 require_positive("muzzle_velocity", self.inputs.muzzle_velocity)?;
1598 require_positive("cd_scale", self.inputs.cd_scale)?;
1604
1605 require_finite("muzzle_angle", self.inputs.muzzle_angle)?;
1606 require_finite("azimuth_angle", self.inputs.azimuth_angle)?;
1607 require_finite("shooting_angle", self.inputs.shooting_angle)?;
1608 require_finite("cant_angle", self.inputs.cant_angle)?;
1609 require_finite("muzzle_height", self.inputs.muzzle_height)?;
1610
1611 for (name, value) in [
1615 ("zero_poi_vertical_m", self.inputs.zero_poi_vertical_m),
1616 ("zero_poi_horizontal_m", self.inputs.zero_poi_horizontal_m),
1617 ] {
1618 require_finite(name, value)?;
1619 if value.abs() >= 1.0 {
1620 return Err(BallisticsError::from(format!(
1621 "{name} must be smaller than 1.0 m in magnitude (it is a linear POI \
1622 offset at the zero range, in meters)"
1623 )));
1624 }
1625 }
1626
1627 require_finite("sight_offset_lateral_m", self.inputs.sight_offset_lateral_m)?;
1631 if self.inputs.sight_offset_lateral_m.abs() >= 0.5 {
1632 return Err(BallisticsError::from(
1633 "sight_offset_lateral_m must be smaller than 0.5 m in magnitude (it is \
1634 the lateral sight-to-bore mount offset, in meters)",
1635 ));
1636 }
1637
1638 if !(self.inputs.ground_threshold.is_finite()
1641 || self.inputs.ground_threshold == f64::NEG_INFINITY)
1642 {
1643 return Err(BallisticsError::from(
1644 "ground_threshold must be finite or negative infinity",
1645 ));
1646 }
1647
1648 match &self.wind_sock {
1649 Some(wind_sock) => wind_sock
1650 .validate_segments()
1651 .map_err(BallisticsError::from)?,
1652 None => {
1653 require_finite("wind.speed", self.wind.speed)?;
1654 require_finite("wind.direction", self.wind.direction)?;
1655 require_finite("wind.vertical_speed", self.wind.vertical_speed)?;
1656 }
1657 }
1658
1659 require_finite("atmosphere.temperature", self.atmosphere.temperature)?;
1660 require_finite("atmosphere.pressure", self.atmosphere.pressure)?;
1661 require_finite("atmosphere.humidity", self.atmosphere.humidity)?;
1662 require_finite("atmosphere.altitude", self.atmosphere.altitude)?;
1663
1664 require_positive("max_range", self.max_range)?;
1665 if !self.inputs.use_rk4 || !self.inputs.use_adaptive_rk45 {
1668 require_positive("time_step", self.time_step)?;
1669 }
1670
1671 if self.inputs.enable_trajectory_sampling {
1672 require_finite("sight_height", self.inputs.sight_height)?;
1673 require_positive("sample_interval", self.inputs.sample_interval)?;
1674 projected_sample_count(self.max_range, self.inputs.sample_interval)?;
1675 }
1676
1677 if self.inputs.drops_reference == DropsReference::Target {
1681 require_finite("target_height", self.inputs.target_height)?;
1682 if self.inputs.shooting_angle.cos() <= 1e-9 {
1683 return Err(BallisticsError::from(
1684 "drops reference 'target' is undefined for shooting angles at or beyond 90 degrees",
1685 ));
1686 }
1687 }
1688
1689 if self.inputs.enable_coriolis {
1690 require_finite("shot_azimuth", self.inputs.shot_azimuth)?;
1691 if let Some(latitude) = self.inputs.latitude {
1692 require_finite("latitude", latitude)?;
1693 }
1694 }
1695
1696 Ok(())
1697 }
1698
1699 fn validate_result_sanity(&self, result: &TrajectoryResult) -> Result<(), BallisticsError> {
1706 let require_finite = |name: &str, value: f64| {
1707 if value.is_finite() {
1708 Ok(())
1709 } else {
1710 Err(BallisticsError::from(format!(
1711 "trajectory result contains non-finite {name}"
1712 )))
1713 }
1714 };
1715 let require_non_negative = |name: &str, value: f64| {
1716 if value >= 0.0 {
1717 Ok(())
1718 } else {
1719 Err(BallisticsError::from(format!(
1720 "trajectory result contains non-physical negative {name} ({value})"
1721 )))
1722 }
1723 };
1724 let require_indexed_finite = |collection: &str, index: usize, field: &str, value: f64| {
1725 if value.is_finite() {
1726 Ok(())
1727 } else {
1728 Err(BallisticsError::from(format!(
1729 "trajectory result contains non-finite {collection}[{index}].{field}"
1730 )))
1731 }
1732 };
1733 let require_indexed_non_negative = |collection: &str,
1734 index: usize,
1735 field: &str,
1736 value: f64| {
1737 if value >= 0.0 {
1738 Ok(())
1739 } else {
1740 Err(BallisticsError::from(format!(
1741 "trajectory result contains non-physical negative {collection}[{index}].{field} ({value})"
1742 )))
1743 }
1744 };
1745
1746 require_finite("max_range", result.max_range)?;
1747 require_finite("max_height", result.max_height)?;
1748 require_finite("time_of_flight", result.time_of_flight)?;
1749 require_finite("impact_velocity", result.impact_velocity)?;
1750 require_finite("impact_energy", result.impact_energy)?;
1751 require_finite("projectile_mass_kg", result.projectile_mass_kg)?;
1752 require_finite("line_of_sight_height_m", result.line_of_sight_height_m)?;
1753 require_finite(
1754 "station_speed_of_sound_mps",
1755 result.station_speed_of_sound_mps,
1756 )?;
1757
1758 require_non_negative("max_range", result.max_range)?;
1762 require_non_negative("time_of_flight", result.time_of_flight)?;
1763 require_non_negative("impact_velocity", result.impact_velocity)?;
1764 require_non_negative("impact_energy", result.impact_energy)?;
1765 require_non_negative("projectile_mass_kg", result.projectile_mass_kg)?;
1766 require_non_negative(
1767 "station_speed_of_sound_mps",
1768 result.station_speed_of_sound_mps,
1769 )?;
1770
1771 for (index, point) in result.points.iter().enumerate() {
1772 require_indexed_finite("points", index, "time", point.time)?;
1773 require_indexed_finite("points", index, "position.x", point.position.x)?;
1774 require_indexed_finite("points", index, "position.y", point.position.y)?;
1775 require_indexed_finite("points", index, "position.z", point.position.z)?;
1776 require_indexed_finite(
1777 "points",
1778 index,
1779 "velocity_magnitude",
1780 point.velocity_magnitude,
1781 )?;
1782 require_indexed_finite("points", index, "kinetic_energy", point.kinetic_energy)?;
1783 require_indexed_non_negative("points", index, "time", point.time)?;
1784 require_indexed_non_negative(
1785 "points",
1786 index,
1787 "velocity_magnitude",
1788 point.velocity_magnitude,
1789 )?;
1790 require_indexed_non_negative("points", index, "kinetic_energy", point.kinetic_energy)?;
1791 }
1792
1793 if let Some(samples) = &result.sampled_points {
1794 for (index, sample) in samples.iter().enumerate() {
1795 require_indexed_finite("sampled_points", index, "distance_m", sample.distance_m)?;
1796 require_indexed_finite("sampled_points", index, "drop_m", sample.drop_m)?;
1797 require_indexed_finite(
1798 "sampled_points",
1799 index,
1800 "wind_drift_m",
1801 sample.wind_drift_m,
1802 )?;
1803 require_indexed_finite(
1804 "sampled_points",
1805 index,
1806 "velocity_mps",
1807 sample.velocity_mps,
1808 )?;
1809 require_indexed_finite("sampled_points", index, "energy_j", sample.energy_j)?;
1810 require_indexed_finite("sampled_points", index, "time_s", sample.time_s)?;
1811 }
1812 }
1813
1814 for (name, value) in [
1815 ("min_pitch_damping", result.min_pitch_damping),
1816 ("transonic_mach", result.transonic_mach),
1817 ("max_yaw_angle", result.max_yaw_angle),
1818 ("max_precession_angle", result.max_precession_angle),
1819 ] {
1820 if let Some(value) = value {
1821 require_finite(name, value)?;
1822 }
1823 }
1824
1825 if let Some(state) = result.angular_state {
1826 for (name, value) in [
1827 ("angular_state.pitch_angle", state.pitch_angle),
1828 ("angular_state.yaw_angle", state.yaw_angle),
1829 ("angular_state.pitch_rate", state.pitch_rate),
1830 ("angular_state.yaw_rate", state.yaw_rate),
1831 ("angular_state.precession_angle", state.precession_angle),
1832 ("angular_state.nutation_phase", state.nutation_phase),
1833 ] {
1834 require_finite(name, value)?;
1835 }
1836 }
1837
1838 if let Some(jump) = result.aerodynamic_jump {
1839 for (name, value) in [
1840 ("aerodynamic_jump.vertical_jump_moa", jump.vertical_jump_moa),
1841 (
1842 "aerodynamic_jump.horizontal_jump_moa",
1843 jump.horizontal_jump_moa,
1844 ),
1845 ("aerodynamic_jump.jump_angle_rad", jump.jump_angle_rad),
1846 (
1847 "aerodynamic_jump.magnus_component_moa",
1848 jump.magnus_component_moa,
1849 ),
1850 ("aerodynamic_jump.yaw_component_moa", jump.yaw_component_moa),
1851 (
1852 "aerodynamic_jump.stabilization_factor",
1853 jump.stabilization_factor,
1854 ),
1855 ] {
1856 require_finite(name, value)?;
1857 }
1858 }
1859
1860 Ok(())
1861 }
1862
1863 fn validate_integration_state(
1876 &self,
1877 position: &Vector3<f64>,
1878 velocity: &Vector3<f64>,
1879 time: f64,
1880 ) -> Result<(), BallisticsError> {
1881 if !(position.iter().all(|value| value.is_finite())
1882 && velocity.iter().all(|value| value.is_finite())
1883 && time.is_finite())
1884 {
1885 return Err(BallisticsError::from(
1886 "trajectory integration produced a non-finite state (often from physically \
1887 extreme inputs — e.g. an absurd bore/muzzle height placing the launch far \
1888 from sea level, or a degenerate atmosphere; check those inputs, or set \
1889 --altitude explicitly)",
1890 ));
1891 }
1892
1893 let speed = velocity.magnitude();
1894 let budget = self.speed_budget(time);
1895 if speed > budget {
1896 return Err(BallisticsError::from(format!(
1897 "trajectory integration diverged: speed {speed:.3e} m/s at t={time:.6}s exceeds \
1898 the physical budget of {budget:.3e} m/s"
1899 )));
1900 }
1901 Ok(())
1902 }
1903
1904 fn speed_budget(&self, time: f64) -> f64 {
1909 let scalar_wind = self.wind.speed.abs() + self.wind.vertical_speed.abs();
1910 let wind_bound = match &self.wind_sock {
1911 Some(sock) => scalar_wind.max(sock.max_speed_mps()),
1912 None => scalar_wind,
1913 };
1914 2.0 * (self.inputs.muzzle_velocity + wind_bound + 10.0)
1915 + crate::constants::G_ACCEL_MPS2 * time
1916 }
1917
1918 fn push_trajectory_point(
1922 &self,
1923 points: &mut Vec<TrajectoryPoint>,
1924 point: TrajectoryPoint,
1925 collapsed_steps: &mut usize,
1926 ) -> Result<(), BallisticsError> {
1927 if points.len() + *collapsed_steps >= self.max_trajectory_points {
1932 return Err(BallisticsError::from(format!(
1933 "trajectory point limit of {} exceeded",
1934 self.max_trajectory_points
1935 )));
1936 }
1937 if let Some(last) = points.last_mut() {
1948 if point.position.x == last.position.x {
1949 *last = point;
1950 *collapsed_steps += 1;
1951 return Ok(());
1952 }
1953 }
1954 points.push(point);
1955 Ok(())
1956 }
1957
1958 pub fn set_wind_segments(&mut self, segments: Vec<crate::wind::WindSegment>) {
1965 self.wind_sock = if segments.is_empty() {
1966 None
1967 } else {
1968 Some(crate::wind::WindSock::new(segments))
1969 };
1970 }
1971
1972 pub fn set_atmo_segments(&mut self, segments: Vec<crate::atmosphere::AtmoSegment>) {
1981 self.atmo_sock = if segments.is_empty() {
1982 None
1983 } else {
1984 Some(crate::atmosphere::AtmoSock::new(segments))
1985 };
1986 }
1987
1988 fn launch_angles_from(
1997 &self,
1998 aj: Option<&crate::aerodynamic_jump::AerodynamicJumpComponents>,
1999 ) -> (f64, f64) {
2000 let (mut elev, mut azim) = (self.inputs.muzzle_angle, self.inputs.azimuth_angle);
2001 if self.inputs.cant_angle != 0.0 {
2008 let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
2009 let (e0, a0) = (elev, azim);
2010 elev = e0 * cos_c - a0 * sin_c;
2011 azim = a0 * cos_c + e0 * sin_c;
2012 }
2013 match aj {
2014 Some(c) => {
2015 const MOA_PER_RAD: f64 = 3437.7467707849;
2017 (
2018 elev + c.vertical_jump_moa / MOA_PER_RAD,
2019 azim + c.horizontal_jump_moa / MOA_PER_RAD,
2020 )
2021 }
2022 None => (elev, azim),
2023 }
2024 }
2025
2026 fn aerodynamic_jump_components(
2034 &self,
2035 ) -> Option<crate::aerodynamic_jump::AerodynamicJumpComponents> {
2036 if !self.inputs.enable_aerodynamic_jump {
2037 return None;
2038 }
2039 let diameter_m = self.inputs.bullet_diameter;
2043 if !(self.inputs.twist_rate.is_finite()
2044 && self.inputs.twist_rate != 0.0
2045 && diameter_m.is_finite()
2046 && diameter_m > 0.0
2047 && self.inputs.bullet_length.is_finite()
2048 && self.inputs.bullet_length > 0.0
2049 && self.inputs.muzzle_velocity.is_finite())
2050 {
2051 return None;
2052 }
2053
2054 let (_, _, temp_c, pressure_hpa) = self.resolved_atmosphere();
2056 let sg = crate::stability::compute_stability_coefficient(
2057 &self.inputs,
2058 (self.atmosphere.altitude, temp_c, pressure_hpa, 0.0),
2059 );
2060 if !(sg.is_finite() && sg > 0.0) {
2061 return None;
2062 }
2063 let length_calibers = self.inputs.bullet_length / diameter_m;
2064
2065 const MS_TO_MPH: f64 = 2.236_936_292_054_4;
2071 let crosswind_from_right_mps = if let Some(sock) = &self.wind_sock {
2072 -sock.vector_for_range_stateless(0.0)[2]
2073 } else {
2074 self.wind.speed * self.wind.direction.sin()
2075 };
2076 let crosswind_from_right_mph = crosswind_from_right_mps * MS_TO_MPH;
2077
2078 let vertical_jump_moa = crate::aerodynamic_jump::litz_crosswind_jump_moa(
2079 sg,
2080 length_calibers,
2081 crosswind_from_right_mph,
2082 self.inputs.is_twist_right,
2083 );
2084 if !vertical_jump_moa.is_finite() {
2085 return None;
2086 }
2087
2088 const MOA_PER_RAD: f64 = 3437.7467707849;
2089 Some(crate::aerodynamic_jump::AerodynamicJumpComponents {
2090 vertical_jump_moa,
2091 horizontal_jump_moa: 0.0,
2093 jump_angle_rad: vertical_jump_moa.abs() / MOA_PER_RAD,
2094 magnus_component_moa: 0.0,
2095 yaw_component_moa: 0.0,
2096 stabilization_factor: (sg / 1.5).clamp(0.0, 1.0),
2097 })
2098 }
2099
2100 fn resolved_atmosphere(&self) -> (f64, f64, f64, f64) {
2101 let (temp_c, pressure_hpa) = match self.station_atmosphere_resolution {
2102 StationAtmosphereResolution::LegacyDefaultSentinels => {
2103 crate::atmosphere::resolve_station_conditions(
2104 self.atmosphere.temperature,
2105 self.atmosphere.pressure,
2106 self.atmosphere.altitude,
2107 )
2108 }
2109 StationAtmosphereResolution::Authoritative => {
2110 (self.atmosphere.temperature, self.atmosphere.pressure)
2111 }
2112 };
2113 let (density, speed_of_sound) = crate::atmosphere::calculate_atmosphere(
2114 self.atmosphere.altitude,
2115 Some(temp_c),
2116 Some(pressure_hpa),
2117 self.atmosphere.humidity,
2118 );
2119 (density, speed_of_sound, temp_c, pressure_hpa)
2120 }
2121
2122 fn precession_nutation_params(
2123 &self,
2124 velocity_mps: f64,
2125 air_density_kg_m3: f64,
2126 speed_of_sound_mps: f64,
2127 ) -> PrecessionNutationParams {
2128 let (spin_inertia, transverse_inertia) = self.precession_nutation_inertias;
2129 let spin_rate_rad_s = if self.inputs.twist_rate > 0.0 {
2130 let velocity_fps = velocity_mps * 3.28084;
2131 let twist_rate_ft = self.inputs.twist_rate / 12.0;
2132 (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI
2133 } else {
2134 0.0
2135 };
2136
2137 PrecessionNutationParams {
2138 mass_kg: self.inputs.bullet_mass,
2139 caliber_m: self.inputs.bullet_diameter,
2140 length_m: self.inputs.bullet_length,
2141 spin_rate_rad_s,
2142 spin_inertia,
2143 transverse_inertia,
2144 velocity_mps,
2145 air_density_kg_m3,
2146 mach: velocity_mps / speed_of_sound_mps,
2147 pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
2148 nutation_damping_factor: 0.05,
2149 }
2150 }
2151
2152 fn append_terminal_endpoint(
2161 &self,
2162 points: &mut Vec<TrajectoryPoint>,
2163 post_position: Vector3<f64>,
2164 post_velocity: Vector3<f64>,
2165 post_time: f64,
2166 max_height: &mut f64,
2167 collapsed_steps: &mut usize,
2168 ) -> Result<TrajectoryTermination, BallisticsError> {
2169 let previous = points
2170 .last()
2171 .cloned()
2172 .ok_or_else(|| BallisticsError::from("No trajectory points generated"))?;
2173
2174 let mut crossings = Vec::with_capacity(3);
2175 if previous.position.x < self.max_range && post_position.x >= self.max_range {
2176 let span = post_position.x - previous.position.x;
2177 if span.is_finite() && span > 0.0 {
2178 crossings.push((
2179 (self.max_range - previous.position.x) / span,
2180 TrajectoryTermination::MaxRange,
2181 ));
2182 }
2183 }
2184 if self.inputs.ground_threshold.is_finite()
2185 && previous.position.y > self.inputs.ground_threshold
2186 && post_position.y <= self.inputs.ground_threshold
2187 {
2188 let span = post_position.y - previous.position.y;
2189 if span.is_finite() && span < 0.0 {
2190 crossings.push((
2191 (self.inputs.ground_threshold - previous.position.y) / span,
2192 TrajectoryTermination::GroundThreshold,
2193 ));
2194 }
2195 }
2196 if previous.time < TRAJECTORY_TIME_LIMIT_S && post_time >= TRAJECTORY_TIME_LIMIT_S {
2197 let span = post_time - previous.time;
2198 if span.is_finite() && span > 0.0 {
2199 crossings.push((
2200 (TRAJECTORY_TIME_LIMIT_S - previous.time) / span,
2201 TrajectoryTermination::TimeLimit,
2202 ));
2203 }
2204 }
2205
2206 let (fraction, termination) = crossings
2207 .into_iter()
2208 .filter(|(fraction, _)| fraction.is_finite() && (0.0..=1.0).contains(fraction))
2209 .min_by(|left, right| {
2210 let priority = |termination: TrajectoryTermination| match termination {
2211 TrajectoryTermination::GroundThreshold => 0,
2212 TrajectoryTermination::MaxRange => 1,
2213 TrajectoryTermination::TimeLimit => 2,
2214 TrajectoryTermination::VelocityFloor => 3,
2215 };
2216 left.0
2217 .total_cmp(&right.0)
2218 .then_with(|| priority(left.1).cmp(&priority(right.1)))
2219 })
2220 .ok_or_else(|| {
2221 BallisticsError::from(
2222 "trajectory integration stopped without crossing a supported boundary",
2223 )
2224 })?;
2225
2226 let mut position = previous.position + (post_position - previous.position) * fraction;
2227 match termination {
2228 TrajectoryTermination::MaxRange => position.x = self.max_range,
2229 TrajectoryTermination::GroundThreshold => {
2230 position.y = self.inputs.ground_threshold;
2231 }
2232 TrajectoryTermination::TimeLimit | TrajectoryTermination::VelocityFloor => {}
2233 }
2234 let velocity_magnitude = previous.velocity_magnitude
2235 + (post_velocity.magnitude() - previous.velocity_magnitude) * fraction;
2236 let mut time = previous.time + (post_time - previous.time) * fraction;
2237 if termination == TrajectoryTermination::TimeLimit {
2238 time = TRAJECTORY_TIME_LIMIT_S;
2239 }
2240 let kinetic_energy =
2241 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2242
2243 if position.y > *max_height {
2244 *max_height = position.y;
2245 }
2246 let terminal_point = TrajectoryPoint {
2247 time,
2248 position,
2249 velocity_magnitude,
2250 kinetic_energy,
2251 drag_coefficient: None,
2252 };
2253 if terminal_point.position.x < previous.position.x {
2254 return Err(BallisticsError::from(
2255 "trajectory terminal state reversed downrange before the crossed boundary",
2256 ));
2257 }
2258 if terminal_point.position.x == previous.position.x {
2259 let last = points.last_mut().ok_or_else(|| {
2264 BallisticsError::from("trajectory points disappeared during terminal finalization")
2265 })?;
2266 *last = terminal_point;
2267 } else {
2268 self.push_trajectory_point(points, terminal_point, collapsed_steps)?;
2269 }
2270 Ok(termination)
2271 }
2272
2273 fn gravity_acceleration(&self) -> Vector3<f64> {
2274 let theta = self.inputs.shooting_angle;
2275 Vector3::new(
2276 -crate::constants::G_ACCEL_MPS2 * theta.sin(),
2277 -crate::constants::G_ACCEL_MPS2 * theta.cos(),
2278 0.0,
2279 )
2280 }
2281
2282 fn get_wind_at_altitude(&self, altitude_m: f64) -> Vector3<f64> {
2283 let model = match self.inputs.wind_shear_model.as_str() {
2298 "logarithmic" => WindShearModel::Logarithmic,
2299 "power_law" | "powerlaw" | "exponential" => WindShearModel::PowerLaw,
2300 "ekman_spiral" | "ekman" => WindShearModel::EkmanSpiral,
2301 "custom_layers" | "custom" => WindShearModel::CustomLayers,
2302 _ => WindShearModel::PowerLaw,
2303 };
2304 let speed_ratio = crate::wind_shear::boundary_layer_speed_ratio(altitude_m, model);
2305
2306 crate::wind::wind_vector(self.wind.speed, self.wind.direction, 0.0) * speed_ratio
2313 + Vector3::new(0.0, self.wind.vertical_speed, 0.0)
2314 }
2315
2316 pub fn solve(&self) -> Result<TrajectoryResult, BallisticsError> {
2317 self.validate_for_solve()?;
2318 let mut result = if self.inputs.use_rk4 {
2319 if self.inputs.use_adaptive_rk45 {
2320 self.solve_rk45()?
2321 } else {
2322 self.solve_rk4()?
2323 }
2324 } else {
2325 self.solve_euler()?
2326 };
2327 self.apply_spin_drift(&mut result);
2328 self.validate_result_sanity(&result)?;
2329 Ok(result)
2330 }
2331
2332 fn apply_spin_drift(&self, result: &mut TrajectoryResult) {
2338 if !self.inputs.use_enhanced_spin_drift {
2339 return;
2340 }
2341 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 {
2345 return;
2346 }
2347
2348 let sg = self.effective_spin_drift_sg();
2355
2356 for p in result.points.iter_mut() {
2357 if p.time <= 0.0 {
2358 continue;
2359 }
2360 p.position.z +=
2362 crate::spin_drift::litz_drift_meters(sg, p.time, self.inputs.is_twist_right);
2363 }
2364
2365 if let Some(samples) = result.sampled_points.as_mut() {
2369 for s in samples.iter_mut() {
2370 if s.time_s <= 0.0 {
2371 continue;
2372 }
2373 s.wind_drift_m +=
2374 crate::spin_drift::litz_drift_meters(sg, s.time_s, self.inputs.is_twist_right);
2375 }
2376 }
2377 }
2378
2379 fn effective_spin_drift_sg(&self) -> f64 {
2384 let (_, _, temp_c, press_hpa) = self.resolved_atmosphere();
2385 crate::spin_drift::effective_sg_from_inputs(&self.inputs, temp_c, press_hpa)
2386 }
2387
2388 fn initial_position(&self) -> Vector3<f64> {
2406 if self.inputs.cant_angle == 0.0 && self.inputs.sight_offset_lateral_m == 0.0 {
2407 return Vector3::new(0.0, self.inputs.muzzle_height, 0.0);
2408 }
2409 let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
2410 let sh = self.inputs.sight_height;
2411 let off = self.inputs.sight_offset_lateral_m;
2412 Vector3::new(
2413 0.0,
2414 self.inputs.muzzle_height + sh * (1.0 - cos_c) + off * sin_c,
2415 -sh * sin_c - off * cos_c,
2416 )
2417 }
2418
2419 fn build_sampled_points(
2426 &self,
2427 points: &[TrajectoryPoint],
2428 max_height: f64,
2429 transonic_distances: Vec<f64>,
2430 mach_transitions: &MachTransitionTracker,
2431 ) -> Result<Option<Vec<TrajectorySample>>, BallisticsError> {
2432 if !self.inputs.enable_trajectory_sampling {
2433 return Ok(None);
2434 }
2435
2436 let last_point = points.last().ok_or("No trajectory points generated")?;
2437 let trajectory_data = TrajectoryData {
2438 times: points.iter().map(|p| p.time).collect(),
2439 positions: points.iter().map(|p| p.position).collect(),
2440 velocities: points
2441 .iter()
2442 .map(|p| {
2443 Vector3::new(0.0, 0.0, p.velocity_magnitude)
2445 })
2446 .collect(),
2447 transonic_distances, mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2449 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2450 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2451 };
2452
2453 let sight_position_m = self.inputs.muzzle_height + self.inputs.sight_height;
2458 let target_reference = self.inputs.drops_reference == DropsReference::Target;
2459 let target_vertical_height_m = if target_reference && self.inputs.target_height != 0.0 {
2466 self.inputs.target_height
2467 } else {
2468 sight_position_m
2469 };
2470 let outputs = TrajectoryOutputs {
2471 target_distance_horiz_m: last_point.position.x, target_vertical_height_m,
2473 time_of_flight_s: last_point.time,
2474 max_ord_dist_horiz_m: max_height,
2475 sight_height_m: sight_position_m,
2476 };
2477
2478 let mut samples = sample_trajectory(
2480 &trajectory_data,
2481 &outputs,
2482 self.inputs.sample_interval,
2483 self.inputs.bullet_mass,
2484 )?;
2485
2486 if target_reference {
2493 let cos_theta = self.inputs.shooting_angle.cos();
2494 for sample in &mut samples {
2495 sample.drop_m /= cos_theta;
2496 }
2497 }
2498 Ok(Some(samples))
2499 }
2500
2501 fn solve_euler(&self) -> Result<TrajectoryResult, BallisticsError> {
2502 let mut time = 0.0;
2504 let mut position = self.initial_position();
2508 let aj_components = self.aerodynamic_jump_components();
2514 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2515 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2516 let mut velocity = Vector3::new(
2517 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2521
2522 let mut points = Vec::new();
2523 let mut collapsed_steps = 0;
2526 let mut max_height = position.y;
2527 let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
2533 let mut mach_transitions = MachTransitionTracker::default();
2534
2535 let mut angular_state = if self.inputs.enable_precession_nutation {
2537 Some(AngularState {
2538 pitch_angle: 0.001, yaw_angle: 0.001,
2540 pitch_rate: 0.0,
2541 yaw_rate: 0.0,
2542 precession_angle: 0.0,
2543 nutation_phase: 0.0,
2544 })
2545 } else {
2546 None
2547 };
2548 let mut max_yaw_angle = 0.0;
2549 let mut max_precession_angle = 0.0;
2550
2551 let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2553 self.resolved_atmosphere();
2554 let base_ratio = air_density / 1.225;
2559
2560 let wind_vector = crate::wind::wind_vector(
2566 self.wind.speed,
2567 self.wind.direction,
2568 self.wind.vertical_speed,
2569 );
2570
2571 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2574 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2575 );
2576
2577 while position.x < self.max_range
2579 && position.y > self.inputs.ground_threshold
2580 && time < TRAJECTORY_TIME_LIMIT_S
2581 {
2582 let velocity_magnitude = velocity.magnitude();
2584 let kinetic_energy =
2585 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2586
2587 self.push_trajectory_point(
2588 &mut points,
2589 TrajectoryPoint {
2590 time,
2591 position,
2592 velocity_magnitude,
2593 kinetic_energy,
2594 drag_coefficient: None,
2595 },
2596 &mut collapsed_steps,
2597 )?;
2598
2599 {
2602 let mach_here = if speed_of_sound > 0.0 {
2603 velocity_magnitude / speed_of_sound
2604 } else {
2605 0.0
2606 };
2607 mach_transitions.record_downward_crossings(
2608 mach_here,
2609 position.x,
2610 &mut transonic_distances,
2611 );
2612 }
2613
2614 if position.y > max_height {
2616 max_height = position.y;
2617 }
2618
2619 if self.inputs.enable_pitch_damping {
2621 let mach = velocity_magnitude / speed_of_sound;
2622
2623 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
2625 transonic_mach = Some(mach);
2626 }
2627
2628 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
2630
2631 if pitch_damping < min_pitch_damping {
2633 min_pitch_damping = pitch_damping;
2634 }
2635 }
2636
2637 if self.inputs.enable_precession_nutation {
2639 if let Some(ref mut state) = angular_state {
2640 let velocity_magnitude = velocity.magnitude();
2641 let params = self.precession_nutation_params(
2642 velocity_magnitude,
2643 air_density,
2644 speed_of_sound,
2645 );
2646
2647 *state = calculate_combined_angular_motion(
2649 ¶ms,
2650 state,
2651 time,
2652 self.time_step,
2653 0.001, );
2655
2656 if state.yaw_angle.abs() > max_yaw_angle {
2658 max_yaw_angle = state.yaw_angle.abs();
2659 }
2660 if state.precession_angle.abs() > max_precession_angle {
2661 max_precession_angle = state.precession_angle.abs();
2662 }
2663 }
2664 }
2665
2666 let acceleration = self.calculate_acceleration(
2673 &position,
2674 &velocity,
2675 &wind_vector,
2676 (resolved_temp_c, resolved_press_hpa, base_ratio),
2677 );
2678
2679 velocity += acceleration * self.time_step;
2681 position += velocity * self.time_step;
2682 time += self.time_step;
2683 self.validate_integration_state(&position, &velocity, time)?;
2684 }
2685
2686 let termination = self.append_terminal_endpoint(
2687 &mut points,
2688 position,
2689 velocity,
2690 time,
2691 &mut max_height,
2692 &mut collapsed_steps,
2693 )?;
2694
2695 self.annotate_drag_coefficients(&mut points, speed_of_sound);
2700
2701 let last_point = points.last().ok_or("No trajectory points generated")?;
2702
2703 let sampled_points =
2705 self.build_sampled_points(&points, max_height, transonic_distances, &mach_transitions)?;
2706
2707 Ok(TrajectoryResult {
2708 max_range: last_point.position.x, max_height,
2710 time_of_flight: last_point.time,
2711 impact_velocity: last_point.velocity_magnitude,
2712 impact_energy: last_point.kinetic_energy,
2713 projectile_mass_kg: self.inputs.bullet_mass,
2714 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
2715 station_speed_of_sound_mps: speed_of_sound,
2716 termination,
2717 points,
2718 sampled_points,
2719 min_pitch_damping: if self.inputs.enable_pitch_damping {
2720 Some(min_pitch_damping)
2721 } else {
2722 None
2723 },
2724 transonic_mach,
2725 angular_state,
2726 max_yaw_angle: if self.inputs.enable_precession_nutation {
2727 Some(max_yaw_angle)
2728 } else {
2729 None
2730 },
2731 max_precession_angle: if self.inputs.enable_precession_nutation {
2732 Some(max_precession_angle)
2733 } else {
2734 None
2735 },
2736 aerodynamic_jump: aj_components,
2737 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2738 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2739 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2740 })
2741 }
2742
2743 fn solve_rk4(&self) -> Result<TrajectoryResult, BallisticsError> {
2744 let mut time = 0.0;
2746 let mut position = self.initial_position();
2751
2752 let aj_components = self.aerodynamic_jump_components();
2758 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2759 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2760 let mut velocity = Vector3::new(
2761 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2765
2766 let mut points = Vec::new();
2767 let mut collapsed_steps = 0;
2770 let mut max_height = position.y;
2771 let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
2777 let mut mach_transitions = MachTransitionTracker::default();
2778
2779 let mut angular_state = if self.inputs.enable_precession_nutation {
2781 Some(AngularState {
2782 pitch_angle: 0.001, yaw_angle: 0.001,
2784 pitch_rate: 0.0,
2785 yaw_rate: 0.0,
2786 precession_angle: 0.0,
2787 nutation_phase: 0.0,
2788 })
2789 } else {
2790 None
2791 };
2792 let mut max_yaw_angle = 0.0;
2793 let mut max_precession_angle = 0.0;
2794
2795 let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2797 self.resolved_atmosphere();
2798 let base_ratio = air_density / 1.225;
2803
2804 let wind_vector = crate::wind::wind_vector(
2810 self.wind.speed,
2811 self.wind.direction,
2812 self.wind.vertical_speed,
2813 );
2814
2815 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2818 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2819 );
2820
2821 while position.x < self.max_range
2823 && position.y > self.inputs.ground_threshold
2824 && time < TRAJECTORY_TIME_LIMIT_S
2825 {
2826 let velocity_magnitude = velocity.magnitude();
2828 let kinetic_energy =
2829 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2830
2831 self.push_trajectory_point(
2832 &mut points,
2833 TrajectoryPoint {
2834 time,
2835 position,
2836 velocity_magnitude,
2837 kinetic_energy,
2838 drag_coefficient: None,
2839 },
2840 &mut collapsed_steps,
2841 )?;
2842
2843 {
2846 let mach_here = if speed_of_sound > 0.0 {
2847 velocity_magnitude / speed_of_sound
2848 } else {
2849 0.0
2850 };
2851 mach_transitions.record_downward_crossings(
2852 mach_here,
2853 position.x,
2854 &mut transonic_distances,
2855 );
2856 }
2857
2858 if position.y > max_height {
2859 max_height = position.y;
2860 }
2861
2862 if self.inputs.enable_pitch_damping {
2864 let mach = velocity_magnitude / speed_of_sound;
2865
2866 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
2868 transonic_mach = Some(mach);
2869 }
2870
2871 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
2873
2874 if pitch_damping < min_pitch_damping {
2876 min_pitch_damping = pitch_damping;
2877 }
2878 }
2879
2880 if self.inputs.enable_precession_nutation {
2882 if let Some(ref mut state) = angular_state {
2883 let velocity_magnitude = velocity.magnitude();
2884 let params = self.precession_nutation_params(
2885 velocity_magnitude,
2886 air_density,
2887 speed_of_sound,
2888 );
2889
2890 *state = calculate_combined_angular_motion(
2892 ¶ms,
2893 state,
2894 time,
2895 self.time_step,
2896 0.001, );
2898
2899 if state.yaw_angle.abs() > max_yaw_angle {
2901 max_yaw_angle = state.yaw_angle.abs();
2902 }
2903 if state.precession_angle.abs() > max_precession_angle {
2904 max_precession_angle = state.precession_angle.abs();
2905 }
2906 }
2907 }
2908
2909 let dt = self.time_step;
2911
2912 let acc1 = self.calculate_acceleration(
2914 &position,
2915 &velocity,
2916 &wind_vector,
2917 (resolved_temp_c, resolved_press_hpa, base_ratio),
2918 );
2919
2920 let pos2 = position + velocity * (dt * 0.5);
2922 let vel2 = velocity + acc1 * (dt * 0.5);
2923 let acc2 = self.calculate_acceleration(
2924 &pos2,
2925 &vel2,
2926 &wind_vector,
2927 (resolved_temp_c, resolved_press_hpa, base_ratio),
2928 );
2929
2930 let pos3 = position + vel2 * (dt * 0.5);
2932 let vel3 = velocity + acc2 * (dt * 0.5);
2933 let acc3 = self.calculate_acceleration(
2934 &pos3,
2935 &vel3,
2936 &wind_vector,
2937 (resolved_temp_c, resolved_press_hpa, base_ratio),
2938 );
2939
2940 let pos4 = position + vel3 * dt;
2942 let vel4 = velocity + acc3 * dt;
2943 let acc4 = self.calculate_acceleration(
2944 &pos4,
2945 &vel4,
2946 &wind_vector,
2947 (resolved_temp_c, resolved_press_hpa, base_ratio),
2948 );
2949
2950 position += (velocity + vel2 * 2.0 + vel3 * 2.0 + vel4) * (dt / 6.0);
2952 velocity += (acc1 + acc2 * 2.0 + acc3 * 2.0 + acc4) * (dt / 6.0);
2953 time += dt;
2954 self.validate_integration_state(&position, &velocity, time)?;
2955 }
2956
2957 let termination = self.append_terminal_endpoint(
2958 &mut points,
2959 position,
2960 velocity,
2961 time,
2962 &mut max_height,
2963 &mut collapsed_steps,
2964 )?;
2965
2966 self.annotate_drag_coefficients(&mut points, speed_of_sound);
2971
2972 let last_point = points.last().ok_or("No trajectory points generated")?;
2973
2974 let sampled_points =
2976 self.build_sampled_points(&points, max_height, transonic_distances, &mach_transitions)?;
2977
2978 Ok(TrajectoryResult {
2979 max_range: last_point.position.x, max_height,
2981 time_of_flight: last_point.time,
2982 impact_velocity: last_point.velocity_magnitude,
2983 impact_energy: last_point.kinetic_energy,
2984 projectile_mass_kg: self.inputs.bullet_mass,
2985 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
2986 station_speed_of_sound_mps: speed_of_sound,
2987 termination,
2988 points,
2989 sampled_points,
2990 min_pitch_damping: if self.inputs.enable_pitch_damping {
2991 Some(min_pitch_damping)
2992 } else {
2993 None
2994 },
2995 transonic_mach,
2996 angular_state,
2997 max_yaw_angle: if self.inputs.enable_precession_nutation {
2998 Some(max_yaw_angle)
2999 } else {
3000 None
3001 },
3002 max_precession_angle: if self.inputs.enable_precession_nutation {
3003 Some(max_precession_angle)
3004 } else {
3005 None
3006 },
3007 aerodynamic_jump: aj_components,
3008 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
3009 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
3010 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
3011 })
3012 }
3013
3014 fn solve_rk45(&self) -> Result<TrajectoryResult, BallisticsError> {
3015 let mut time = 0.0;
3017 let mut position = self.initial_position();
3021
3022 let aj_components = self.aerodynamic_jump_components();
3028 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
3029 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
3030 let mut velocity = Vector3::new(
3031 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
3035
3036 let mut points = Vec::new();
3037 let mut collapsed_steps = 0;
3040 let mut max_height = position.y;
3041 let mut dt = 0.001; let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
3046 self.resolved_atmosphere();
3047 let base_ratio = air_density / 1.225;
3052 let wind_vector = crate::wind::wind_vector(
3057 self.wind.speed,
3058 self.wind.direction,
3059 self.wind.vertical_speed,
3060 );
3061
3062 let mut transonic_distances: Vec<f64> = Vec::new();
3064 let mut mach_transitions = MachTransitionTracker::default();
3065
3066 let mut min_pitch_damping = f64::INFINITY;
3071 let mut transonic_mach: Option<f64> = None;
3072 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
3073 self.inputs.bullet_model.as_deref().unwrap_or("default"),
3074 );
3075 let mut angular_state = if self.inputs.enable_precession_nutation {
3076 Some(AngularState {
3077 pitch_angle: 0.001,
3078 yaw_angle: 0.001,
3079 pitch_rate: 0.0,
3080 yaw_rate: 0.0,
3081 precession_angle: 0.0,
3082 nutation_phase: 0.0,
3083 })
3084 } else {
3085 None
3086 };
3087 let mut max_yaw_angle = 0.0;
3088 let mut max_precession_angle = 0.0;
3089
3090 while position.x < self.max_range
3091 && position.y > self.inputs.ground_threshold
3092 && time < TRAJECTORY_TIME_LIMIT_S
3093 {
3094 let velocity_magnitude = velocity.magnitude();
3096 let kinetic_energy = 0.5 * self.inputs.bullet_mass * velocity_magnitude.powi(2);
3097
3098 self.push_trajectory_point(
3099 &mut points,
3100 TrajectoryPoint {
3101 time,
3102 position,
3103 velocity_magnitude,
3104 kinetic_energy,
3105 drag_coefficient: None,
3106 },
3107 &mut collapsed_steps,
3108 )?;
3109
3110 {
3113 let mach_here = if speed_of_sound > 0.0 {
3114 velocity_magnitude / speed_of_sound
3115 } else {
3116 0.0
3117 };
3118 mach_transitions.record_downward_crossings(
3119 mach_here,
3120 position.x,
3121 &mut transonic_distances,
3122 );
3123 }
3124
3125 if position.y > max_height {
3126 max_height = position.y;
3127 }
3128
3129 if self.inputs.enable_pitch_damping {
3132 let mach = velocity_magnitude / speed_of_sound;
3133 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
3134 transonic_mach = Some(mach);
3135 }
3136 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
3137 if pitch_damping < min_pitch_damping {
3138 min_pitch_damping = pitch_damping;
3139 }
3140 }
3141
3142 let accepted_step = self.adaptive_rk45_step(
3145 &position,
3146 &velocity,
3147 dt,
3148 &wind_vector,
3149 (resolved_temp_c, resolved_press_hpa, base_ratio),
3150 );
3151 debug_assert!(
3152 accepted_step.error <= RK45_TOLERANCE || accepted_step.used_dt <= RK45_MIN_DT
3153 );
3154
3155 if self.inputs.enable_precession_nutation {
3159 if let Some(ref mut state) = angular_state {
3160 let params = self.precession_nutation_params(
3161 velocity_magnitude,
3162 air_density,
3163 speed_of_sound,
3164 );
3165
3166 *state = calculate_combined_angular_motion(
3167 ¶ms,
3168 state,
3169 time,
3170 accepted_step.used_dt,
3171 0.001,
3172 );
3173
3174 if state.yaw_angle.abs() > max_yaw_angle {
3175 max_yaw_angle = state.yaw_angle.abs();
3176 }
3177 if state.precession_angle.abs() > max_precession_angle {
3178 max_precession_angle = state.precession_angle.abs();
3179 }
3180 }
3181 }
3182
3183 position = accepted_step.position;
3184 velocity = accepted_step.velocity;
3185 time += accepted_step.used_dt;
3186 self.validate_integration_state(&position, &velocity, time)?;
3187
3188 dt = accepted_step.next_dt;
3190 }
3191
3192 if points.is_empty() {
3194 return Err(BallisticsError::from("No trajectory points calculated"));
3195 }
3196
3197 let termination = self.append_terminal_endpoint(
3199 &mut points,
3200 position,
3201 velocity,
3202 time,
3203 &mut max_height,
3204 &mut collapsed_steps,
3205 )?;
3206
3207 self.annotate_drag_coefficients(&mut points, speed_of_sound);
3209
3210 let last_point = points.last().unwrap();
3211
3212 let sampled_points =
3214 self.build_sampled_points(&points, max_height, transonic_distances, &mach_transitions)?;
3215
3216 Ok(TrajectoryResult {
3217 max_range: last_point.position.x, max_height,
3219 time_of_flight: last_point.time,
3220 impact_velocity: last_point.velocity_magnitude,
3221 impact_energy: last_point.kinetic_energy,
3222 projectile_mass_kg: self.inputs.bullet_mass,
3223 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
3224 station_speed_of_sound_mps: speed_of_sound,
3225 termination,
3226 points,
3227 sampled_points,
3228 min_pitch_damping: if self.inputs.enable_pitch_damping {
3229 Some(min_pitch_damping)
3230 } else {
3231 None
3232 },
3233 transonic_mach,
3234 angular_state,
3235 max_yaw_angle: if self.inputs.enable_precession_nutation {
3236 Some(max_yaw_angle)
3237 } else {
3238 None
3239 },
3240 max_precession_angle: if self.inputs.enable_precession_nutation {
3241 Some(max_precession_angle)
3242 } else {
3243 None
3244 },
3245 aerodynamic_jump: aj_components,
3246 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
3247 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
3248 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
3249 })
3250 }
3251
3252 fn adaptive_rk45_step(
3253 &self,
3254 position: &Vector3<f64>,
3255 velocity: &Vector3<f64>,
3256 initial_dt: f64,
3257 wind_vector: &Vector3<f64>,
3258 resolved_atmo: (f64, f64, f64),
3259 ) -> Rk45AcceptedStep {
3260 let mut trial_dt = initial_dt;
3261
3262 loop {
3263 let trial = self.rk45_step(
3264 position,
3265 velocity,
3266 trial_dt,
3267 wind_vector,
3268 RK45_TOLERANCE,
3269 resolved_atmo,
3270 );
3271 let next_dt = if trial.suggested_dt.is_finite() {
3276 (RK45_SAFETY_FACTOR * trial.suggested_dt).clamp(RK45_MIN_DT, RK45_MAX_DT)
3277 } else {
3278 RK45_MIN_DT
3279 };
3280
3281 if trial.error <= RK45_TOLERANCE || trial_dt <= RK45_MIN_DT {
3282 return Rk45AcceptedStep {
3283 position: trial.position,
3284 velocity: trial.velocity,
3285 used_dt: trial_dt,
3286 next_dt,
3287 error: trial.error,
3288 };
3289 }
3290
3291 trial_dt = next_dt;
3292 }
3293 }
3294
3295 fn rk45_step(
3296 &self,
3297 position: &Vector3<f64>,
3298 velocity: &Vector3<f64>,
3299 dt: f64,
3300 wind_vector: &Vector3<f64>,
3301 tolerance: f64,
3302 resolved_atmo: (f64, f64, f64), ) -> Rk45Trial {
3304 const A21: f64 = 1.0 / 5.0;
3306 const A31: f64 = 3.0 / 40.0;
3307 const A32: f64 = 9.0 / 40.0;
3308 const A41: f64 = 44.0 / 45.0;
3309 const A42: f64 = -56.0 / 15.0;
3310 const A43: f64 = 32.0 / 9.0;
3311 const A51: f64 = 19372.0 / 6561.0;
3312 const A52: f64 = -25360.0 / 2187.0;
3313 const A53: f64 = 64448.0 / 6561.0;
3314 const A54: f64 = -212.0 / 729.0;
3315 const A61: f64 = 9017.0 / 3168.0;
3316 const A62: f64 = -355.0 / 33.0;
3317 const A63: f64 = 46732.0 / 5247.0;
3318 const A64: f64 = 49.0 / 176.0;
3319 const A65: f64 = -5103.0 / 18656.0;
3320 const A71: f64 = 35.0 / 384.0;
3321 const A73: f64 = 500.0 / 1113.0;
3322 const A74: f64 = 125.0 / 192.0;
3323 const A75: f64 = -2187.0 / 6784.0;
3324 const A76: f64 = 11.0 / 84.0;
3325
3326 const B1: f64 = 35.0 / 384.0;
3328 const B3: f64 = 500.0 / 1113.0;
3329 const B4: f64 = 125.0 / 192.0;
3330 const B5: f64 = -2187.0 / 6784.0;
3331 const B6: f64 = 11.0 / 84.0;
3332
3333 const B1_ERR: f64 = 5179.0 / 57600.0;
3335 const B3_ERR: f64 = 7571.0 / 16695.0;
3336 const B4_ERR: f64 = 393.0 / 640.0;
3337 const B5_ERR: f64 = -92097.0 / 339200.0;
3338 const B6_ERR: f64 = 187.0 / 2100.0;
3339 const B7_ERR: f64 = 1.0 / 40.0;
3340
3341 let k1_v = self.calculate_acceleration(position, velocity, wind_vector, resolved_atmo);
3343 let k1_p = *velocity;
3344
3345 let p2 = position + dt * A21 * k1_p;
3346 let v2 = velocity + dt * A21 * k1_v;
3347 let k2_v = self.calculate_acceleration(&p2, &v2, wind_vector, resolved_atmo);
3348 let k2_p = v2;
3349
3350 let p3 = position + dt * (A31 * k1_p + A32 * k2_p);
3351 let v3 = velocity + dt * (A31 * k1_v + A32 * k2_v);
3352 let k3_v = self.calculate_acceleration(&p3, &v3, wind_vector, resolved_atmo);
3353 let k3_p = v3;
3354
3355 let p4 = position + dt * (A41 * k1_p + A42 * k2_p + A43 * k3_p);
3356 let v4 = velocity + dt * (A41 * k1_v + A42 * k2_v + A43 * k3_v);
3357 let k4_v = self.calculate_acceleration(&p4, &v4, wind_vector, resolved_atmo);
3358 let k4_p = v4;
3359
3360 let p5 = position + dt * (A51 * k1_p + A52 * k2_p + A53 * k3_p + A54 * k4_p);
3361 let v5 = velocity + dt * (A51 * k1_v + A52 * k2_v + A53 * k3_v + A54 * k4_v);
3362 let k5_v = self.calculate_acceleration(&p5, &v5, wind_vector, resolved_atmo);
3363 let k5_p = v5;
3364
3365 let p6 = position + dt * (A61 * k1_p + A62 * k2_p + A63 * k3_p + A64 * k4_p + A65 * k5_p);
3366 let v6 = velocity + dt * (A61 * k1_v + A62 * k2_v + A63 * k3_v + A64 * k4_v + A65 * k5_v);
3367 let k6_v = self.calculate_acceleration(&p6, &v6, wind_vector, resolved_atmo);
3368 let k6_p = v6;
3369
3370 let p7 = position + dt * (A71 * k1_p + A73 * k3_p + A74 * k4_p + A75 * k5_p + A76 * k6_p);
3371 let v7 = velocity + dt * (A71 * k1_v + A73 * k3_v + A74 * k4_v + A75 * k5_v + A76 * k6_v);
3372 let k7_v = self.calculate_acceleration(&p7, &v7, wind_vector, resolved_atmo);
3373 let k7_p = v7;
3374
3375 let new_pos = position + dt * (B1 * k1_p + B3 * k3_p + B4 * k4_p + B5 * k5_p + B6 * k6_p);
3377 let new_vel = velocity + dt * (B1 * k1_v + B3 * k3_v + B4 * k4_v + B5 * k5_v + B6 * k6_v);
3378
3379 let pos_err = position
3381 + dt * (B1_ERR * k1_p
3382 + B3_ERR * k3_p
3383 + B4_ERR * k4_p
3384 + B5_ERR * k5_p
3385 + B6_ERR * k6_p
3386 + B7_ERR * k7_p);
3387 let vel_err = velocity
3388 + dt * (B1_ERR * k1_v
3389 + B3_ERR * k3_v
3390 + B4_ERR * k4_v
3391 + B5_ERR * k5_v
3392 + B6_ERR * k6_v
3393 + B7_ERR * k7_v);
3394
3395 let error = cli_rk45_error_norm(position, velocity, &new_pos, &new_vel, &pos_err, &vel_err);
3397
3398 let dt_new = if error < tolerance {
3400 dt * (tolerance / error).powf(0.2).min(2.0)
3401 } else {
3402 dt * (tolerance / error).powf(0.25).max(0.1)
3403 };
3404
3405 Rk45Trial {
3406 position: new_pos,
3407 velocity: new_vel,
3408 suggested_dt: dt_new,
3409 error,
3410 }
3411 }
3412
3413 fn apply_cluster_bc_correction(&self, base_bc: f64, velocity_fps: f64) -> f64 {
3414 if let Some(ref cluster_bc) = self.cluster_bc {
3415 cluster_bc.apply_correction_for_drag_model(
3416 base_bc,
3417 self.inputs.caliber_inches,
3418 self.inputs.weight_grains,
3419 velocity_fps,
3420 self.inputs.bc_type,
3421 )
3422 } else {
3423 base_bc
3424 }
3425 }
3426
3427 fn calculate_acceleration(
3428 &self,
3429 position: &Vector3<f64>,
3430 velocity: &Vector3<f64>,
3431 wind_vector: &Vector3<f64>,
3432 resolved_atmo: (f64, f64, f64), ) -> Vector3<f64> {
3434 let actual_wind = if let Some(ref sock) = self.wind_sock {
3440 sock.vector_for_range_stateless(position.x)
3441 } else if self.inputs.enable_wind_shear {
3442 self.get_wind_at_altitude(position.y)
3443 } else {
3444 *wind_vector
3445 };
3446 let actual_wind =
3447 crate::derivatives::level_vector_to_shot_frame(actual_wind, self.inputs.shooting_angle);
3448
3449 let relative_velocity = velocity - actual_wind;
3450 let velocity_magnitude = relative_velocity.magnitude();
3451
3452 if velocity_magnitude < 0.001 {
3453 return self.gravity_acceleration();
3454 }
3455
3456 let (base_temp_c, base_press_hpa, station_ratio) = resolved_atmo;
3467
3468 let (drag_base_temp_c, drag_base_press_hpa, drag_base_ratio, drag_humidity_percent) =
3476 if let Some(ref sock) = self.atmo_sock {
3477 let (zone_temp_c, zone_press_hpa, zone_humidity) = sock.atmo_for_range(position.x);
3478 let zone_base_ratio = crate::atmosphere::calculate_air_density_cimp(
3479 zone_temp_c,
3480 zone_press_hpa,
3481 zone_humidity,
3482 ) / 1.225;
3483 (zone_temp_c, zone_press_hpa, zone_base_ratio, zone_humidity)
3484 } else {
3485 (
3486 base_temp_c,
3487 base_press_hpa,
3488 station_ratio,
3489 self.atmosphere.humidity,
3490 )
3491 };
3492 let local_alt = crate::atmosphere::shot_frame_altitude(
3493 self.atmosphere.altitude,
3494 position.x,
3495 position.y,
3496 self.inputs.shooting_angle,
3497 );
3498 let (air_density, speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
3499 local_alt,
3500 self.atmosphere.altitude,
3501 drag_base_temp_c,
3502 drag_base_press_hpa,
3503 drag_base_ratio,
3504 drag_humidity_percent,
3505 );
3506
3507 let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
3511
3512 let velocity_fps = velocity_magnitude * 3.28084;
3514
3515 let cd_to_retard = crate::constants::CD_TO_RETARD;
3520 let standard_factor = cd * cd_to_retard;
3521 let density_scale = air_density / 1.225; let a_drag_ft_s2 =
3525 (velocity_fps * velocity_fps) * standard_factor * density_scale / retard_denom;
3526 let a_drag_m_s2 = a_drag_ft_s2 * 0.3048; let drag_acceleration = -a_drag_m_s2 * (relative_velocity / velocity_magnitude);
3530
3531 let mut accel = drag_acceleration + self.gravity_acceleration();
3534
3535 if self.inputs.enable_coriolis {
3538 if let Some(lat_deg) = self.inputs.latitude {
3539 let omega_earth = 7.2921159e-5_f64; let lat = lat_deg.to_radians();
3541 let az = self.inputs.shot_azimuth; let omega = Vector3::new(
3548 omega_earth * lat.cos() * az.cos(), omega_earth * lat.sin(), -omega_earth * lat.cos() * az.sin(), );
3552 let omega = crate::derivatives::level_vector_to_shot_frame(
3553 omega,
3554 self.inputs.shooting_angle,
3555 );
3556 accel += -2.0 * omega.cross(velocity);
3561 }
3562 }
3563
3564 if self.inputs.enable_magnus
3571 && !self.inputs.use_enhanced_spin_drift
3572 && self.inputs.bullet_diameter > 0.0
3573 && self.inputs.twist_rate > 0.0
3574 {
3575 let diameter_m = self.inputs.bullet_diameter;
3576 let (spin_rad_s, spin_param) = crate::spin_drift::calculate_magnus_spin_state(
3577 self.inputs.muzzle_velocity,
3578 velocity_magnitude,
3579 self.inputs.twist_rate,
3580 diameter_m,
3581 );
3582 let mach = velocity_magnitude / speed_of_sound;
3584
3585 let d_in = self.inputs.bullet_diameter / 0.0254;
3587 let m_gr = self.inputs.bullet_mass / crate::constants::GRAINS_TO_KG;
3588 let l_in = if self.inputs.bullet_length > 0.0 {
3589 self.inputs.bullet_length / 0.0254
3590 } else {
3591 let est_m = crate::stability::estimate_bullet_length_m(
3593 self.inputs.bullet_diameter,
3594 self.inputs.bullet_mass,
3595 );
3596 if est_m > 0.0 {
3597 est_m / 0.0254
3598 } else {
3599 4.5 * d_in
3600 }
3601 };
3602 let sg = crate::spin_drift::calculate_dynamic_stability(
3606 m_gr,
3607 velocity_magnitude,
3608 spin_rad_s,
3609 d_in,
3610 l_in,
3611 air_density,
3612 );
3613
3614 let (yaw_rad, _) = crate::spin_drift::calculate_yaw_of_repose(
3616 sg,
3617 velocity_magnitude,
3618 spin_rad_s,
3619 0.0, 0.0, air_density,
3622 d_in,
3623 l_in,
3624 m_gr,
3625 mach,
3626 "match",
3627 false,
3628 );
3629
3630 let c_np = crate::derivatives::calculate_magnus_moment_coefficient(mach);
3632 let area = std::f64::consts::PI * (diameter_m / 2.0).powi(2);
3633 let magnus_force = 0.5
3634 * air_density
3635 * velocity_magnitude.powi(2)
3636 * area
3637 * c_np
3638 * spin_param
3639 * yaw_rad.sin();
3640
3641 if magnus_force.abs() > 1e-12 {
3645 if let Some(dir) = crate::derivatives::yaw_of_repose_magnus_direction(
3646 relative_velocity,
3647 self.gravity_acceleration(),
3648 self.inputs.is_twist_right,
3649 ) {
3650 accel += (magnus_force / self.inputs.bullet_mass) * dir;
3651 }
3652 }
3653 }
3654
3655 accel
3656 }
3657
3658 fn drag_terms(&self, velocity_magnitude: f64, speed_of_sound: f64) -> (f64, f64) {
3670 let cd = self.calculate_drag_coefficient(velocity_magnitude, speed_of_sound);
3671
3672 let velocity_fps = velocity_magnitude * 3.28084;
3673
3674 let (base_bc, bc_from_segments) = if let Some(segments) = self
3679 .inputs
3680 .bc_segments_data
3681 .as_ref()
3682 .filter(|segments| self.inputs.use_bc_segments && !segments.is_empty())
3683 {
3684 (
3686 crate::bc_estimation::velocity_segment_bc(
3687 velocity_fps,
3688 segments,
3689 self.inputs.bc_value,
3690 ),
3691 true,
3692 )
3693 } else if let Some(segments) = self
3694 .inputs
3695 .bc_segments
3696 .as_ref()
3697 .filter(|segments| !segments.is_empty())
3698 {
3699 (
3700 crate::derivatives::interpolated_bc(
3701 velocity_magnitude / speed_of_sound,
3702 segments,
3703 Some(&self.inputs),
3704 ),
3705 true,
3706 )
3707 } else {
3708 (self.inputs.bc_value, false)
3709 };
3710
3711 let effective_bc = if bc_from_segments {
3716 base_bc
3717 } else {
3718 self.apply_cluster_bc_correction(base_bc, velocity_fps)
3719 };
3720 let effective_bc = effective_bc.max(1e-6);
3723
3724 let retard_denom = if self.inputs.custom_drag_table.is_some() {
3729 self.inputs.custom_drag_denominator(effective_bc)
3730 } else {
3731 effective_bc
3732 };
3733
3734 (cd, retard_denom)
3735 }
3736
3737 pub fn effective_drag_coefficient(
3758 &self,
3759 velocity_magnitude: f64,
3760 speed_of_sound: f64,
3761 ) -> Option<f64> {
3762 if !velocity_magnitude.is_finite() || speed_of_sound <= 1e-9 {
3763 return None;
3764 }
3765 let sectional_density = self.inputs.sectional_density_lb_in2()?;
3766 let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
3767 if retard_denom <= 0.0 {
3768 return None;
3769 }
3770 let effective = cd * sectional_density / retard_denom;
3771 effective.is_finite().then_some(effective)
3772 }
3773
3774 fn annotate_drag_coefficients(&self, points: &mut [TrajectoryPoint], speed_of_sound: f64) {
3784 for point in points.iter_mut() {
3785 point.drag_coefficient =
3786 self.effective_drag_coefficient(point.velocity_magnitude, speed_of_sound);
3787 }
3788 }
3789
3790 fn calculate_drag_coefficient(&self, velocity: f64, speed_of_sound: f64) -> f64 {
3791 let mach = velocity / speed_of_sound;
3792
3793 if let Some(ref table) = self.inputs.custom_drag_table {
3797 return table.interpolate(mach) * self.inputs.cd_scale;
3802 }
3803
3804 crate::drag::get_drag_coefficient(mach, &self.inputs.bc_type)
3807 }
3808}
3809
3810#[derive(Debug, Clone)]
3812pub struct MonteCarloParams {
3813 pub num_simulations: usize,
3814 pub velocity_std_dev: f64,
3815 pub angle_std_dev: f64,
3816 pub bc_std_dev: f64,
3817 pub wind_speed_std_dev: f64,
3818 pub target_distance: Option<f64>,
3819 pub base_wind_speed: f64,
3820 pub base_wind_direction: f64,
3821 pub azimuth_std_dev: f64, }
3823
3824impl Default for MonteCarloParams {
3825 fn default() -> Self {
3826 Self {
3827 num_simulations: 1000,
3828 velocity_std_dev: 1.0,
3829 angle_std_dev: 0.001,
3830 bc_std_dev: 0.01,
3831 wind_speed_std_dev: 1.0,
3832 target_distance: None,
3833 base_wind_speed: 0.0,
3834 base_wind_direction: 0.0,
3835 azimuth_std_dev: 0.001, }
3837 }
3838}
3839
3840#[derive(Debug, Clone)]
3842pub struct MonteCarloResults {
3843 pub ranges: Vec<f64>,
3844 pub impact_velocities: Vec<f64>,
3845 pub impact_positions: Vec<Vector3<f64>>,
3851}
3852
3853pub const DEFAULT_HIT_RADIUS_M: f64 = 0.3;
3856
3857pub const TARGET_NOT_REACHED_SENTINEL_M: f64 = -1.0e9;
3863
3864impl MonteCarloResults {
3865 pub fn position_reached_target(position: &Vector3<f64>) -> bool {
3867 position.iter().all(|component| component.is_finite())
3868 && position.y != TARGET_NOT_REACHED_SENTINEL_M
3869 }
3870
3871 pub fn target_arrival_count(&self) -> usize {
3873 self.impact_positions
3874 .iter()
3875 .filter(|position| Self::position_reached_target(position))
3876 .count()
3877 }
3878
3879 pub fn target_shortfall_fraction(&self) -> f64 {
3882 if self.impact_positions.is_empty() {
3883 return 0.0;
3884 }
3885 (self.impact_positions.len() - self.target_arrival_count()) as f64
3886 / self.impact_positions.len() as f64
3887 }
3888
3889 pub fn target_plane_cep(&self) -> Option<f64> {
3895 let mut radial_misses: Vec<f64> = self
3896 .impact_positions
3897 .iter()
3898 .filter(|position| Self::position_reached_target(position))
3899 .map(Vector3::norm)
3900 .filter(|miss| miss.is_finite())
3901 .collect();
3902 radial_misses.sort_by(f64::total_cmp);
3903 if radial_misses.is_empty() {
3904 None
3905 } else {
3906 Some(radial_misses[radial_misses.len() / 2])
3907 }
3908 }
3909
3910 pub fn hit_probability(&self, hit_radius_m: f64) -> f64 {
3919 if self.impact_positions.is_empty() {
3920 return 0.0;
3921 }
3922 let hits = self
3923 .impact_positions
3924 .iter()
3925 .filter(|position| Self::position_is_hit(position, hit_radius_m))
3926 .count();
3927 hits as f64 / self.impact_positions.len() as f64
3928 }
3929
3930 pub fn position_is_hit(position: &Vector3<f64>, hit_radius_m: f64) -> bool {
3942 Self::position_reached_target(position) && position.norm() < hit_radius_m
3943 }
3944
3945 pub fn hit_probability_wilson(
3967 &self,
3968 hit_radius_m: f64,
3969 level: ConfidenceLevel,
3970 ) -> (f64, (f64, f64), u64) {
3971 let trials = self.impact_positions.len() as u64;
3972 let hits = self
3973 .impact_positions
3974 .iter()
3975 .filter(|position| Self::position_is_hit(position, hit_radius_m))
3976 .count() as u64;
3977 (
3978 self.hit_probability(hit_radius_m),
3979 wilson_interval(hits, trials, level),
3980 trials,
3981 )
3982 }
3983
3984 pub fn rect_hit_probability(&self, width_m: f64, height_m: f64) -> f64 {
3996 let dimensions_invalid =
3997 width_m.is_nan() || width_m <= 0.0 || height_m.is_nan() || height_m <= 0.0;
3998 if self.impact_positions.is_empty() || dimensions_invalid {
3999 return 0.0;
4000 }
4001 let half_width = width_m / 2.0;
4002 let half_height = height_m / 2.0;
4003 let hits = self
4004 .impact_positions
4005 .iter()
4006 .filter(|position| {
4007 Self::position_reached_target(position)
4008 && position.z.abs() <= half_width
4009 && position.y.abs() <= half_height
4010 })
4011 .count();
4012 hits as f64 / self.impact_positions.len() as f64
4013 }
4014}
4015
4016fn wind_from_signed_speed_sample(
4017 signed_speed: f64,
4018 sampled_direction: f64,
4019 vertical_speed: f64,
4020) -> WindConditions {
4021 if signed_speed < 0.0 {
4026 WindConditions {
4027 speed: -signed_speed,
4028 direction: sampled_direction + std::f64::consts::PI,
4029 vertical_speed,
4030 }
4031 } else {
4032 WindConditions {
4033 speed: signed_speed,
4034 direction: sampled_direction,
4035 vertical_speed,
4036 }
4037 }
4038}
4039
4040struct MonteCarloWindSampler {
4041 speed: rand_distr::Normal<f64>,
4042 direction: rand_distr::Normal<f64>,
4043 vertical_speed: f64,
4045}
4046
4047impl MonteCarloWindSampler {
4048 fn new(
4049 base_wind: &WindConditions,
4050 wind_speed_std_dev: f64,
4051 wind_direction_std_dev: f64,
4052 ) -> Result<Self, BallisticsError> {
4053 use rand_distr::Normal;
4054
4055 if !wind_direction_std_dev.is_finite() || wind_direction_std_dev < 0.0 {
4056 return Err("Wind direction standard deviation must be finite and non-negative".into());
4057 }
4058
4059 let speed = Normal::new(base_wind.speed, wind_speed_std_dev)
4060 .map_err(|e| format!("Invalid wind speed distribution: {e}"))?;
4061 let direction = Normal::new(base_wind.direction, wind_direction_std_dev)
4062 .map_err(|e| format!("Invalid wind direction distribution: {e}"))?;
4063 Ok(Self {
4064 speed,
4065 direction,
4066 vertical_speed: base_wind.vertical_speed,
4067 })
4068 }
4069
4070 fn sample<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> WindConditions {
4071 use rand_distr::Distribution;
4072
4073 wind_from_signed_speed_sample(
4074 self.speed.sample(rng),
4075 self.direction.sample(rng),
4076 self.vertical_speed,
4077 )
4078 }
4079}
4080
4081#[derive(Debug, Clone, Copy)]
4084struct TrialOutcome {
4085 range: f64,
4088 impact_velocity: f64,
4090 impact_position: Vector3<f64>,
4093}
4094
4095struct MonteCarloTrialSampler {
4115 base_inputs: BallisticInputs,
4116 atmosphere: AtmosphericConditions,
4117 solver_max_range: f64,
4118 target_distance: f64,
4121 baseline_at_target: Vector3<f64>,
4124 velocity_delta_dist: rand_distr::Normal<f64>,
4125 angle_dist: rand_distr::Normal<f64>,
4126 bc_dist: rand_distr::Normal<f64>,
4127 wind_sampler: MonteCarloWindSampler,
4128 azimuth_dist: rand_distr::Normal<f64>,
4129}
4130
4131impl MonteCarloTrialSampler {
4132 fn new(
4136 base_inputs: BallisticInputs,
4137 base_wind: &WindConditions,
4138 params: &MonteCarloParams,
4139 wind_direction_std_dev: f64,
4140 ) -> Result<Self, BallisticsError> {
4141 use rand_distr::Normal;
4142
4143 let atmosphere = AtmosphericConditions {
4144 temperature: base_inputs.temperature,
4145 pressure: base_inputs.pressure,
4146 humidity: base_inputs.humidity_percent(),
4147 altitude: base_inputs.altitude,
4148 };
4149 let target_hint = params
4150 .target_distance
4151 .unwrap_or(base_inputs.target_distance);
4152 let solver_max_range = target_hint.max(1000.0) * 2.0;
4153
4154 let mut baseline_solver =
4156 TrajectorySolver::new(base_inputs.clone(), base_wind.clone(), atmosphere.clone());
4157 baseline_solver.set_max_range(solver_max_range);
4158 let baseline_result = baseline_solver.solve()?;
4159
4160 let target_distance = params.target_distance.unwrap_or(baseline_result.max_range);
4162
4163 let baseline_at_target = baseline_result
4165 .position_at_range(target_distance)
4166 .ok_or("Could not interpolate baseline at target distance")?;
4167
4168 let velocity_delta_dist = Normal::new(0.0, params.velocity_std_dev)
4173 .map_err(|e| format!("Invalid velocity distribution: {}", e))?;
4174 let angle_dist = Normal::new(base_inputs.muzzle_angle, params.angle_std_dev)
4175 .map_err(|e| format!("Invalid angle distribution: {}", e))?;
4176 let bc_dist = Normal::new(base_inputs.bc_value, params.bc_std_dev)
4177 .map_err(|e| format!("Invalid BC distribution: {}", e))?;
4178 let wind_sampler = MonteCarloWindSampler::new(
4181 base_wind,
4182 params.wind_speed_std_dev,
4183 wind_direction_std_dev,
4184 )?;
4185 let azimuth_dist = Normal::new(base_inputs.azimuth_angle, params.azimuth_std_dev)
4186 .map_err(|e| format!("Invalid azimuth distribution: {}", e))?;
4187
4188 Ok(Self {
4189 base_inputs,
4190 atmosphere,
4191 solver_max_range,
4192 target_distance,
4193 baseline_at_target,
4194 velocity_delta_dist,
4195 angle_dist,
4196 bc_dist,
4197 wind_sampler,
4198 azimuth_dist,
4199 })
4200 }
4201
4202 fn sample_one_trial<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> Option<TrialOutcome> {
4212 use rand_distr::Distribution;
4213
4214 let mut inputs = self.base_inputs.clone();
4216 let muzzle_velocity_delta = self.velocity_delta_dist.sample(&mut *rng);
4217 inputs.muzzle_angle = self.angle_dist.sample(&mut *rng);
4218 inputs.bc_value = self.bc_dist.sample(&mut *rng).max(0.01);
4219 inputs.azimuth_angle = self.azimuth_dist.sample(&mut *rng); let wind = self.wind_sampler.sample(&mut *rng);
4223
4224 let mut solver = TrajectorySolver::new(inputs, wind, self.atmosphere.clone());
4228 solver.inputs.muzzle_velocity =
4229 (solver.inputs.muzzle_velocity + muzzle_velocity_delta).max(0.0);
4230 solver.set_max_range(self.solver_max_range);
4231 let result = solver.solve().ok()?;
4233
4234 let impact_position = if result.max_range < self.target_distance {
4240 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)
4243 } else {
4244 let pos_at_target = result.position_at_range(self.target_distance)?;
4246 Vector3::new(
4251 0.0,
4252 pos_at_target.y - self.baseline_at_target.y,
4253 pos_at_target.z - self.baseline_at_target.z,
4254 )
4255 };
4256
4257 Some(TrialOutcome {
4258 range: result.max_range,
4259 impact_velocity: result.impact_velocity,
4260 impact_position,
4261 })
4262 }
4263}
4264
4265pub fn run_monte_carlo(
4267 base_inputs: BallisticInputs,
4268 params: MonteCarloParams,
4269) -> Result<MonteCarloResults, BallisticsError> {
4270 run_monte_carlo_with_direction_std_dev(base_inputs, params, 0.0)
4271}
4272
4273pub fn run_monte_carlo_with_direction_std_dev(
4278 base_inputs: BallisticInputs,
4279 params: MonteCarloParams,
4280 wind_direction_std_dev: f64,
4281) -> Result<MonteCarloResults, BallisticsError> {
4282 let base_wind = WindConditions {
4283 speed: params.base_wind_speed,
4284 direction: params.base_wind_direction,
4285 vertical_speed: 0.0,
4286 };
4287 run_monte_carlo_with_wind_and_direction_std_dev(
4288 base_inputs,
4289 base_wind,
4290 params,
4291 wind_direction_std_dev,
4292 )
4293}
4294
4295pub fn run_monte_carlo_with_wind(
4297 base_inputs: BallisticInputs,
4298 base_wind: WindConditions,
4299 params: MonteCarloParams,
4300) -> Result<MonteCarloResults, BallisticsError> {
4301 run_monte_carlo_with_wind_and_direction_std_dev(base_inputs, base_wind, params, 0.0)
4302}
4303
4304pub fn run_monte_carlo_with_wind_and_direction_std_dev(
4309 base_inputs: BallisticInputs,
4310 base_wind: WindConditions,
4311 params: MonteCarloParams,
4312 wind_direction_std_dev: f64,
4313) -> Result<MonteCarloResults, BallisticsError> {
4314 let mut rng = rand::rng();
4315 run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
4316 base_inputs,
4317 base_wind,
4318 params,
4319 wind_direction_std_dev,
4320 &mut rng,
4321 )
4322}
4323
4324pub fn run_monte_carlo_with_wind_and_direction_std_dev_seeded(
4331 base_inputs: BallisticInputs,
4332 base_wind: WindConditions,
4333 params: MonteCarloParams,
4334 wind_direction_std_dev: f64,
4335 seed: u64,
4336) -> Result<MonteCarloResults, BallisticsError> {
4337 use rand::{rngs::StdRng, SeedableRng};
4338 let mut rng = StdRng::seed_from_u64(seed);
4339 run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
4340 base_inputs,
4341 base_wind,
4342 params,
4343 wind_direction_std_dev,
4344 &mut rng,
4345 )
4346}
4347
4348fn run_monte_carlo_with_wind_and_direction_std_dev_using_rng<R: rand::Rng + ?Sized>(
4349 base_inputs: BallisticInputs,
4350 base_wind: WindConditions,
4351 params: MonteCarloParams,
4352 wind_direction_std_dev: f64,
4353 rng: &mut R,
4354) -> Result<MonteCarloResults, BallisticsError> {
4355 let mut ranges = Vec::new();
4356 let mut impact_velocities = Vec::new();
4357 let mut impact_positions = Vec::new();
4358
4359 let sampler =
4360 MonteCarloTrialSampler::new(base_inputs, &base_wind, ¶ms, wind_direction_std_dev)?;
4361
4362 for _ in 0..params.num_simulations {
4363 if let Some(outcome) = sampler.sample_one_trial(rng) {
4367 ranges.push(outcome.range);
4368 impact_velocities.push(outcome.impact_velocity);
4369 impact_positions.push(outcome.impact_position);
4370 }
4371 }
4372
4373 if ranges.is_empty() {
4374 return Err("No successful simulations".into());
4375 }
4376
4377 Ok(MonteCarloResults {
4378 ranges,
4379 impact_velocities,
4380 impact_positions,
4381 })
4382}
4383
4384pub const MC_ADAPTIVE_SCHEMA_VERSION_V1: u32 = 1;
4386
4387pub const MC_ADAPTIVE_METHOD_V1: &str = "anytime_beta_binomial_mixture_cs_v1";
4394
4395pub const MC_ADAPTIVE_ASSUMPTIONS_V1: [&str; 4] = [
4409 "Sampling uncertainty only: intervals cover Monte Carlo sampling error, not model error in the trajectory solver or its inputs.",
4410 "Anytime-valid stopping: the beta-binomial mixture confidence sequence keeps its coverage guarantee despite stopping the moment the target half-width is met.",
4411 "Input dispersions are the independent normal distributions declared in MonteCarloParams; correlations between inputs are not modeled.",
4412 "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.",
4413];
4414
4415#[derive(Debug, Clone)]
4422pub struct McConvergence {
4423 pub level: ConfidenceLevel,
4425 pub target_half_width: f64,
4428 pub min_samples: u64,
4432 pub max_samples: u64,
4435 pub batch_size: u64,
4438}
4439
4440impl Default for McConvergence {
4441 fn default() -> Self {
4442 Self {
4443 level: ConfidenceLevel::P95,
4444 target_half_width: 0.02,
4445 min_samples: 1_000,
4446 max_samples: 100_000,
4447 batch_size: 500,
4448 }
4449 }
4450}
4451
4452impl McConvergence {
4453 pub fn validate(&self) -> Result<(), String> {
4462 if !self.target_half_width.is_finite() || self.target_half_width <= 0.0 {
4463 return Err(format!(
4464 "McConvergence.target_half_width must be a finite value greater than zero (got {})",
4465 self.target_half_width
4466 ));
4467 }
4468 if self.batch_size == 0 {
4469 return Err("McConvergence.batch_size must be greater than zero".to_string());
4470 }
4471 if self.max_samples == 0 {
4472 return Err("McConvergence.max_samples must be greater than zero".to_string());
4473 }
4474 if self.max_samples < self.min_samples {
4475 return Err(format!(
4476 "McConvergence.max_samples ({}) must be at least McConvergence.min_samples ({})",
4477 self.max_samples, self.min_samples
4478 ));
4479 }
4480 Ok(())
4481 }
4482}
4483
4484#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize)]
4486#[serde(rename_all = "snake_case")]
4487pub enum McStopReason {
4488 TargetHalfWidthMet,
4491 MaxSamplesReached,
4494}
4495
4496#[derive(Debug, Clone, serde::Serialize)]
4503pub struct AdaptiveMcReportV1 {
4504 pub schema_version: u32,
4506 pub method: String,
4508 pub assumptions: Vec<String>,
4510 pub confidence_percent: u32,
4512 pub hit_probability: f64,
4515 pub ci_low: f64,
4517 pub ci_high: f64,
4519 pub samples: u64,
4521 pub attempts: u64,
4532 pub arrivals: u64,
4547 pub stop_reason: McStopReason,
4549 pub hit_radius_m: f64,
4551 pub target_distance_m: f64,
4553 pub mean_impact_velocity_mps: f64,
4558 pub std_impact_velocity_mps: f64,
4560 pub mean_drop_at_target_m: f64,
4568 pub std_drop_at_target_m: f64,
4571 pub mean_wind_drift_at_target_m: f64,
4574 pub std_wind_drift_at_target_m: f64,
4577}
4578
4579pub fn run_monte_carlo_adaptive_seeded(
4654 base_inputs: &BallisticInputs,
4655 base_wind: &WindConditions,
4656 params: &MonteCarloParams,
4657 convergence: &McConvergence,
4658 hit_radius_m: f64,
4659 seed: u64,
4660) -> Result<AdaptiveMcReportV1, String> {
4661 use rand::{rngs::StdRng, SeedableRng};
4662
4663 convergence.validate()?;
4664
4665 let sampler = MonteCarloTrialSampler::new(base_inputs.clone(), base_wind, params, 0.0)
4668 .map_err(|e| e.to_string())?;
4669
4670 let mut rng = StdRng::seed_from_u64(seed);
4671 let mut hits_cs = BernoulliConfidenceSequence::new(convergence.level);
4672 let mut impact_velocity = Welford::new();
4673 let mut drop_at_target = Welford::new();
4674 let mut drift_at_target = Welford::new();
4675
4676 let mut attempts: u64 = 0;
4677 let mut stop_reason = McStopReason::MaxSamplesReached;
4678
4679 while attempts < convergence.max_samples {
4680 let batch = convergence
4683 .batch_size
4684 .min(convergence.max_samples - attempts);
4685 let mut batch_hits: u64 = 0;
4686 let mut batch_trials: u64 = 0;
4687
4688 for _ in 0..batch {
4689 attempts += 1;
4690 let Some(outcome) = sampler.sample_one_trial(&mut rng) else {
4691 continue; };
4693 batch_trials += 1;
4694 if MonteCarloResults::position_is_hit(&outcome.impact_position, hit_radius_m) {
4695 batch_hits += 1;
4696 }
4697 if MonteCarloResults::position_reached_target(&outcome.impact_position) {
4701 impact_velocity.push(outcome.impact_velocity);
4702 drop_at_target.push(outcome.impact_position.y);
4703 drift_at_target.push(outcome.impact_position.z);
4704 }
4705 }
4706
4707 hits_cs.update_batch(batch_hits, batch_trials);
4708
4709 if hits_cs.trials() >= convergence.min_samples
4710 && hits_cs.half_width() <= convergence.target_half_width
4711 {
4712 stop_reason = McStopReason::TargetHalfWidthMet;
4713 break;
4714 }
4715 }
4716
4717 let samples = hits_cs.trials();
4718 if samples == 0 {
4719 return Err("No successful simulations".to_string());
4720 }
4721 let (ci_low, ci_high) = hits_cs.bounds();
4722
4723 Ok(AdaptiveMcReportV1 {
4724 schema_version: MC_ADAPTIVE_SCHEMA_VERSION_V1,
4725 method: MC_ADAPTIVE_METHOD_V1.to_string(),
4726 assumptions: MC_ADAPTIVE_ASSUMPTIONS_V1
4727 .iter()
4728 .map(|s| s.to_string())
4729 .collect(),
4730 confidence_percent: convergence.level.as_percent(),
4731 hit_probability: hits_cs.successes() as f64 / samples as f64,
4732 ci_low,
4733 ci_high,
4734 samples,
4735 attempts,
4736 arrivals: drop_at_target.count(),
4740 stop_reason,
4741 hit_radius_m,
4742 target_distance_m: sampler.target_distance,
4743 mean_impact_velocity_mps: impact_velocity.mean(),
4744 std_impact_velocity_mps: impact_velocity.sample_std(),
4745 mean_drop_at_target_m: drop_at_target.mean(),
4746 std_drop_at_target_m: drop_at_target.sample_std(),
4747 mean_wind_drift_at_target_m: drift_at_target.mean(),
4748 std_wind_drift_at_target_m: drift_at_target.sample_std(),
4749 })
4750}
4751
4752pub fn calculate_zero_angle(
4754 inputs: BallisticInputs,
4755 target_distance: f64,
4756 target_height: f64,
4757) -> Result<f64, BallisticsError> {
4758 calculate_zero_angle_with_conditions(
4759 inputs,
4760 target_distance,
4761 target_height,
4762 WindConditions::default(),
4763 AtmosphericConditions::default(),
4764 )
4765}
4766
4767pub fn calculate_zero_angle_with_conditions(
4768 inputs: BallisticInputs,
4769 target_distance: f64,
4770 target_height: f64,
4771 wind: WindConditions,
4772 atmosphere: AtmosphericConditions,
4773) -> Result<f64, BallisticsError> {
4774 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
4775 solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
4776}
4777
4778pub fn calculate_zero_angle_with_resolved_conditions(
4784 inputs: BallisticInputs,
4785 target_distance: f64,
4786 target_height: f64,
4787 wind: WindConditions,
4788 atmosphere: AtmosphericConditions,
4789) -> Result<f64, BallisticsError> {
4790 let mut solver =
4791 TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
4792 solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
4793}
4794
4795pub const ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M: f64 = 2000.0;
4803
4804pub fn calculate_zero_range_from_angle_with_conditions(
4819 inputs: BallisticInputs,
4820 zero_angle_rad: f64,
4821 target_height: f64,
4822 wind: WindConditions,
4823 atmosphere: AtmosphericConditions,
4824) -> Result<ZeroCrossings, BallisticsError> {
4825 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
4826 solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
4827 solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
4828}
4829
4830pub fn calculate_zero_range_from_angle_with_resolved_conditions(
4836 inputs: BallisticInputs,
4837 zero_angle_rad: f64,
4838 target_height: f64,
4839 wind: WindConditions,
4840 atmosphere: AtmosphericConditions,
4841) -> Result<ZeroCrossings, BallisticsError> {
4842 let mut solver =
4843 TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
4844 solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
4845 solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
4846}
4847
4848#[derive(Debug, Clone, Copy, PartialEq, Eq)]
4850pub enum BcFitMode {
4851 Drop,
4853 Velocity,
4856}
4857
4858#[derive(Debug, Clone, Copy)]
4860pub struct BcEstimate {
4861 pub bc: f64,
4863 pub rms_error: f64,
4865 pub drag_model: DragModel,
4867 pub mode: BcFitMode,
4869 pub at_bound: bool,
4873}
4874
4875fn fit_value_at(
4883 points: &[TrajectoryPoint],
4884 target_dist: f64,
4885 mode: BcFitMode,
4886 drop_offset: f64,
4887) -> Option<f64> {
4888 let val = |p: &TrajectoryPoint| match mode {
4889 BcFitMode::Drop => drop_offset - p.position.y,
4890 BcFitMode::Velocity => p.velocity_magnitude,
4891 };
4892 for i in 0..points.len() {
4893 if points[i].position.x >= target_dist {
4894 if i == 0 {
4895 return Some(val(&points[0]));
4896 }
4897 let p1 = &points[i - 1];
4898 let p2 = &points[i];
4899 let dx = p2.position.x - p1.position.x;
4900 if dx.abs() < 1e-9 {
4901 return Some(val(p2));
4902 }
4903 let t = (target_dist - p1.position.x) / dx;
4904 return Some(val(p1) + t * (val(p2) - val(p1)));
4905 }
4906 }
4907 None
4908}
4909
4910fn fit_residual_sse(
4911 trajectory: &[TrajectoryPoint],
4912 observations: &[(f64, f64)],
4913 mode: BcFitMode,
4914 drop_offset: f64,
4915) -> Option<f64> {
4916 if observations.is_empty() {
4917 return None;
4918 }
4919 let mut total = 0.0;
4920 for (target_dist, target_val) in observations {
4921 let value = fit_value_at(trajectory, *target_dist, mode, drop_offset)?;
4924 let error = value - target_val;
4925 total += error * error;
4926 }
4927 Some(total)
4928}
4929
4930#[allow(clippy::too_many_arguments)] pub fn estimate_bc_fit(
4948 velocity: f64,
4949 mass: f64,
4950 diameter: f64,
4951 points: &[(f64, f64)],
4952 drag_model: DragModel,
4953 mode: BcFitMode,
4954 atmosphere: AtmosphericConditions,
4955 zero_range: Option<f64>,
4956 sight_height: f64,
4957) -> Result<BcEstimate, BallisticsError> {
4958 if points.is_empty() {
4959 return Err(BallisticsError::from(
4960 "No data points provided for BC estimation.".to_string(),
4961 ));
4962 }
4963 let max_dist = points.iter().map(|(d, _)| *d).fold(0.0_f64, f64::max);
4964 let drop_offset = if zero_range.is_some() {
4967 sight_height
4968 } else {
4969 0.0
4970 };
4971
4972 let sse = |bc_value: f64| -> Option<f64> {
4974 let mut inputs = BallisticInputs {
4975 muzzle_velocity: velocity,
4976 bc_value,
4977 bc_type: drag_model,
4978 bullet_mass: mass,
4979 bullet_diameter: diameter,
4980 sight_height,
4981 ..Default::default()
4982 };
4983 if let Some(zr) = zero_range {
4986 let za = calculate_zero_angle_with_conditions(
4992 inputs.clone(),
4993 zr,
4994 sight_height,
4995 WindConditions::default(),
4996 atmosphere.clone(),
4997 )
4998 .ok()?;
4999 inputs.muzzle_angle = za;
5000 }
5001 let mut solver =
5002 TrajectorySolver::new(inputs, WindConditions::default(), atmosphere.clone());
5003 solver.set_max_range(max_dist * 1.5);
5004 let result = solver.solve().ok()?;
5005 fit_residual_sse(&result.points, points, mode, drop_offset)
5006 };
5007
5008 let (bc_min, bc_max) = match drag_model {
5012 DragModel::G7 => (0.05, 0.70),
5013 _ => (0.10, 1.20),
5014 };
5015
5016 let mut best_bc = f64::NAN;
5018 let mut best_sse = f64::MAX;
5019 let mut bc = bc_min;
5020 while bc <= bc_max + 1e-9 {
5021 if let Some(s) = sse(bc) {
5022 if s < best_sse {
5023 best_sse = s;
5024 best_bc = bc;
5025 }
5026 }
5027 bc += 0.01;
5028 }
5029 if !best_bc.is_finite() {
5030 return Err(BallisticsError::from(
5031 "Unable to estimate BC from provided data. Check that the values and units are correct."
5032 .to_string(),
5033 ));
5034 }
5035
5036 let lo = (best_bc - 0.01).max(bc_min);
5038 let hi = (best_bc + 0.01).min(bc_max);
5039 let mut bc = lo;
5040 while bc <= hi + 1e-9 {
5041 if let Some(s) = sse(bc) {
5042 if s < best_sse {
5043 best_sse = s;
5044 best_bc = bc;
5045 }
5046 }
5047 bc += 0.001;
5048 }
5049
5050 let at_bound = best_bc <= bc_min + 0.011 || best_bc >= bc_max - 0.011;
5053 let rms_error = (best_sse / points.len() as f64).sqrt();
5056 Ok(BcEstimate {
5057 bc: best_bc,
5058 rms_error,
5059 drag_model,
5060 mode,
5061 at_bound,
5062 })
5063}
5064
5065pub fn estimate_bc_from_trajectory(
5068 velocity: f64,
5069 mass: f64,
5070 diameter: f64,
5071 points: &[(f64, f64)], ) -> Result<f64, BallisticsError> {
5073 estimate_bc_fit(
5074 velocity,
5075 mass,
5076 diameter,
5077 points,
5078 DragModel::G1,
5079 BcFitMode::Drop,
5080 AtmosphericConditions::default(),
5081 None,
5082 0.05,
5083 )
5084 .map(|e| e.bc)
5085}
5086
5087use rand;
5089use rand_distr;
5090
5091#[cfg(test)]
5092mod mba737_powder_resolution_tests {
5093 use super::*;
5094
5095 #[test]
5096 fn linear_model_cold_powder_subtracts() {
5097 let v = resolve_powder_adjusted_velocity(823.0, 11.1, true, 0.5486, 21.1, None, None);
5099 assert!((v - (823.0 + 0.5486 * (11.1 - 21.1))).abs() < 1e-12);
5100 assert!(v < 823.0);
5101 }
5102
5103 #[test]
5104 fn linear_model_hot_powder_adds() {
5105 let v = resolve_powder_adjusted_velocity(823.0, 31.1, true, 0.5486, 21.1, None, None);
5106 assert!((v - (823.0 + 0.5486 * 10.0)).abs() < 1e-12);
5107 }
5108
5109 #[test]
5110 fn disabled_flag_is_passthrough() {
5111 let v = resolve_powder_adjusted_velocity(823.0, 40.0, false, 0.5486, 21.1, None, None);
5112 assert_eq!(v, 823.0);
5113 }
5114
5115 #[test]
5116 fn curve_overrides_linear_and_interpolates_at_powder_temp() {
5117 let curve = [(4.4, 798.6), (21.1, 823.0), (37.8, 841.2)];
5118 let v = resolve_powder_adjusted_velocity(
5120 823.0,
5121 30.0,
5122 true,
5123 99.0,
5124 21.1,
5125 Some(&curve),
5126 Some(4.4),
5127 );
5128 assert!((v - 798.6).abs() < 1e-9);
5129 }
5130
5131 #[test]
5132 fn curve_falls_back_to_ambient_and_clamps() {
5133 let curve = [(4.4, 798.6), (37.8, 841.2)];
5134 let v = resolve_powder_adjusted_velocity(823.0, -40.0, true, 1.0, 21.1, Some(&curve), None);
5136 assert!((v - 798.6).abs() < 1e-9);
5137 let v_hot =
5138 resolve_powder_adjusted_velocity(823.0, 60.0, true, 1.0, 21.1, Some(&curve), None);
5139 assert!((v_hot - 841.2).abs() < 1e-9);
5140 }
5141
5142 #[test]
5143 fn empty_curve_suppresses_linear_fallback() {
5144 let v = resolve_powder_adjusted_velocity(823.0, 40.0, true, 0.5486, 21.1, Some(&[]), None);
5146 assert_eq!(v, 823.0);
5147 }
5148
5149 #[test]
5150 fn sweep_huge_range_errors_instead_of_overflowing() {
5151 assert!(parse_powder_sweep("0:1e20:1").is_err());
5154 assert!(parse_powder_sweep("0:1e308:1e-3").is_err());
5155 }
5156
5157 #[test]
5158 fn sweep_fractional_step_keeps_end_row() {
5159 let rows = parse_powder_sweep("0:0.3:0.1").unwrap();
5161 assert_eq!(rows.len(), 4);
5162 assert!((rows[3] - 0.3).abs() < 1e-9);
5163 }
5164
5165 #[test]
5166 fn solver_and_helper_agree_on_linear_model() {
5167 let inputs = BallisticInputs {
5169 use_powder_sensitivity: true,
5170 powder_temp_sensitivity: 0.5486,
5171 powder_temp: 21.1,
5172 temperature: 4.4,
5173 ..Default::default()
5174 };
5175 let expected = resolve_powder_adjusted_velocity(
5176 inputs.muzzle_velocity,
5177 inputs.temperature,
5178 true,
5179 0.5486,
5180 21.1,
5181 None,
5182 None,
5183 );
5184 let solver = TrajectorySolver::new(
5185 inputs,
5186 WindConditions::default(),
5187 AtmosphericConditions::default(),
5188 );
5189 assert!((solver.inputs.muzzle_velocity - expected).abs() < 1e-12);
5190 }
5191}
5192
5193#[cfg(test)]
5194mod mba1302_solver_seam_tests {
5195 use super::*;
5196 use crate::wind::WindSegment;
5197
5198 #[test]
5199 fn authoritative_station_atmosphere_preserves_explicit_standard_values_at_altitude() {
5200 let atmosphere = AtmosphericConditions {
5201 temperature: 15.0,
5202 pressure: 1013.25,
5203 humidity: 50.0,
5204 altitude: 2_000.0,
5205 };
5206 let legacy = TrajectorySolver::new(
5207 BallisticInputs::default(),
5208 WindConditions::default(),
5209 atmosphere.clone(),
5210 );
5211 let authoritative = TrajectorySolver::new_with_resolved_station_atmosphere(
5212 BallisticInputs::default(),
5213 WindConditions::default(),
5214 atmosphere,
5215 );
5216
5217 let (legacy_density, _, legacy_temp_c, legacy_pressure_hpa) = legacy.resolved_atmosphere();
5218 let (authoritative_density, _, authoritative_temp_c, authoritative_pressure_hpa) =
5219 authoritative.resolved_atmosphere();
5220 let (icao_temp_k, icao_pressure_pa) =
5221 crate::atmosphere::calculate_icao_standard_atmosphere(2_000.0);
5222 let (expected_authoritative_density, _) =
5223 crate::atmosphere::calculate_atmosphere(2_000.0, Some(15.0), Some(1013.25), 50.0);
5224
5225 assert!((legacy_temp_c - (icao_temp_k - 273.15)).abs() < 1e-12);
5226 assert!((legacy_pressure_hpa - icao_pressure_pa / 100.0).abs() < 1e-12);
5227 assert_eq!(authoritative_temp_c.to_bits(), 15.0_f64.to_bits());
5228 assert_eq!(authoritative_pressure_hpa.to_bits(), 1013.25_f64.to_bits());
5229 assert_eq!(
5230 authoritative_density.to_bits(),
5231 expected_authoritative_density.to_bits()
5232 );
5233 assert!(
5234 (authoritative_density - legacy_density).abs() > 0.1,
5235 "explicit standard values at altitude must differ from ICAO-at-altitude: explicit={authoritative_density}, ICAO={legacy_density}"
5236 );
5237 }
5238
5239 #[test]
5248 fn precomputed_absolute_resolution_via_authoritative_matches_legacy_new() {
5249 for (temperature, pressure, altitude) in [
5250 (15.0, 1013.25, 0.0), (15.0, 1013.25, 2000.0), (-5.0, 850.0, 2000.0), (22.0, 950.0, 500.0),
5254 ] {
5255 let atmosphere = AtmosphericConditions {
5256 temperature,
5257 pressure,
5258 humidity: 50.0,
5259 altitude,
5260 };
5261 let legacy = TrajectorySolver::new(
5262 BallisticInputs::default(),
5263 WindConditions::default(),
5264 atmosphere.clone(),
5265 );
5266
5267 let (resolved_temp_c, resolved_pressure_hpa) =
5268 crate::atmosphere::resolve_station_conditions_with_pressure_mode(
5269 temperature,
5270 pressure,
5271 altitude,
5272 crate::atmosphere::PressureReferenceMode::Absolute,
5273 );
5274 let precomputed_atmosphere = AtmosphericConditions {
5275 temperature: resolved_temp_c,
5276 pressure: resolved_pressure_hpa,
5277 humidity: 50.0,
5278 altitude,
5279 };
5280 let precomputed = TrajectorySolver::new_with_resolved_station_atmosphere(
5281 BallisticInputs::default(),
5282 WindConditions::default(),
5283 precomputed_atmosphere,
5284 );
5285
5286 let (legacy_density, legacy_sos, legacy_temp_c, legacy_pressure_hpa) =
5287 legacy.resolved_atmosphere();
5288 let (pre_density, pre_sos, pre_temp_c, pre_pressure_hpa) =
5289 precomputed.resolved_atmosphere();
5290
5291 assert_eq!(
5292 legacy_temp_c.to_bits(),
5293 pre_temp_c.to_bits(),
5294 "temperature=({temperature}, {pressure}, {altitude})"
5295 );
5296 assert_eq!(
5297 legacy_pressure_hpa.to_bits(),
5298 pre_pressure_hpa.to_bits(),
5299 "pressure=({temperature}, {pressure}, {altitude})"
5300 );
5301 assert_eq!(legacy_density.to_bits(), pre_density.to_bits());
5302 assert_eq!(legacy_sos.to_bits(), pre_sos.to_bits());
5303 }
5304 }
5305
5306 fn configured_euler_zero(vertical_wind_mps: f64, time_step_s: f64) -> TrajectorySolver {
5307 let inputs = BallisticInputs {
5308 muzzle_velocity: 800.0,
5309 bc_value: 0.5,
5310 bc_type: DragModel::G7,
5311 bullet_mass: 0.0109,
5312 bullet_diameter: 0.00782,
5313 bullet_length: 0.0309,
5314 sight_height: 0.05,
5315 ground_threshold: -100.0,
5316 use_rk4: false,
5317 use_adaptive_rk45: false,
5318 ..BallisticInputs::default()
5319 };
5320 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(
5321 inputs,
5322 WindConditions::default(),
5323 AtmosphericConditions::default(),
5324 );
5325 solver.set_max_range(300.0);
5326 solver.set_time_step(time_step_s);
5327 if vertical_wind_mps != 0.0 {
5328 solver.set_wind_segments(vec![WindSegment {
5329 speed_kmh: 0.0,
5330 angle_deg: 0.0,
5331 until_m: 400.0,
5332 vertical_mps: vertical_wind_mps,
5333 }]);
5334 }
5335 solver
5336 }
5337
5338 #[test]
5339 fn inclined_shot_zeroes_like_a_level_rifle() {
5340 const ZERO_DISTANCE_M: f64 = 91.44; const SIGHT_HEIGHT_M: f64 = 0.0381; let inputs = BallisticInputs {
5349 bc_value: 0.5,
5350 bullet_mass: 150.0 * 0.06479891 / 1000.0,
5351 muzzle_velocity: 2700.0 * 0.3048,
5352 sight_height: SIGHT_HEIGHT_M,
5353 ..Default::default()
5354 };
5355
5356 let mut level = inputs.clone();
5357 level.shooting_angle = 0.0;
5358 let level_angle = TrajectorySolver::new(level, Default::default(), Default::default())
5359 .find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
5360 .expect("level zero must solve");
5361
5362 let mut inclined = inputs;
5363 inclined.shooting_angle = 5.71_f64.to_radians();
5364 let inclined_angle =
5365 TrajectorySolver::new(inclined, Default::default(), Default::default())
5366 .find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
5367 .expect("MBA-1412: a 5.71 deg incline at a 100 yd zero must be solvable");
5368
5369 assert!(
5370 (inclined_angle - level_angle).abs() < 1e-9,
5371 "zeroing is level-rifle sight geometry; incline must not move the solved zero: \
5372 level={level_angle}, inclined={inclined_angle}"
5373 );
5374 }
5375
5376 #[test]
5377 fn configured_zero_keeps_segments_method_and_time_step_then_sets_base_angle() {
5378 const TARGET_DISTANCE_M: f64 = 150.0;
5379 const TARGET_HEIGHT_M: f64 = 0.05;
5380
5381 let mut segmented = configured_euler_zero(-10.0, 0.02);
5384 let coarse_height = segmented
5385 .zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
5386 .expect("coarse configured trial")
5387 .expect("coarse trial reaches target");
5388 let mut fine = segmented.clone();
5389 fine.set_time_step(0.001);
5390 let fine_height = fine
5391 .zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
5392 .expect("fine configured trial")
5393 .expect("fine trial reaches target");
5394 assert!(
5395 (coarse_height - fine_height).abs() > 1e-5,
5396 "configured Euler step must affect zero trials: coarse={coarse_height}, fine={fine_height}"
5397 );
5398
5399 let segmented_angle = segmented
5400 .calculate_and_set_zero_angle(
5401 TARGET_DISTANCE_M,
5402 TARGET_HEIGHT_M,
5403 ZeroTargetFrame::SightLine,
5404 )
5405 .expect("segmented zero");
5406 assert_eq!(
5407 segmented.inputs.muzzle_angle.to_bits(),
5408 segmented_angle.to_bits(),
5409 "successful zero must install its angle on the configured solver"
5410 );
5411 assert_eq!(segmented.time_step.to_bits(), 0.02_f64.to_bits());
5412 assert_eq!(segmented.max_range.to_bits(), 300.0_f64.to_bits());
5413 assert!(segmented.wind_sock.is_some());
5414 assert_eq!(
5415 segmented.station_atmosphere_resolution,
5416 StationAtmosphereResolution::Authoritative
5417 );
5418 let zero_height = segmented
5419 .zero_trial_height_at(
5420 segmented_angle,
5421 TARGET_DISTANCE_M,
5422 ZeroTargetFrame::SightLine,
5423 )
5424 .expect("verify segmented zero")
5425 .expect("zeroed trial reaches target");
5426 assert!(
5427 (zero_height - TARGET_HEIGHT_M).abs() < 0.0001,
5428 "configured zero missed target: height={zero_height}"
5429 );
5430
5431 let mut calm = configured_euler_zero(0.0, 0.02);
5432 let calm_angle = calm
5433 .calculate_and_set_zero_angle(
5434 TARGET_DISTANCE_M,
5435 TARGET_HEIGHT_M,
5436 ZeroTargetFrame::SightLine,
5437 )
5438 .expect("calm zero");
5439 assert!(
5440 (segmented_angle - calm_angle).abs() > 1e-5,
5441 "segmented vertical wind must participate in zero trials: segmented={segmented_angle}, calm={calm_angle}"
5442 );
5443 }
5444}
5445
5446#[cfg(test)]
5447mod result_sanity_tests {
5448 use super::*;
5449
5450 fn default_solver() -> TrajectorySolver {
5451 TrajectorySolver::new(
5452 BallisticInputs::default(),
5453 WindConditions::default(),
5454 AtmosphericConditions::default(),
5455 )
5456 }
5457
5458 fn minimal_result() -> TrajectoryResult {
5459 TrajectoryResult {
5460 max_range: 100.0,
5461 max_height: 1.0,
5462 time_of_flight: 0.5,
5463 impact_velocity: 700.0,
5464 impact_energy: 2450.0,
5465 projectile_mass_kg: 0.01,
5466 line_of_sight_height_m: 1.5,
5467 station_speed_of_sound_mps: 340.0,
5468 termination: TrajectoryTermination::MaxRange,
5469 points: vec![],
5470 sampled_points: None,
5471 min_pitch_damping: None,
5472 transonic_mach: None,
5473 angular_state: None,
5474 max_yaw_angle: None,
5475 max_precession_angle: None,
5476 aerodynamic_jump: None,
5477 mach_1_2_distance_m: None,
5478 mach_1_0_distance_m: None,
5479 mach_0_9_distance_m: None,
5480 }
5481 }
5482
5483 #[test]
5484 fn mba1293_negative_scalars_fail_the_result_postcondition() {
5485 let solver = default_solver();
5486 solver
5487 .validate_result_sanity(&minimal_result())
5488 .expect("a sane result must pass");
5489
5490 for (name, mutate) in [
5491 (
5492 "max_range",
5493 (|r| r.max_range = -50.588) as fn(&mut TrajectoryResult),
5494 ),
5495 ("time_of_flight", |r| r.time_of_flight = -1.0),
5496 ("impact_velocity", |r| r.impact_velocity = -700.0),
5497 ("impact_energy", |r| r.impact_energy = -1.0),
5498 ] {
5499 let mut result = minimal_result();
5500 mutate(&mut result);
5501 let error = solver
5502 .validate_result_sanity(&result)
5503 .expect_err("negative scalar must fail");
5504 assert!(
5505 error.to_string().contains(name),
5506 "error for {name} did not name the field: {error}"
5507 );
5508 }
5509 }
5510
5511 #[test]
5512 fn mba1293_speed_budget_bounds_legitimate_states_and_rejects_divergence() {
5513 let solver = default_solver();
5514 let mv = solver.inputs.muzzle_velocity;
5515
5516 let position = Vector3::new(10.0, 0.0, 0.0);
5518 solver
5519 .validate_integration_state(&position, &Vector3::new(mv, 0.0, 0.0), 0.01)
5520 .expect("muzzle-speed state must pass");
5521
5522 let error = solver
5524 .validate_integration_state(&position, &Vector3::new(-13.0 * mv, 0.0, 0.0), 0.01)
5525 .expect_err("13x muzzle speed must fail the budget");
5526 assert!(error.to_string().contains("diverged"), "{error}");
5527
5528 let after_fall = mv + crate::constants::G_ACCEL_MPS2 * 60.0;
5530 solver
5531 .validate_integration_state(&position, &Vector3::new(0.0, -after_fall, 0.0), 60.0)
5532 .expect("gravity-accelerated speed within g*t must pass");
5533 }
5534}
5535
5536#[cfg(test)]
5537mod trajectory_point_budget_tests {
5538 use super::*;
5539 use crate::MAX_TRAJECTORY_SAMPLES;
5540
5541 fn solver_with_budget(
5542 use_rk4: bool,
5543 use_adaptive_rk45: bool,
5544 point_budget: usize,
5545 max_range: f64,
5546 ) -> TrajectorySolver {
5547 let inputs = BallisticInputs {
5548 use_rk4,
5549 use_adaptive_rk45,
5550 ground_threshold: f64::NEG_INFINITY,
5551 ..BallisticInputs::default()
5552 };
5553 let mut solver = TrajectorySolver::new(
5554 inputs,
5555 WindConditions::default(),
5556 AtmosphericConditions::default(),
5557 );
5558 solver.max_trajectory_points = point_budget;
5559 solver.set_max_range(max_range);
5560 solver.set_time_step(0.001);
5561 solver
5562 }
5563
5564 #[test]
5565 fn mba1283_every_solver_errors_instead_of_exceeding_point_budget() {
5566 for (mode, use_rk4, use_adaptive_rk45) in [
5567 ("Euler", false, false),
5568 ("RK4", true, false),
5569 ("RK45", true, true),
5570 ] {
5571 let error = solver_with_budget(use_rk4, use_adaptive_rk45, 3, 10.0)
5572 .solve()
5573 .expect_err("a solve requiring more than three points must fail");
5574 assert!(
5575 error.to_string().contains("point limit of 3"),
5576 "unexpected {mode} point-budget error: {error}"
5577 );
5578 }
5579 }
5580
5581 #[test]
5582 fn mba1283_interpolated_endpoint_counts_toward_point_budget() {
5583 for (mode, use_rk4, use_adaptive_rk45) in [
5584 ("Euler", false, false),
5585 ("RK4", true, false),
5586 ("RK45", true, true),
5587 ] {
5588 let result = solver_with_budget(use_rk4, use_adaptive_rk45, 2, 0.1)
5589 .solve()
5590 .expect("the initial point plus exact endpoint fit a two-point budget");
5591 assert_eq!(result.points.len(), 2, "unexpected {mode} point count");
5592
5593 let error = solver_with_budget(use_rk4, use_adaptive_rk45, 1, 0.1)
5594 .solve()
5595 .expect_err("the exact endpoint must not exceed a one-point budget");
5596 assert!(
5597 error.to_string().contains("point limit of 1"),
5598 "unexpected {mode} endpoint-budget error: {error}"
5599 );
5600 }
5601 }
5602
5603 #[test]
5611 fn mba1599_collapsed_stall_steps_count_toward_point_budget() {
5612 for (mode, use_rk4, use_adaptive_rk45, bc_value) in [
5613 ("Euler", false, false, 2.0e-4),
5614 ("RK4", true, false, 1.0e-4),
5615 ("RK45", true, true, 1.0e-4),
5616 ] {
5617 let stalled = |point_budget: usize| {
5618 let inputs = BallisticInputs {
5619 bc_value,
5620 bc_type: DragModel::G7,
5621 bullet_mass: 0.01134,
5622 bullet_diameter: 0.00782,
5623 muzzle_velocity: 823.0,
5624 use_rk4,
5625 use_adaptive_rk45,
5626 ground_threshold: -100.0,
5627 ..BallisticInputs::default()
5628 };
5629 let mut solver = TrajectorySolver::new(
5630 inputs,
5631 WindConditions::default(),
5632 AtmosphericConditions::default(),
5633 );
5634 solver.max_trajectory_points = point_budget;
5635 solver.set_max_range(300.0);
5636 solver.set_time_step(0.001);
5637 solver.solve()
5638 };
5639
5640 let flown = stalled(MAX_TRAJECTORY_POINTS)
5641 .unwrap_or_else(|error| panic!("{mode}: the stall ends inside the cap: {error}"));
5642 assert_eq!(
5643 flown.termination,
5644 TrajectoryTermination::GroundThreshold,
5645 "{mode}"
5646 );
5647 let stored = flown.points.len();
5648 let budget = stored + 1;
5649 let Err(error) = stalled(budget) else {
5650 panic!("{mode}: the flight fit a {budget}-point budget with {stored} stored points, so the collapsed steps were not counted");
5651 };
5652 assert!(
5653 error
5654 .to_string()
5655 .contains(&format!("trajectory point limit of {budget} exceeded")),
5656 "unexpected {mode} stall-budget error: {error}"
5657 );
5658 }
5659 }
5660
5661 #[test]
5665 fn mba1599_point_budget_fails_on_the_same_call_with_or_without_collapsing() {
5666 let solver = solver_with_budget(true, false, 4, 10.0);
5667 let point = |x: f64| TrajectoryPoint {
5668 time: x,
5669 position: Vector3::new(x, 0.0, 0.0),
5670 velocity_magnitude: 1.0,
5671 kinetic_energy: 0.0,
5672 drag_coefficient: None,
5673 };
5674 let assert_over_budget = |error: BallisticsError| {
5675 assert!(
5676 error
5677 .to_string()
5678 .contains("trajectory point limit of 4 exceeded"),
5679 "unexpected point-budget error: {error}"
5680 );
5681 };
5682 for distances in [
5683 [0.0, 1.0, 2.0, 3.0],
5684 [0.0, 1.0, 1.0, 2.0],
5685 [0.0, 0.0, 0.0, 0.0],
5686 ] {
5687 let mut points = Vec::new();
5688 let mut collapsed_steps = 0;
5689 for x in distances {
5690 solver
5691 .push_trajectory_point(&mut points, point(x), &mut collapsed_steps)
5692 .expect("four calls fit a four-point budget");
5693 }
5694 assert_eq!(points.len() + collapsed_steps, 4, "{distances:?}");
5695
5696 for x in [distances[3], 5.0] {
5698 let error = solver
5699 .push_trajectory_point(
5700 &mut points.clone(),
5701 point(x),
5702 &mut collapsed_steps.clone(),
5703 )
5704 .expect_err("a fifth call exceeds a four-point budget");
5705 assert_over_budget(error);
5706 }
5707 let error = solver
5709 .append_terminal_endpoint(
5710 &mut points.clone(),
5711 Vector3::new(20.0, 0.0, 0.0),
5712 Vector3::new(1.0, 0.0, 0.0),
5713 10.0,
5714 &mut 0.0,
5715 &mut collapsed_steps.clone(),
5716 )
5717 .expect_err("the terminal point is a fifth call");
5718 assert_over_budget(error);
5719 }
5720 }
5721
5722 #[test]
5723 fn mba1299_every_solver_preflights_the_sample_budget() {
5724 for (mode, use_rk4, use_adaptive_rk45) in [
5725 ("Euler", false, false),
5726 ("RK4", true, false),
5727 ("RK45", true, true),
5728 ] {
5729 let inputs = BallisticInputs {
5730 use_rk4,
5731 use_adaptive_rk45,
5732 enable_trajectory_sampling: true,
5733 sample_interval: 1.0,
5734 ground_threshold: f64::NEG_INFINITY,
5735 ..BallisticInputs::default()
5736 };
5737 let mut solver = TrajectorySolver::new(
5738 inputs,
5739 WindConditions::default(),
5740 AtmosphericConditions::default(),
5741 );
5742 solver.set_max_range(MAX_TRAJECTORY_SAMPLES as f64);
5743 solver.max_trajectory_points = 0;
5746
5747 let error = solver
5748 .solve()
5749 .expect_err("an over-limit sample grid must fail before integration");
5750 assert!(
5751 error
5752 .to_string()
5753 .contains("trajectory sample limit of 250000 exceeded"),
5754 "unexpected {mode} sample-budget error: {error}"
5755 );
5756 }
5757 }
5758
5759 #[test]
5760 fn mba1299_normal_sampling_does_not_change_solver_results() {
5761 for (mode, use_rk4, use_adaptive_rk45) in [
5762 ("Euler", false, false),
5763 ("RK4", true, false),
5764 ("RK45", true, true),
5765 ] {
5766 let solve = |enable_trajectory_sampling| {
5767 let inputs = BallisticInputs {
5768 use_rk4,
5769 use_adaptive_rk45,
5770 enable_trajectory_sampling,
5771 sample_interval: 0.5,
5772 ground_threshold: f64::NEG_INFINITY,
5773 ..BallisticInputs::default()
5774 };
5775 let mut solver = TrajectorySolver::new(
5776 inputs,
5777 WindConditions::default(),
5778 AtmosphericConditions::default(),
5779 );
5780 solver.set_max_range(2.0);
5781 solver.solve().expect("normal short-range solve")
5782 };
5783
5784 let baseline = solve(false);
5785 let sampled = solve(true);
5786 for (field, left, right) in [
5787 ("max_range", baseline.max_range, sampled.max_range),
5788 ("max_height", baseline.max_height, sampled.max_height),
5789 (
5790 "time_of_flight",
5791 baseline.time_of_flight,
5792 sampled.time_of_flight,
5793 ),
5794 (
5795 "impact_velocity",
5796 baseline.impact_velocity,
5797 sampled.impact_velocity,
5798 ),
5799 (
5800 "impact_energy",
5801 baseline.impact_energy,
5802 sampled.impact_energy,
5803 ),
5804 ] {
5805 assert_eq!(
5806 left.to_bits(),
5807 right.to_bits(),
5808 "{mode} sampling changed {field}"
5809 );
5810 }
5811 assert_eq!(baseline.points.len(), sampled.points.len());
5812 for (index, (left, right)) in baseline.points.iter().zip(&sampled.points).enumerate() {
5813 assert_eq!(
5814 left.time.to_bits(),
5815 right.time.to_bits(),
5816 "{mode} point {index}"
5817 );
5818 assert_eq!(
5819 left.position.map(f64::to_bits),
5820 right.position.map(f64::to_bits),
5821 "{mode} point {index} position"
5822 );
5823 assert_eq!(
5824 left.velocity_magnitude.to_bits(),
5825 right.velocity_magnitude.to_bits(),
5826 "{mode} point {index} velocity"
5827 );
5828 assert_eq!(
5829 left.kinetic_energy.to_bits(),
5830 right.kinetic_energy.to_bits(),
5831 "{mode} point {index} energy"
5832 );
5833 }
5834 assert!(baseline.sampled_points.is_none());
5835 let samples = sampled
5836 .sampled_points
5837 .expect("sampling-enabled solve should return observations");
5838 assert_eq!(
5839 samples
5840 .iter()
5841 .map(|sample| sample.distance_m)
5842 .collect::<Vec<_>>(),
5843 vec![0.0, 0.5, 1.0, 1.5, 2.0],
5844 "{mode} normal sampling grid changed"
5845 );
5846 }
5847 }
5848}
5849
5850#[cfg(test)]
5851mod monte_carlo_result_tests {
5852 use super::*;
5853
5854 fn make_results(impact_positions: Vec<Vector3<f64>>) -> MonteCarloResults {
5855 let count = impact_positions.len();
5856 MonteCarloResults {
5857 ranges: vec![500.0; count],
5858 impact_velocities: vec![300.0; count],
5859 impact_positions,
5860 }
5861 }
5862
5863 #[test]
5864 fn target_plane_cep_excludes_shortfall_markers() {
5865 let mut positions: Vec<Vector3<f64>> = (1..=5)
5866 .map(|radius| Vector3::new(0.0, radius as f64, 0.0))
5867 .collect();
5868 positions.extend((0..5).map(|_| Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)));
5869 let results = make_results(positions);
5870
5871 assert_eq!(results.target_arrival_count(), 5);
5872 assert_eq!(results.target_shortfall_fraction(), 0.5);
5873 assert_eq!(results.target_plane_cep(), Some(3.0));
5874
5875 let one_shortfall = make_results(vec![
5876 Vector3::new(0.0, 1.0, 0.0),
5877 Vector3::new(0.0, 2.0, 0.0),
5878 Vector3::new(0.0, 3.0, 0.0),
5879 Vector3::new(0.0, 4.0, 0.0),
5880 Vector3::new(0.0, 5.0, 0.0),
5881 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5882 ]);
5883 assert_eq!(one_shortfall.target_plane_cep(), Some(3.0));
5884 }
5885
5886 #[test]
5887 fn all_shortfalls_have_no_cep_but_still_count_as_misses() {
5888 let all_shortfalls = make_results(vec![
5889 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5890 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5891 ]);
5892 assert_eq!(all_shortfalls.target_arrival_count(), 0);
5893 assert_eq!(all_shortfalls.target_shortfall_fraction(), 1.0);
5894 assert_eq!(all_shortfalls.target_plane_cep(), None);
5895 assert_eq!(all_shortfalls.hit_probability(0.3), 0.0);
5896
5897 let one_hit_one_shortfall = make_results(vec![
5898 Vector3::new(0.0, 0.1, 0.0),
5899 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5900 ]);
5901 assert_eq!(one_hit_one_shortfall.hit_probability(0.3), 0.5);
5902 }
5903
5904 #[test]
5906 fn rect_hit_probability_checks_independent_axis_halves() {
5907 let results = make_results(vec![
5908 Vector3::new(0.0, 0.1, 0.1),
5910 Vector3::new(0.0, 0.0, 0.2),
5912 Vector3::new(0.0, 0.0, 0.201),
5914 Vector3::new(0.0, 0.301, 0.0),
5916 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5918 ]);
5919 assert!((results.rect_hit_probability(0.4, 0.6) - 0.4).abs() < 1e-12);
5921 }
5922
5923 #[test]
5924 fn rect_hit_probability_matches_circular_hit_probability_for_a_centered_hit() {
5925 let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
5926 assert_eq!(results.rect_hit_probability(0.5, 0.5), 1.0);
5927 assert_eq!(results.hit_probability(0.3), 1.0);
5928 }
5929
5930 #[test]
5931 fn rect_hit_probability_is_zero_for_empty_or_nonpositive_dimensions() {
5932 let empty = make_results(vec![]);
5933 assert_eq!(empty.rect_hit_probability(1.0, 1.0), 0.0);
5934
5935 let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
5936 assert_eq!(results.rect_hit_probability(0.0, 1.0), 0.0);
5937 assert_eq!(results.rect_hit_probability(1.0, 0.0), 0.0);
5938 assert_eq!(results.rect_hit_probability(-1.0, 1.0), 0.0);
5939 }
5940}
5941
5942#[cfg(test)]
5943mod monte_carlo_seeded_tests {
5944 use super::*;
5945
5946 fn seeded_test_fixture() -> (BallisticInputs, WindConditions) {
5955 (
5956 BallisticInputs {
5957 muzzle_velocity: 800.0,
5958 ..BallisticInputs::default()
5959 },
5960 WindConditions::default(),
5961 )
5962 }
5963
5964 fn loose_params() -> MonteCarloParams {
5974 MonteCarloParams {
5975 num_simulations: 1, velocity_std_dev: 3.0,
5977 angle_std_dev: 3.5e-4,
5978 bc_std_dev: 0.01,
5979 wind_speed_std_dev: 1.0,
5980 target_distance: Some(300.0),
5981 base_wind_speed: 0.0,
5982 base_wind_direction: 0.0,
5983 azimuth_std_dev: 3.5e-4,
5984 }
5985 }
5986
5987 fn mixed_arrival_params() -> MonteCarloParams {
5995 MonteCarloParams {
5996 num_simulations: 1,
5997 target_distance: Some(1920.0),
5998 ..MonteCarloParams::default()
5999 }
6000 }
6001
6002 #[test]
6018 fn legacy_seeded_estimates_are_pinned_bit_for_bit() {
6019 let (inputs, wind) = seeded_test_fixture();
6020 let params = MonteCarloParams {
6021 num_simulations: 200,
6022 target_distance: Some(500.0),
6023 ..MonteCarloParams::default()
6024 };
6025
6026 let results = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
6027 inputs,
6028 wind,
6029 params,
6030 0.01,
6031 0x1352_5EED,
6032 )
6033 .expect("seeded legacy run");
6034
6035 assert_eq!(results.ranges.len(), 200, "ranges length");
6040 assert_eq!(
6041 results.impact_velocities.len(),
6042 200,
6043 "impact_velocities length"
6044 );
6045 assert_eq!(
6046 results.impact_positions.len(),
6047 200,
6048 "impact_positions length"
6049 );
6050
6051 assert_eq!(
6053 results.hit_probability(DEFAULT_HIT_RADIUS_M).to_bits(),
6054 0.14_f64.to_bits(),
6055 "hit_probability = {:?}",
6056 results.hit_probability(DEFAULT_HIT_RADIUS_M)
6057 );
6058
6059 let expected_ranges: [f64; 3] = [1907.972891143359, 1936.408435469319, 1912.8150447617645];
6060 let expected_velocities: [f64; 3] =
6061 [238.6187151542299, 239.91923651600106, 243.14112455427164];
6062 let expected_positions: [(f64, f64, f64); 3] = [
6063 (0.0, -0.0643556039548101, 0.7344970252014579),
6064 (0.0, 0.5769422971539162, 0.27227201756386726),
6065 (0.0, -0.7440425792472842, 0.1541804446282822),
6066 ];
6067
6068 for (i, expected) in expected_ranges.iter().enumerate() {
6069 assert_eq!(
6070 results.ranges[i].to_bits(),
6071 expected.to_bits(),
6072 "ranges[{i}] = {:?}, pinned {expected:?}",
6073 results.ranges[i]
6074 );
6075 }
6076 for (i, expected) in expected_velocities.iter().enumerate() {
6077 assert_eq!(
6078 results.impact_velocities[i].to_bits(),
6079 expected.to_bits(),
6080 "impact_velocities[{i}] = {:?}, pinned {expected:?}",
6081 results.impact_velocities[i]
6082 );
6083 }
6084 for (i, (x, y, z)) in expected_positions.iter().enumerate() {
6085 let actual = results.impact_positions[i];
6086 assert_eq!(
6087 actual.x.to_bits(),
6088 x.to_bits(),
6089 "impact_positions[{i}].x = {:?}",
6090 actual.x
6091 );
6092 assert_eq!(
6093 actual.y.to_bits(),
6094 y.to_bits(),
6095 "impact_positions[{i}].y = {:?}",
6096 actual.y
6097 );
6098 assert_eq!(
6099 actual.z.to_bits(),
6100 z.to_bits(),
6101 "impact_positions[{i}].z = {:?}",
6102 actual.z
6103 );
6104 }
6105 }
6106
6107 #[test]
6108 fn seeded_runs_are_deterministic_and_match_the_using_rng_path() {
6109 let inputs = BallisticInputs {
6110 muzzle_velocity: 800.0,
6111 ..BallisticInputs::default()
6112 };
6113 let params = MonteCarloParams {
6114 num_simulations: 64,
6115 target_distance: Some(200.0),
6116 ..MonteCarloParams::default()
6117 };
6118
6119 let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
6120 inputs.clone(),
6121 WindConditions::default(),
6122 params.clone(),
6123 0.01,
6124 42,
6125 )
6126 .expect("seeded run a");
6127 let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
6128 inputs,
6129 WindConditions::default(),
6130 params,
6131 0.01,
6132 42,
6133 )
6134 .expect("seeded run b");
6135
6136 assert_eq!(a.ranges.len(), b.ranges.len());
6137 for (ra, rb) in a.ranges.iter().zip(b.ranges.iter()) {
6138 assert_eq!(ra.to_bits(), rb.to_bits());
6139 }
6140 for (pa, pb) in a.impact_positions.iter().zip(b.impact_positions.iter()) {
6141 assert_eq!(pa.x.to_bits(), pb.x.to_bits());
6142 assert_eq!(pa.y.to_bits(), pb.y.to_bits());
6143 assert_eq!(pa.z.to_bits(), pb.z.to_bits());
6144 }
6145 }
6146
6147 #[test]
6148 fn different_seeds_generally_produce_different_draws() {
6149 let inputs = BallisticInputs {
6150 muzzle_velocity: 800.0,
6151 ..BallisticInputs::default()
6152 };
6153 let params = MonteCarloParams {
6154 num_simulations: 32,
6155 velocity_std_dev: 5.0,
6156 target_distance: Some(200.0),
6157 ..MonteCarloParams::default()
6158 };
6159
6160 let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
6161 inputs.clone(),
6162 WindConditions::default(),
6163 params.clone(),
6164 0.0,
6165 1,
6166 )
6167 .expect("seeded run a");
6168 let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
6169 inputs,
6170 WindConditions::default(),
6171 params,
6172 0.0,
6173 2,
6174 )
6175 .expect("seeded run b");
6176
6177 assert_ne!(a.impact_velocities, b.impact_velocities);
6178 }
6179
6180 #[test]
6181 fn adaptive_stops_at_target_half_width_on_an_easy_case() {
6182 let (inputs, wind) = seeded_test_fixture();
6183 let conv = McConvergence {
6184 target_half_width: 0.05,
6185 ..Default::default()
6186 };
6187 let r = run_monte_carlo_adaptive_seeded(
6188 &inputs,
6189 &wind,
6190 &loose_params(),
6191 &conv,
6192 DEFAULT_HIT_RADIUS_M,
6193 0x1352_ADA9,
6194 )
6195 .unwrap();
6196
6197 assert_eq!(r.stop_reason, McStopReason::TargetHalfWidthMet);
6198 assert!(
6199 (r.ci_high - r.ci_low) / 2.0 <= 0.05 + 1e-12,
6200 "half-width {} exceeds the requested 0.05",
6201 (r.ci_high - r.ci_low) / 2.0
6202 );
6203 assert!(r.samples >= conv.min_samples, "stopped below min_samples");
6204 assert!(r.samples < conv.max_samples, "did not actually stop early");
6205 assert!(r.samples.is_multiple_of(conv.batch_size) || r.samples == conv.min_samples);
6206 assert!(r.ci_low <= r.hit_probability && r.hit_probability <= r.ci_high);
6207
6208 assert_eq!(r.hit_radius_m, DEFAULT_HIT_RADIUS_M);
6210 assert_eq!(r.target_distance_m, 300.0);
6211 assert_eq!(r.confidence_percent, 95);
6212 assert!(
6214 r.samples > 1,
6215 "params.num_simulations must be ignored by the adaptive driver"
6216 );
6217 assert!(
6220 r.mean_impact_velocity_mps > 0.0,
6221 "no impact velocity accumulated"
6222 );
6223 assert!(
6224 r.std_drop_at_target_m > 0.0 && r.std_wind_drift_at_target_m > 0.0,
6225 "dispersion collapsed: drop sd {} drift sd {}",
6226 r.std_drop_at_target_m,
6227 r.std_wind_drift_at_target_m
6228 );
6229 }
6230
6231 #[test]
6232 fn adaptive_caps_at_max_samples_on_an_impossible_target() {
6233 let (inputs, wind) = seeded_test_fixture();
6234 let conv = McConvergence {
6235 target_half_width: 1e-6,
6236 max_samples: 3_000,
6237 batch_size: 500,
6238 min_samples: 1_000,
6239 level: ConfidenceLevel::P95,
6240 };
6241 let r = run_monte_carlo_adaptive_seeded(
6242 &inputs,
6243 &wind,
6244 &loose_params(),
6245 &conv,
6246 DEFAULT_HIT_RADIUS_M,
6247 7,
6248 )
6249 .unwrap();
6250
6251 assert_eq!(r.stop_reason, McStopReason::MaxSamplesReached);
6252 assert_eq!(r.samples, 3_000);
6253 }
6254
6255 #[test]
6264 fn adaptive_stops_between_the_floor_and_the_ceiling() {
6265 let (inputs, wind) = seeded_test_fixture();
6266 let conv = McConvergence {
6267 level: ConfidenceLevel::P95,
6268 target_half_width: 0.03,
6269 min_samples: 0,
6270 max_samples: 10_000,
6271 batch_size: 100,
6272 };
6273 let r = run_monte_carlo_adaptive_seeded(
6274 &inputs,
6275 &wind,
6276 &loose_params(),
6277 &conv,
6278 DEFAULT_HIT_RADIUS_M,
6279 0x1352_5A1D,
6280 )
6281 .unwrap();
6282
6283 assert_eq!(r.stop_reason, McStopReason::TargetHalfWidthMet);
6284 assert!(r.samples > 0);
6285 assert!(
6286 r.samples.is_multiple_of(conv.batch_size),
6287 "samples {} is not a whole number of batches",
6288 r.samples
6289 );
6290 assert!(
6291 r.samples > conv.min_samples,
6292 "stopped on the floor, not on the data"
6293 );
6294 assert!(
6295 r.samples < conv.max_samples,
6296 "ran to the ceiling, so nothing adaptive was exercised"
6297 );
6298 assert!(
6299 r.samples > conv.batch_size,
6300 "stopped on the very first look ({} samples); the multi-batch path is untested",
6301 r.samples
6302 );
6303 assert!((r.ci_high - r.ci_low) / 2.0 <= 0.03 + 1e-12);
6304 assert!(r.ci_low <= r.hit_probability && r.hit_probability <= r.ci_high);
6305 assert_eq!(r.attempts, r.samples, "no trial should have been dropped");
6306 }
6307
6308 #[test]
6316 fn adaptive_runs_a_truncated_final_batch_up_to_max_samples() {
6317 let (inputs, wind) = seeded_test_fixture();
6318 let conv = McConvergence {
6319 level: ConfidenceLevel::P95,
6320 target_half_width: 1e-6,
6321 min_samples: 0,
6322 max_samples: 750,
6323 batch_size: 500,
6324 };
6325 let r = run_monte_carlo_adaptive_seeded(
6326 &inputs,
6327 &wind,
6328 &loose_params(),
6329 &conv,
6330 DEFAULT_HIT_RADIUS_M,
6331 0x1352_7B10,
6332 )
6333 .unwrap();
6334
6335 assert_eq!(r.stop_reason, McStopReason::MaxSamplesReached);
6336 assert_eq!(
6337 r.samples, 750,
6338 "the 250-trial final batch did not run, or was not truncated"
6339 );
6340 assert_eq!(r.attempts, 750);
6341 assert!(!r.samples.is_multiple_of(conv.batch_size));
6342 }
6343
6344 #[test]
6345 fn adaptive_is_deterministic_for_a_seed() {
6346 let (inputs, wind) = seeded_test_fixture();
6347 let conv = McConvergence::default();
6348 let a = run_monte_carlo_adaptive_seeded(
6349 &inputs,
6350 &wind,
6351 &loose_params(),
6352 &conv,
6353 DEFAULT_HIT_RADIUS_M,
6354 99,
6355 )
6356 .unwrap();
6357 let b = run_monte_carlo_adaptive_seeded(
6358 &inputs,
6359 &wind,
6360 &loose_params(),
6361 &conv,
6362 DEFAULT_HIT_RADIUS_M,
6363 99,
6364 )
6365 .unwrap();
6366
6367 assert_eq!(a.hit_probability.to_bits(), b.hit_probability.to_bits());
6368 assert_eq!(a.samples, b.samples);
6369 assert_eq!(a.ci_low.to_bits(), b.ci_low.to_bits());
6370 assert_eq!(a.ci_high.to_bits(), b.ci_high.to_bits());
6371 assert_eq!(
6374 a.mean_impact_velocity_mps.to_bits(),
6375 b.mean_impact_velocity_mps.to_bits()
6376 );
6377 assert_eq!(
6378 a.std_drop_at_target_m.to_bits(),
6379 b.std_drop_at_target_m.to_bits()
6380 );
6381 }
6382
6383 #[test]
6384 fn adaptive_report_carries_schema_method_and_all_four_assumptions() {
6385 let (inputs, wind) = seeded_test_fixture();
6386 let conv = McConvergence {
6390 min_samples: 0,
6391 max_samples: 50,
6392 batch_size: 50,
6393 target_half_width: 1.0,
6394 level: ConfidenceLevel::P90,
6395 };
6396 let r = run_monte_carlo_adaptive_seeded(
6397 &inputs,
6398 &wind,
6399 &mixed_arrival_params(),
6400 &conv,
6401 DEFAULT_HIT_RADIUS_M,
6402 0x1352_D0C5,
6403 )
6404 .unwrap();
6405
6406 assert_eq!(r.schema_version, MC_ADAPTIVE_SCHEMA_VERSION_V1);
6407 assert_eq!(r.schema_version, 1);
6408 assert_eq!(r.method, "anytime_beta_binomial_mixture_cs_v1");
6409 assert_eq!(r.confidence_percent, 90);
6410
6411 assert_eq!(r.attempts, 50, "one full batch was drawn");
6417 assert_eq!(r.samples, 50, "no trial was dropped by the solver");
6418 assert!(
6419 r.arrivals > 0 && r.arrivals < r.samples,
6420 "fixture must split the run: arrivals {} of samples {}",
6421 r.arrivals,
6422 r.samples
6423 );
6424 assert!(
6427 r.arrivals >= 2,
6428 "arrivals {} too few for a sample sd",
6429 r.arrivals
6430 );
6431 assert!(r.std_drop_at_target_m > 0.0 && r.std_impact_velocity_mps > 0.0);
6432 assert!(r.attempts >= r.samples && r.samples >= r.arrivals);
6434
6435 assert_eq!(r.assumptions.len(), 4, "exactly four assumptions expected");
6441 assert_eq!(
6442 r.assumptions[0],
6443 "Sampling uncertainty only: intervals cover Monte Carlo sampling error, not model error in the trajectory solver or its inputs."
6444 );
6445 assert_eq!(
6446 r.assumptions[1],
6447 "Anytime-valid stopping: the beta-binomial mixture confidence sequence keeps its coverage guarantee despite stopping the moment the target half-width is met."
6448 );
6449 assert_eq!(
6450 r.assumptions[2],
6451 "Input dispersions are the independent normal distributions declared in MonteCarloParams; correlations between inputs are not modeled."
6452 );
6453 assert_eq!(
6454 r.assumptions[3],
6455 "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."
6456 );
6457
6458 assert_eq!(
6460 serde_json::to_string(&McStopReason::TargetHalfWidthMet).unwrap(),
6461 "\"target_half_width_met\""
6462 );
6463 assert_eq!(
6464 serde_json::to_string(&McStopReason::MaxSamplesReached).unwrap(),
6465 "\"max_samples_reached\""
6466 );
6467 }
6468
6469 #[test]
6470 fn adaptive_rejects_nonsense_convergence() {
6471 let (inputs, wind) = seeded_test_fixture();
6472 let run = |conv: McConvergence| {
6473 run_monte_carlo_adaptive_seeded(
6474 &inputs,
6475 &wind,
6476 &loose_params(),
6477 &conv,
6478 DEFAULT_HIT_RADIUS_M,
6479 1,
6480 )
6481 .unwrap_err()
6482 };
6483
6484 for bad_width in [0.0, -0.01, f64::NAN] {
6485 let err = run(McConvergence {
6486 target_half_width: bad_width,
6487 ..Default::default()
6488 });
6489 assert!(
6490 err.contains("target_half_width"),
6491 "error must name the field, got: {err}"
6492 );
6493 }
6494
6495 let err = run(McConvergence {
6496 batch_size: 0,
6497 ..Default::default()
6498 });
6499 assert!(err.contains("batch_size"), "got: {err}");
6500
6501 let err = run(McConvergence {
6502 min_samples: 5_000,
6503 max_samples: 1_000,
6504 ..Default::default()
6505 });
6506 assert!(err.contains("max_samples"), "got: {err}");
6507 assert!(err.contains("min_samples"), "got: {err}");
6508
6509 let err = run(McConvergence {
6510 max_samples: 0,
6511 min_samples: 0,
6512 ..Default::default()
6513 });
6514 assert!(err.contains("max_samples"), "got: {err}");
6515
6516 let err = McConvergence {
6518 batch_size: 0,
6519 ..Default::default()
6520 }
6521 .validate()
6522 .unwrap_err();
6523 assert!(err.contains("batch_size"));
6524 }
6525
6526 #[test]
6527 fn wilson_companion_matches_hit_probability_and_wilson_interval() {
6528 let (inputs, wind) = seeded_test_fixture();
6529 let params = MonteCarloParams {
6530 num_simulations: 128,
6531 target_distance: Some(500.0),
6532 ..MonteCarloParams::default()
6533 };
6534 let results = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
6535 inputs,
6536 wind,
6537 params,
6538 0.01,
6539 0x1352_C0DE,
6540 )
6541 .expect("seeded legacy run");
6542
6543 for level in [
6544 ConfidenceLevel::P90,
6545 ConfidenceLevel::P95,
6546 ConfidenceLevel::P99,
6547 ] {
6548 let (p_hat, (lo, hi), n) = results.hit_probability_wilson(DEFAULT_HIT_RADIUS_M, level);
6549
6550 assert_eq!(
6553 p_hat.to_bits(),
6554 results.hit_probability(DEFAULT_HIT_RADIUS_M).to_bits()
6555 );
6556 assert_eq!(n, results.impact_positions.len() as u64);
6557
6558 let hits = results
6562 .impact_positions
6563 .iter()
6564 .filter(|p| MonteCarloResults::position_is_hit(p, DEFAULT_HIT_RADIUS_M))
6565 .count() as u64;
6566 let (want_lo, want_hi) = wilson_interval(hits, n, level);
6567 assert_eq!(lo.to_bits(), want_lo.to_bits());
6568 assert_eq!(hi.to_bits(), want_hi.to_bits());
6569 assert!(lo <= p_hat && p_hat <= hi, "interval excludes p_hat");
6570 }
6571
6572 let empty = MonteCarloResults {
6574 ranges: Vec::new(),
6575 impact_velocities: Vec::new(),
6576 impact_positions: Vec::new(),
6577 };
6578 assert_eq!(
6579 empty.hit_probability_wilson(DEFAULT_HIT_RADIUS_M, ConfidenceLevel::P95),
6580 (0.0, (0.0, 1.0), 0)
6581 );
6582 }
6583}
6584
6585#[cfg(test)]
6586mod monte_carlo_powder_curve_tests {
6587 use super::*;
6588 use rand::{rngs::StdRng, SeedableRng};
6589
6590 #[test]
6591 fn powder_curve_preserves_sampled_muzzle_velocity_dispersion() {
6592 let inputs = BallisticInputs {
6593 muzzle_velocity: 700.0,
6594 powder_temp_curve: Some(vec![(15.0, 800.0)]),
6595 powder_curve_temp_c: Some(15.0),
6596 ..BallisticInputs::default()
6597 };
6598 let params = MonteCarloParams {
6599 num_simulations: 16,
6600 velocity_std_dev: 20.0,
6601 angle_std_dev: 1e-12,
6602 bc_std_dev: 1e-12,
6603 wind_speed_std_dev: 1e-12,
6604 target_distance: Some(100.0),
6605 azimuth_std_dev: 1e-12,
6606 ..MonteCarloParams::default()
6607 };
6608
6609 let mut rng = StdRng::seed_from_u64(0x5EED_1176);
6610 let results = run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
6611 inputs,
6612 WindConditions::default(),
6613 params,
6614 0.0,
6615 &mut rng,
6616 )
6617 .expect("Monte Carlo solve");
6618 let min_velocity = results
6619 .impact_velocities
6620 .iter()
6621 .copied()
6622 .fold(f64::INFINITY, f64::min);
6623 let max_velocity = results
6624 .impact_velocities
6625 .iter()
6626 .copied()
6627 .fold(f64::NEG_INFINITY, f64::max);
6628
6629 assert!(
6630 max_velocity - min_velocity > 1.0,
6631 "20 m/s muzzle spread collapsed after curve resolution: impact-velocity span={} m/s",
6632 max_velocity - min_velocity
6633 );
6634 }
6635}
6636
6637#[cfg(test)]
6638mod monte_carlo_wind_sampling_tests {
6639 use super::*;
6640 use rand::{rngs::StdRng, SeedableRng};
6641
6642 #[test]
6643 fn wind_speed_sigma_does_not_change_seeded_direction_draws() {
6644 let base_wind = WindConditions {
6645 speed: 100.0,
6646 direction: 0.37,
6647 vertical_speed: 0.0,
6648 };
6649 let narrow_speed = MonteCarloWindSampler::new(&base_wind, 0.5, 0.2).unwrap();
6650 let wide_speed = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
6651 let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
6652 let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
6653 let mut speed_changed = false;
6654
6655 for _ in 0..32 {
6656 let narrow = narrow_speed.sample(&mut narrow_rng);
6657 let wide = wide_speed.sample(&mut wide_rng);
6658 assert!(narrow.speed > 0.0 && wide.speed > 0.0);
6659 assert_eq!(narrow.direction.to_bits(), wide.direction.to_bits());
6660 speed_changed |= narrow.speed.to_bits() != wide.speed.to_bits();
6661 }
6662 assert!(
6663 speed_changed,
6664 "different speed sigmas must still vary speed draws"
6665 );
6666 }
6667
6668 #[test]
6669 fn zero_direction_sigma_has_no_angular_jitter() {
6670 let base_wind = WindConditions {
6671 speed: 100.0,
6672 direction: 0.37,
6673 vertical_speed: 0.0,
6674 };
6675 let sampler = MonteCarloWindSampler::new(&base_wind, 4.0, 0.0).unwrap();
6676 let mut rng = StdRng::seed_from_u64(0x5EED_1223);
6677 let mut speed_changed = false;
6678
6679 for _ in 0..32 {
6680 let wind = sampler.sample(&mut rng);
6681 speed_changed |= wind.speed.to_bits() != base_wind.speed.to_bits();
6682 assert_eq!(wind.direction.to_bits(), base_wind.direction.to_bits());
6683 }
6684 assert!(speed_changed, "speed uncertainty should remain active");
6685 }
6686
6687 #[test]
6688 fn direction_sigma_controls_seeded_angular_spread_in_radians() {
6689 let base_wind = WindConditions {
6690 speed: 100.0,
6691 direction: 0.37,
6692 vertical_speed: 0.0,
6693 };
6694 let narrow = MonteCarloWindSampler::new(&base_wind, 4.0, 0.1).unwrap();
6695 let wide = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
6696 let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
6697 let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
6698 let mut nonzero_direction_draw = false;
6699
6700 for _ in 0..32 {
6701 let narrow_wind = narrow.sample(&mut narrow_rng);
6702 let wide_wind = wide.sample(&mut wide_rng);
6703 assert_eq!(narrow_wind.speed.to_bits(), wide_wind.speed.to_bits());
6704
6705 let narrow_delta = narrow_wind.direction - base_wind.direction;
6706 let wide_delta = wide_wind.direction - base_wind.direction;
6707 assert!((wide_delta - 2.0 * narrow_delta).abs() < 1e-12);
6708 nonzero_direction_draw |= narrow_delta.abs() > 1e-6;
6709 }
6710 assert!(
6711 nonzero_direction_draw,
6712 "positive radians sigma must vary direction"
6713 );
6714 }
6715
6716 #[test]
6717 fn direction_sigma_rejects_negative_or_nonfinite_values() {
6718 let base_wind = WindConditions::default();
6719 for sigma in [-0.1, f64::NAN, f64::INFINITY] {
6720 assert!(MonteCarloWindSampler::new(&base_wind, 1.0, sigma).is_err());
6721 }
6722 }
6723
6724 #[test]
6725 fn base_vertical_wind_rides_into_every_mc_sample() {
6726 use rand::SeedableRng;
6730 let base_wind = WindConditions {
6731 vertical_speed: 4.2,
6732 ..Default::default()
6733 };
6734 let sampler = MonteCarloWindSampler::new(&base_wind, 1.0, 0.2).unwrap();
6735 let mut rng = rand::rngs::StdRng::seed_from_u64(7);
6736 for _ in 0..32 {
6737 let w = sampler.sample(&mut rng);
6738 assert_eq!(w.vertical_speed, 4.2);
6739 }
6740 }
6741
6742 #[test]
6743 fn negative_speed_sample_reverses_wind_direction() {
6744 let direction = 0.25;
6745 let signed_speed = -2.5;
6746 let wind = wind_from_signed_speed_sample(signed_speed, direction, 0.0);
6747 let positive_wind = wind_from_signed_speed_sample(2.5, direction, 0.0);
6748
6749 assert_eq!(wind.speed, 2.5);
6750 assert!(
6751 (wind.direction - (direction + std::f64::consts::PI)).abs() < f64::EPSILON,
6752 "negative speed must reverse direction by pi: got {}",
6753 wind.direction
6754 );
6755 assert_eq!(positive_wind.speed, 2.5);
6756 assert_eq!(positive_wind.direction, direction);
6757
6758 let normalized_x = -wind.speed * wind.direction.cos();
6759 let normalized_z = -wind.speed * wind.direction.sin();
6760 let signed_x = -signed_speed * direction.cos();
6761 let signed_z = -signed_speed * direction.sin();
6762 assert!((normalized_x - signed_x).abs() < 1e-12);
6763 assert!((normalized_z - signed_z).abs() < 1e-12);
6764 }
6765}
6766
6767#[cfg(test)]
6768mod bc_fit_objective_tests {
6769 use super::*;
6770
6771 fn velocity_point(range_m: f64, velocity_mps: f64) -> TrajectoryPoint {
6772 TrajectoryPoint {
6773 time: 0.0,
6774 position: Vector3::new(range_m, 0.0, 0.0),
6775 velocity_magnitude: velocity_mps,
6776 kinetic_energy: 0.0,
6777 drag_coefficient: None,
6778 }
6779 }
6780
6781 #[test]
6782 fn candidate_that_misses_an_observation_has_no_score() {
6783 let trajectory = vec![velocity_point(0.0, 800.0), velocity_point(100.0, 700.0)];
6784 let observations = vec![(50.0, 750.0), (150.0, 600.0)];
6785
6786 assert!(
6787 fit_residual_sse(&trajectory, &observations, BcFitMode::Velocity, 0.0).is_none(),
6788 "a candidate that reaches only one of two observations must not compete on partial SSE"
6789 );
6790
6791 let complete_observations = vec![(50.0, 740.0), (100.0, 680.0)];
6792 assert_eq!(
6793 fit_residual_sse(
6794 &trajectory,
6795 &complete_observations,
6796 BcFitMode::Velocity,
6797 0.0,
6798 ),
6799 Some(500.0)
6800 );
6801 }
6802}
6803
6804#[cfg(test)]
6805mod cluster_bc_reference_space_tests {
6806 use super::*;
6807
6808 fn acceleration_at_1100_fps(inputs: BallisticInputs) -> Vector3<f64> {
6809 let solver = TrajectorySolver::new(
6810 inputs,
6811 WindConditions::default(),
6812 AtmosphericConditions::default(),
6813 );
6814 let position = Vector3::zeros();
6815 let velocity = Vector3::new(1100.0 / 3.28084, 0.0, 0.0);
6816 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6817 solver.calculate_acceleration(
6818 &position,
6819 &velocity,
6820 &Vector3::zeros(),
6821 (temp_c, pressure_hpa, density / 1.225),
6822 )
6823 }
6824
6825 #[test]
6826 fn solver_passes_g7_reference_model_to_cluster_classifier() {
6827 let inputs = BallisticInputs {
6828 bc_value: 0.190,
6829 bc_type: DragModel::G7,
6830 bullet_mass: 77.0 * crate::constants::GRAINS_TO_KG,
6831 bullet_diameter: 0.224 * 0.0254,
6832 use_cluster_bc: true,
6833 ..BallisticInputs::default()
6834 };
6835
6836 let solver = TrajectorySolver::new(
6837 inputs,
6838 WindConditions::default(),
6839 AtmosphericConditions::default(),
6840 );
6841 let corrected = solver.apply_cluster_bc_correction(0.190, 2800.0);
6842
6843 assert!(
6844 (corrected / 0.190 - 1.004).abs() < 1e-12,
6845 "solver selected the wrong G7 cluster multiplier: {}",
6846 corrected / 0.190
6847 );
6848 }
6849
6850 #[test]
6851 fn velocity_bc_segments_are_not_cluster_corrected_twice() {
6852 let segmented_clustered = BallisticInputs {
6853 bc_value: 0.5,
6854 bc_type: DragModel::G7,
6855 use_bc_segments: true,
6856 bc_segments_data: Some(vec![
6857 crate::BCSegmentData {
6858 velocity_min: 0.0,
6859 velocity_max: 1_600.0,
6860 bc_value: 0.4,
6861 },
6862 crate::BCSegmentData {
6863 velocity_min: 1_600.0,
6864 velocity_max: 5_000.0,
6865 bc_value: 0.45,
6866 },
6867 ]),
6868 use_cluster_bc: true,
6869 ..BallisticInputs::default()
6870 };
6871 let mut segmented_only = segmented_clustered.clone();
6872 segmented_only.use_cluster_bc = false;
6873 let mut constant_clustered = segmented_clustered.clone();
6874 constant_clustered.bc_value = 0.4;
6875 constant_clustered.bc_segments_data = None;
6876
6877 let stacked = acceleration_at_1100_fps(segmented_clustered);
6878 let segment_only = acceleration_at_1100_fps(segmented_only);
6879 let cluster_only = acceleration_at_1100_fps(constant_clustered);
6880
6881 assert!(
6882 (stacked.x - segment_only.x).abs() < 1e-12,
6883 "segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
6884 stacked.x,
6885 segment_only.x
6886 );
6887 assert!(
6888 (cluster_only.x - segment_only.x).abs() > 1e-6,
6889 "cluster correction must remain active for a constant BC"
6890 );
6891 }
6892
6893 #[test]
6894 fn mach_bc_segments_are_not_cluster_corrected_twice() {
6895 let mach_segmented_clustered = BallisticInputs {
6896 bc_value: 0.5,
6897 bc_type: DragModel::G7,
6898 use_bc_segments: false,
6899 bc_segments: Some(vec![(0.5, 0.3), (1.5, 0.5)]),
6900 use_cluster_bc: true,
6901 ..BallisticInputs::default()
6902 };
6903 let mut mach_segmented_only = mach_segmented_clustered.clone();
6904 mach_segmented_only.use_cluster_bc = false;
6905
6906 let stacked = acceleration_at_1100_fps(mach_segmented_clustered);
6907 let segment_only = acceleration_at_1100_fps(mach_segmented_only);
6908
6909 assert!(
6910 (stacked.x - segment_only.x).abs() < 1e-12,
6911 "Mach segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
6912 stacked.x,
6913 segment_only.x
6914 );
6915 }
6916}
6917
6918#[cfg(test)]
6919mod velocity_bc_flag_tests {
6920 use super::*;
6921
6922 fn acceleration_at_600_mps(inputs: BallisticInputs) -> Vector3<f64> {
6923 let solver = TrajectorySolver::new(
6924 inputs,
6925 WindConditions::default(),
6926 AtmosphericConditions::default(),
6927 );
6928 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6929 solver.calculate_acceleration(
6930 &Vector3::zeros(),
6931 &Vector3::new(600.0, 0.0, 0.0),
6932 &Vector3::zeros(),
6933 (temp_c, pressure_hpa, density / 1.225),
6934 )
6935 }
6936
6937 #[test]
6938 fn velocity_bc_data_requires_opt_in_in_trajectory_solver() {
6939 let scalar_inputs = BallisticInputs {
6940 bc_value: 0.5,
6941 bc_type: DragModel::G7,
6942 ..BallisticInputs::default()
6943 };
6944 let mut disabled_inputs = scalar_inputs.clone();
6945 disabled_inputs.bc_segments_data = Some(vec![crate::BCSegmentData {
6946 velocity_min: 0.0,
6947 velocity_max: 4_000.0,
6948 bc_value: 0.46,
6949 }]);
6950 disabled_inputs.use_bc_segments = false;
6951 let mut enabled_inputs = disabled_inputs.clone();
6952 enabled_inputs.use_bc_segments = true;
6953 let mut mach_only_inputs = scalar_inputs.clone();
6954 mach_only_inputs.bc_segments = Some(vec![(0.0, 0.4), (3.0, 0.4)]);
6955 let mut disabled_with_both = mach_only_inputs.clone();
6956 disabled_with_both.bc_segments_data = disabled_inputs.bc_segments_data.clone();
6957
6958 let scalar = acceleration_at_600_mps(scalar_inputs);
6959 let disabled = acceleration_at_600_mps(disabled_inputs);
6960 let enabled = acceleration_at_600_mps(enabled_inputs);
6961 let mach_only = acceleration_at_600_mps(mach_only_inputs);
6962 let disabled_with_both = acceleration_at_600_mps(disabled_with_both);
6963
6964 assert_eq!(
6965 disabled.x.to_bits(),
6966 scalar.x.to_bits(),
6967 "a populated velocity table must not change drag while use_bc_segments is false"
6968 );
6969 assert!(
6970 enabled.x < disabled.x - 1.0,
6971 "enabling the lower BC table must increase drag: disabled ax={} enabled ax={}",
6972 disabled.x,
6973 enabled.x
6974 );
6975 assert_eq!(
6976 disabled_with_both.x.to_bits(),
6977 mach_only.x.to_bits(),
6978 "disabling velocity data must fall through to an explicit Mach table"
6979 );
6980 }
6981}
6982
6983#[cfg(test)]
6984mod mach_bc_segment_tests {
6985 use super::*;
6986
6987 #[test]
6988 fn trajectory_solver_interpolates_explicit_mach_bc_segments() {
6989 let segmented_inputs = BallisticInputs {
6990 bc_value: 0.8,
6991 use_bc_segments: false,
6992 bc_segments: Some(vec![(1.0, 0.2), (2.0, 0.4)]),
6993 bc_segments_data: None,
6994 ..BallisticInputs::default()
6995 };
6996
6997 let mut expected_inputs = segmented_inputs.clone();
6998 expected_inputs.bc_value = 0.3;
6999 expected_inputs.bc_segments = None;
7000
7001 let atmosphere = AtmosphericConditions::default();
7002 let segmented_solver = TrajectorySolver::new(
7003 segmented_inputs,
7004 WindConditions::default(),
7005 atmosphere.clone(),
7006 );
7007 let expected_solver =
7008 TrajectorySolver::new(expected_inputs, WindConditions::default(), atmosphere);
7009 let position = Vector3::zeros();
7010 let (density, _, temp_c, pressure_hpa) = segmented_solver.resolved_atmosphere();
7011 let (_, local_speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
7012 segmented_solver.atmosphere.altitude,
7013 segmented_solver.atmosphere.altitude,
7014 temp_c,
7015 pressure_hpa,
7016 density / 1.225,
7017 segmented_solver.atmosphere.humidity,
7018 );
7019 let velocity = Vector3::new(1.5 * local_speed_of_sound, 0.0, 0.0);
7020 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
7021
7022 let segmented_acceleration = segmented_solver.calculate_acceleration(
7023 &position,
7024 &velocity,
7025 &Vector3::zeros(),
7026 resolved_atmo,
7027 );
7028 let expected_acceleration = expected_solver.calculate_acceleration(
7029 &position,
7030 &velocity,
7031 &Vector3::zeros(),
7032 resolved_atmo,
7033 );
7034
7035 assert!(
7036 (segmented_acceleration.x - expected_acceleration.x).abs() < 1e-12,
7037 "Mach 1.5 must interpolate BC 0.3: segmented ax={} expected ax={}",
7038 segmented_acceleration.x,
7039 expected_acceleration.x
7040 );
7041 }
7042}
7043
7044#[cfg(test)]
7045mod custom_drag_table_validation_tests {
7046 use super::*;
7047
7048 #[test]
7049 fn solve_accepts_zero_bc_when_custom_table_present() {
7050 let inputs = BallisticInputs {
7051 bc_value: 0.0, bullet_mass: 0.0106,
7053 bullet_diameter: 0.00782,
7054 muzzle_velocity: 850.0,
7055 custom_drag_table: Some(crate::drag::DragTable::new(
7056 vec![0.5, 1.0, 2.0, 3.0],
7057 vec![0.23, 0.40, 0.30, 0.26],
7058 )),
7059 ..BallisticInputs::default()
7060 };
7061 let solver = TrajectorySolver::new(
7062 inputs,
7063 WindConditions::default(),
7064 AtmosphericConditions::default(),
7065 );
7066 assert!(solver.solve().is_ok());
7068 }
7069
7070 #[test]
7071 fn solve_still_requires_bc_without_table() {
7072 let inputs = BallisticInputs {
7073 bc_value: 0.0,
7074 bullet_mass: 0.0106,
7075 bullet_diameter: 0.00782,
7076 muzzle_velocity: 850.0,
7077 ..BallisticInputs::default()
7078 };
7079 let solver = TrajectorySolver::new(
7080 inputs,
7081 WindConditions::default(),
7082 AtmosphericConditions::default(),
7083 );
7084 assert!(solver.solve().is_err());
7085 }
7086}
7087
7088#[cfg(test)]
7090mod cd_scale_tests {
7091 use super::*;
7092
7093 fn deck() -> crate::drag::DragTable {
7094 crate::drag::DragTable::new(vec![0.5, 1.0, 2.0, 3.0], vec![0.23, 0.40, 0.30, 0.26])
7095 }
7096
7097 fn deck_inputs(cd_scale: f64) -> BallisticInputs {
7098 BallisticInputs {
7099 bullet_mass: 0.0106,
7100 bullet_diameter: 0.00782,
7101 muzzle_velocity: 850.0,
7102 custom_drag_table: Some(deck()),
7103 cd_scale,
7104 ..BallisticInputs::default()
7105 }
7106 }
7107
7108 #[test]
7109 fn default_cd_scale_is_one() {
7110 assert_eq!(BallisticInputs::default().cd_scale, 1.0);
7111 }
7112
7113 #[test]
7116 fn cd_scale_absent_is_byte_identical_to_explicit_one() {
7117 let omitted = BallisticInputs {
7118 bullet_mass: 0.0106,
7119 bullet_diameter: 0.00782,
7120 muzzle_velocity: 850.0,
7121 custom_drag_table: Some(deck()),
7122 ..BallisticInputs::default()
7123 };
7124 let explicit = BallisticInputs {
7125 cd_scale: 1.0,
7126 ..omitted.clone()
7127 };
7128
7129 let solver_omitted = TrajectorySolver::new(
7130 omitted,
7131 WindConditions::default(),
7132 AtmosphericConditions::default(),
7133 );
7134 let solver_explicit = TrajectorySolver::new(
7135 explicit,
7136 WindConditions::default(),
7137 AtmosphericConditions::default(),
7138 );
7139
7140 let cd_omitted = solver_omitted.calculate_drag_coefficient(700.0, 340.0);
7141 let cd_explicit = solver_explicit.calculate_drag_coefficient(700.0, 340.0);
7142 assert_eq!(
7143 cd_omitted.to_bits(),
7144 cd_explicit.to_bits(),
7145 "default cd_scale must be bit-identical to an explicit 1.0"
7146 );
7147
7148 let result = solver_omitted.solve();
7151 assert!(
7152 result.is_ok(),
7153 "existing custom-deck solves must pass unchanged"
7154 );
7155 }
7156
7157 #[test]
7159 fn cd_scale_multiplies_the_interpolated_cd_exactly() {
7160 let velocity = 700.0;
7161 let speed_of_sound = 340.0;
7162 let mach = velocity / speed_of_sound;
7163 let expected_unscaled = deck().interpolate(mach);
7164
7165 for &scale in &[0.90, 1.0, 1.10, 1.5] {
7166 let solver = TrajectorySolver::new(
7167 deck_inputs(scale),
7168 WindConditions::default(),
7169 AtmosphericConditions::default(),
7170 );
7171 let cd = solver.calculate_drag_coefficient(velocity, speed_of_sound);
7172 assert!(
7173 (cd - expected_unscaled * scale).abs() < 1e-12,
7174 "scale={scale}: cd={cd} expected={}",
7175 expected_unscaled * scale
7176 );
7177 }
7178 }
7179
7180 #[test]
7184 fn cd_scale_direction_on_cli_api_solver() {
7185 let solve = |scale: f64| {
7186 TrajectorySolver::new(
7187 deck_inputs(scale),
7188 WindConditions::default(),
7189 AtmosphericConditions::default(),
7190 )
7191 .solve()
7192 .expect("custom-deck solve should succeed")
7193 };
7194
7195 let baseline = solve(1.0);
7196 let scaled_up = solve(1.10);
7197 let scaled_down = solve(0.90);
7198
7199 assert!(
7200 scaled_up.impact_velocity < baseline.impact_velocity,
7201 "cd_scale=1.10 must increase drag -> lower impact velocity: base={} up={}",
7202 baseline.impact_velocity,
7203 scaled_up.impact_velocity
7204 );
7205 assert!(
7206 scaled_down.impact_velocity > baseline.impact_velocity,
7207 "cd_scale=0.90 must decrease drag -> higher impact velocity: base={} down={}",
7208 baseline.impact_velocity,
7209 scaled_down.impact_velocity
7210 );
7211 }
7212
7213 #[test]
7215 fn validate_for_solve_rejects_invalid_cd_scale() {
7216 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
7217 let solver = TrajectorySolver::new(
7218 deck_inputs(bad),
7219 WindConditions::default(),
7220 AtmosphericConditions::default(),
7221 );
7222 assert!(
7223 solver.solve().is_err(),
7224 "cd_scale={bad} must be rejected by validate_for_solve"
7225 );
7226 }
7227 }
7228
7229 #[test]
7236 fn validate_for_solve_rejects_invalid_cd_scale_without_a_custom_drag_table() {
7237 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
7238 let inputs = BallisticInputs {
7239 bc_value: 0.5,
7240 bc_type: crate::DragModel::G1,
7241 bullet_mass: 0.0106,
7242 bullet_diameter: 0.00782,
7243 muzzle_velocity: 850.0,
7244 cd_scale: bad,
7245 ..BallisticInputs::default()
7246 };
7247 assert!(
7248 inputs.custom_drag_table.is_none(),
7249 "precondition: no custom deck"
7250 );
7251 let solver = TrajectorySolver::new(
7252 inputs,
7253 WindConditions::default(),
7254 AtmosphericConditions::default(),
7255 );
7256 assert!(
7257 solver.solve().is_err(),
7258 "cd_scale={bad} must be rejected by validate_for_solve even without a custom \
7259 drag table"
7260 );
7261 }
7262 }
7263
7264 #[test]
7268 fn cd_scale_is_inert_without_a_custom_drag_table() {
7269 let make = |cd_scale: f64| BallisticInputs {
7270 bc_value: 0.5,
7271 bc_type: crate::DragModel::G1,
7272 bullet_mass: 0.0106,
7273 bullet_diameter: 0.00782,
7274 muzzle_velocity: 850.0,
7275 cd_scale,
7276 ..BallisticInputs::default()
7277 };
7278 let solver_neutral = TrajectorySolver::new(
7279 make(1.0),
7280 WindConditions::default(),
7281 AtmosphericConditions::default(),
7282 );
7283 let solver_far = TrajectorySolver::new(
7284 make(1.5),
7285 WindConditions::default(),
7286 AtmosphericConditions::default(),
7287 );
7288 let cd_neutral = solver_neutral.calculate_drag_coefficient(700.0, 340.0);
7289 let cd_far = solver_far.calculate_drag_coefficient(700.0, 340.0);
7290 assert_eq!(
7291 cd_neutral.to_bits(),
7292 cd_far.to_bits(),
7293 "cd_scale must not affect the G-model/BC drag path"
7294 );
7295 }
7296
7297 #[test]
7300 fn cd_scale_shifts_all_three_solver_paths_in_the_same_direction() {
7301 let cli_solve = |scale: f64| {
7303 TrajectorySolver::new(
7304 deck_inputs(scale),
7305 WindConditions::default(),
7306 AtmosphericConditions::default(),
7307 )
7308 .solve()
7309 .expect("cli_api custom-deck solve should succeed")
7310 };
7311 let cli_baseline = cli_solve(1.0);
7312 let cli_scaled = cli_solve(1.10);
7313 assert!(
7314 cli_scaled.impact_velocity < cli_baseline.impact_velocity,
7315 "cli_api: cd_scale=1.10 must lower impact velocity"
7316 );
7317
7318 let derivatives_accel_x = |scale: f64| {
7320 let inputs = deck_inputs(scale);
7321 crate::derivatives::compute_derivatives(
7322 nalgebra::Vector3::zeros(),
7323 nalgebra::Vector3::new(700.0, 0.0, 0.0),
7324 &inputs,
7325 nalgebra::Vector3::zeros(),
7326 (1.225, 340.0, 0.0, 0.0),
7327 inputs.bc_value,
7328 None,
7329 0.0,
7330 None,
7331 )[3]
7332 };
7333 let deriv_baseline = derivatives_accel_x(1.0);
7334 let deriv_scaled = derivatives_accel_x(1.10);
7335 assert!(
7336 deriv_scaled < deriv_baseline,
7337 "derivatives: cd_scale=1.10 must make x-acceleration more negative (more drag): \
7338 base={deriv_baseline} scaled={deriv_scaled}"
7339 );
7340
7341 let fast_final_speed = |scale: f64| {
7343 let inputs = deck_inputs(scale);
7344 let wind_sock = crate::wind::WindSock::new(vec![]);
7345 let params = crate::fast_trajectory::FastIntegrationParams {
7346 horiz: 500.0,
7347 vert: 0.0,
7348 initial_state: [0.0, 0.0, 0.0, 850.0, 0.0, 0.0],
7349 t_span: (0.0, 5.0),
7350 atmo_params: (0.0, 15.0, 1013.25, 1.0),
7351 atmo_sock: None,
7352 };
7353 let solution = crate::fast_trajectory::fast_integrate(&inputs, &wind_sock, params);
7354 assert!(
7355 solution.success,
7356 "fast_integrate must succeed for scale={scale}"
7357 );
7358 let last = solution.t.len() - 1;
7359 let (vx, vy, vz) = (
7360 solution.y[3][last],
7361 solution.y[4][last],
7362 solution.y[5][last],
7363 );
7364 (vx * vx + vy * vy + vz * vz).sqrt()
7365 };
7366 let fast_baseline = fast_final_speed(1.0);
7367 let fast_scaled = fast_final_speed(1.10);
7368 assert!(
7369 fast_scaled < fast_baseline,
7370 "fast_trajectory: cd_scale=1.10 must lower final speed: base={fast_baseline} scaled={fast_scaled}"
7371 );
7372 }
7373}
7374
7375#[cfg(test)]
7376mod humid_local_mach_tests {
7377 use super::*;
7378
7379 fn solver_with_station_humidity(humidity_percent: f64) -> TrajectorySolver {
7380 let inputs = BallisticInputs {
7381 custom_drag_table: Some(crate::drag::DragTable::new(vec![0.5, 1.5], vec![0.1, 1.1])),
7382 ..BallisticInputs::default()
7383 };
7384 TrajectorySolver::new(
7385 inputs,
7386 WindConditions::default(),
7387 AtmosphericConditions {
7388 temperature: 30.0,
7389 pressure: 1013.25,
7390 humidity: humidity_percent,
7391 altitude: 0.0,
7392 },
7393 )
7394 }
7395
7396 fn acceleration(solver: &TrajectorySolver, base_ratio: f64) -> Vector3<f64> {
7397 solver.calculate_acceleration(
7398 &Vector3::zeros(),
7399 &Vector3::new(350.0, 0.0, 0.0),
7400 &Vector3::zeros(),
7401 (30.0, 1013.25, base_ratio),
7402 )
7403 }
7404
7405 #[test]
7406 fn local_mach_uses_station_humidity_when_density_is_held_constant() {
7407 let dry = acceleration(&solver_with_station_humidity(0.0), 1.0);
7408 let humid = acceleration(&solver_with_station_humidity(100.0), 1.0);
7409
7410 assert!(
7411 humid.x > dry.x,
7412 "humid sound speed should lower Mach and drag on the rising test curve: dry ax={} humid ax={}",
7413 dry.x,
7414 humid.x
7415 );
7416 }
7417
7418 #[test]
7419 fn active_atmosphere_zone_uses_zone_humidity_instead_of_station_humidity() {
7420 let zone_humidity = 80.0;
7421 let zone_ratio =
7422 crate::atmosphere::calculate_air_density_cimp(30.0, 1013.25, zone_humidity) / 1.225;
7423 let station_solver = solver_with_station_humidity(zone_humidity);
7424 let mut zoned_solver = solver_with_station_humidity(0.0);
7425 zoned_solver.set_atmo_segments(vec![(30.0, 1013.25, zone_humidity, 1_000.0)]);
7426
7427 let station = acceleration(&station_solver, zone_ratio);
7428 let zoned = acceleration(&zoned_solver, zone_ratio);
7429
7430 assert!(
7431 (zoned - station).norm() < 1e-12,
7432 "active zone T/P/RH should override the station atmosphere: station={station:?} zoned={zoned:?}"
7433 );
7434 }
7435}
7436
7437#[cfg(test)]
7438mod inclined_atmosphere_frame_tests {
7439 use super::*;
7440
7441 fn expected_shot_frame_vector(level: Vector3<f64>, angle: f64) -> Vector3<f64> {
7442 let (sin_angle, cos_angle) = angle.sin_cos();
7443 Vector3::new(
7444 level.x * cos_angle + level.y * sin_angle,
7445 -level.x * sin_angle + level.y * cos_angle,
7446 level.z,
7447 )
7448 }
7449
7450 #[test]
7451 fn inclined_positions_at_same_world_altitude_have_same_solver_acceleration() {
7452 let angle = std::f64::consts::FRAC_PI_6;
7453 let inputs = BallisticInputs {
7454 shooting_angle: angle,
7455 ..BallisticInputs::default()
7456 };
7457 let atmosphere = AtmosphericConditions {
7458 altitude: 100.0,
7459 ..AtmosphericConditions::default()
7460 };
7461 let solver = TrajectorySolver::new(inputs, WindConditions::default(), atmosphere);
7462 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7463 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
7464 let velocity = Vector3::new(600.0, 0.0, 0.0);
7465 let along_slant = Vector3::new(1_000.0, 0.0, 0.0);
7466 let across_slant = Vector3::new(0.0, 500.0 / angle.cos(), 0.0);
7467
7468 let a = solver.calculate_acceleration(
7469 &along_slant,
7470 &velocity,
7471 &Vector3::zeros(),
7472 resolved_atmo,
7473 );
7474 let b = solver.calculate_acceleration(
7475 &across_slant,
7476 &velocity,
7477 &Vector3::zeros(),
7478 resolved_atmo,
7479 );
7480
7481 assert!(
7482 (a - b).norm() < 1e-10,
7483 "solver acceleration differs at equal world altitude: {a:?} vs {b:?}"
7484 );
7485 }
7486
7487 #[test]
7488 fn inclined_headwind_is_rotated_into_solver_frame() {
7489 let angle = std::f64::consts::FRAC_PI_6;
7490 let inputs = BallisticInputs {
7491 shooting_angle: angle,
7492 ..BallisticInputs::default()
7493 };
7494 let solver = TrajectorySolver::new(
7495 inputs,
7496 WindConditions::default(),
7497 AtmosphericConditions::default(),
7498 );
7499 let level_headwind = Vector3::new(-100.0, 0.0, 0.0);
7500 let velocity = expected_shot_frame_vector(level_headwind, angle);
7501 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7502 let actual = solver.calculate_acceleration(
7503 &Vector3::zeros(),
7504 &velocity,
7505 &level_headwind,
7506 (temp_c, pressure_hpa, density / 1.225),
7507 );
7508
7509 assert!(
7510 (actual - solver.gravity_acceleration()).norm() < 1e-12,
7511 "co-moving horizontal wind must leave only shot-frame gravity: {actual:?}"
7512 );
7513 }
7514
7515 #[test]
7516 fn inclined_coriolis_is_rotated_into_solver_frame() {
7517 let angle = std::f64::consts::FRAC_PI_6;
7518 let latitude_deg = 45.0_f64;
7519 let shot_azimuth = 0.4_f64;
7520 let velocity = Vector3::new(600.0, 20.0, 5.0);
7521 let base_inputs = BallisticInputs {
7522 shooting_angle: angle,
7523 latitude: Some(latitude_deg),
7524 shot_azimuth,
7525 ..BallisticInputs::default()
7526 };
7527 let acceleration = |enable_coriolis| {
7528 let mut inputs = base_inputs.clone();
7529 inputs.enable_coriolis = enable_coriolis;
7530 let solver = TrajectorySolver::new(
7531 inputs,
7532 WindConditions::default(),
7533 AtmosphericConditions::default(),
7534 );
7535 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7536 solver.calculate_acceleration(
7537 &Vector3::zeros(),
7538 &velocity,
7539 &Vector3::zeros(),
7540 (temp_c, pressure_hpa, density / 1.225),
7541 )
7542 };
7543
7544 let omega_earth = 7.2921159e-5_f64;
7545 let latitude = latitude_deg.to_radians();
7546 let level_omega = Vector3::new(
7547 omega_earth * latitude.cos() * shot_azimuth.cos(),
7548 omega_earth * latitude.sin(),
7549 -omega_earth * latitude.cos() * shot_azimuth.sin(),
7550 );
7551 let expected = -2.0 * expected_shot_frame_vector(level_omega, angle).cross(&velocity);
7552 let actual = acceleration(true) - acceleration(false);
7553
7554 assert!(
7555 (actual - expected).norm() < 1e-12,
7556 "inclined Coriolis mismatch: actual={actual:?}, expected={expected:?}"
7557 );
7558 }
7559}
7560
7561#[cfg(test)]
7562mod terminal_range_interpolation_tests {
7563 use super::*;
7564
7565 #[test]
7566 fn terminal_finalizer_selects_the_earliest_crossed_boundary() {
7567 let inputs = BallisticInputs {
7568 ground_threshold: 0.0,
7569 ..BallisticInputs::default()
7570 };
7571 let mut solver = TrajectorySolver::new(
7572 inputs,
7573 WindConditions::default(),
7574 AtmosphericConditions::default(),
7575 );
7576 solver.set_max_range(120.0);
7577
7578 let previous_speed = 700.0;
7579 let mut points = vec![TrajectoryPoint {
7580 time: 99.0,
7581 position: Vector3::new(90.0, 1.0, -1.0),
7582 velocity_magnitude: previous_speed,
7583 kinetic_energy: 0.5 * solver.inputs.bullet_mass * previous_speed.powi(2),
7584 drag_coefficient: None,
7585 }];
7586 let mut max_height = 1.0;
7587 let termination = solver
7588 .append_terminal_endpoint(
7589 &mut points,
7590 Vector3::new(130.0, -3.0, 3.0),
7591 Vector3::new(600.0, 0.0, 0.0),
7592 101.0,
7593 &mut max_height,
7594 &mut 0,
7595 )
7596 .expect("the final step brackets supported boundaries");
7597
7598 assert_eq!(termination, TrajectoryTermination::GroundThreshold);
7599 assert_eq!(points.len(), 2);
7600 let terminal = points.last().expect("terminal point");
7601 assert_eq!(terminal.time, 99.5);
7602 assert_eq!(terminal.position, Vector3::new(100.0, 0.0, 0.0));
7603 assert_eq!(terminal.velocity_magnitude, 675.0);
7604 assert_eq!(
7605 terminal.kinetic_energy,
7606 0.5 * solver.inputs.bullet_mass * 675.0_f64.powi(2)
7607 );
7608
7609 solver.set_max_range(100.0);
7611 let mut tied_points = vec![points[0].clone()];
7612 assert_eq!(
7613 solver
7614 .append_terminal_endpoint(
7615 &mut tied_points,
7616 Vector3::new(130.0, -3.0, 3.0),
7617 Vector3::new(600.0, 0.0, 0.0),
7618 101.0,
7619 &mut max_height,
7620 &mut 0,
7621 )
7622 .expect("tied boundaries remain a valid terminal"),
7623 TrajectoryTermination::GroundThreshold
7624 );
7625 }
7626
7627 #[test]
7628 fn sub_ulp_terminal_crossing_replaces_instead_of_duplicating_range() {
7629 let ground_threshold = f64::from_bits(1.0_f64.to_bits() - 1);
7630 let inputs = BallisticInputs {
7631 ground_threshold,
7632 ..BallisticInputs::default()
7633 };
7634 let mut solver = TrajectorySolver::new(
7635 inputs,
7636 WindConditions::default(),
7637 AtmosphericConditions::default(),
7638 );
7639 solver.set_max_range(1_000.0);
7640
7641 let speed = 700.0;
7642 let mut points = vec![TrajectoryPoint {
7643 time: 0.0,
7644 position: Vector3::new(100.0, 1.0, 0.0),
7645 velocity_magnitude: speed,
7646 kinetic_energy: 0.5 * solver.inputs.bullet_mass * speed.powi(2),
7647 drag_coefficient: None,
7648 }];
7649 let mut max_height = 1.0;
7650 let termination = solver
7651 .append_terminal_endpoint(
7652 &mut points,
7653 Vector3::new(101.0, 0.0, 0.0),
7654 Vector3::new(699.0, 0.0, 0.0),
7655 1.0,
7656 &mut max_height,
7657 &mut 0,
7658 )
7659 .expect("sub-ULP ground crossing remains representable as one terminal state");
7660
7661 assert_eq!(termination, TrajectoryTermination::GroundThreshold);
7662 assert_eq!(points.len(), 1);
7663 assert_eq!(points[0].position.x, 100.0);
7664 assert_eq!(points[0].position.y.to_bits(), ground_threshold.to_bits());
7665 assert!(points[0].time > 0.0);
7666 }
7667
7668 #[test]
7674 fn downrange_stall_replaces_instead_of_duplicating_range() {
7675 for (name, use_rk4, use_adaptive_rk45, bc_value) in [
7676 ("Euler", false, false, 2.0e-4),
7677 ("RK4", true, false, 1.0e-4),
7678 ("RK45", true, true, 1.0e-4),
7679 ] {
7680 let inputs = BallisticInputs {
7681 bc_value,
7682 bc_type: DragModel::G7,
7683 bullet_mass: 0.01134,
7684 bullet_diameter: 0.00782,
7685 muzzle_velocity: 823.0,
7686 use_rk4,
7687 use_adaptive_rk45,
7688 ground_threshold: -100.0,
7689 ..BallisticInputs::default()
7690 };
7691 let mut solver = TrajectorySolver::new(
7692 inputs,
7693 WindConditions::default(),
7694 AtmosphericConditions::default(),
7695 );
7696 solver.set_max_range(300.0);
7697 solver.set_time_step(0.001);
7698
7699 let result = solver
7700 .solve()
7701 .unwrap_or_else(|error| panic!("{name}: a stalled flight still solves: {error}"));
7702 assert!(
7703 result
7704 .points
7705 .windows(2)
7706 .all(|pair| pair[1].position.x > pair[0].position.x),
7707 "{name}: stored downrange distances must strictly increase"
7708 );
7709 assert_eq!(
7710 result.termination,
7711 TrajectoryTermination::GroundThreshold,
7712 "{name}"
7713 );
7714 assert!(
7715 result.max_range > 0.0 && result.max_range < 300.0,
7716 "{name}: the stall ends short of max_range, got {} m",
7717 result.max_range
7718 );
7719 let terminal = result.points.last().expect("terminal point");
7720 assert_eq!(terminal.position.y, -100.0, "{name}");
7721 result
7722 .sample_observations(50.0, 100)
7723 .unwrap_or_else(|error| panic!("{name}: checked sampling refused: {error}"));
7724 }
7725 }
7726
7727 #[test]
7728 fn every_solver_appends_an_exact_max_range_endpoint() {
7729 let target_range = 0.1;
7730 let modes = [
7731 ("Euler", false, false),
7732 ("RK4", true, false),
7733 ("RK45", true, true),
7734 ];
7735
7736 for (name, use_rk4, use_adaptive_rk45) in modes {
7737 let inputs = BallisticInputs {
7738 use_rk4,
7739 use_adaptive_rk45,
7740 ground_threshold: f64::NEG_INFINITY,
7741 enable_trajectory_sampling: true,
7742 sample_interval: target_range,
7743 ..BallisticInputs::default()
7744 };
7745 let mut solver = TrajectorySolver::new(
7746 inputs,
7747 WindConditions::default(),
7748 AtmosphericConditions::default(),
7749 );
7750 solver.set_max_range(target_range);
7751
7752 let result = solver.solve().expect("short-range solve should succeed");
7753 let terminal = result.points.last().expect("terminal point is missing");
7754 let muzzle = result.points.first().expect("muzzle point is missing");
7755
7756 assert_eq!(result.termination, TrajectoryTermination::MaxRange);
7757 assert_eq!(
7758 terminal.position.x.to_bits(),
7759 target_range.to_bits(),
7760 "{name} did not terminate exactly at max_range"
7761 );
7762 assert_eq!(result.max_range.to_bits(), target_range.to_bits());
7763 assert!(
7764 result.time_of_flight > 0.0 && result.time_of_flight < solver.time_step,
7765 "{name} terminal time was not interpolated within the crossing step: {}",
7766 result.time_of_flight
7767 );
7768 assert_eq!(result.time_of_flight.to_bits(), terminal.time.to_bits());
7769 assert_eq!(
7770 result.impact_velocity.to_bits(),
7771 terminal.velocity_magnitude.to_bits()
7772 );
7773 assert_eq!(
7774 result.impact_energy.to_bits(),
7775 terminal.kinetic_energy.to_bits()
7776 );
7777 let expected_energy = 0.5 * solver.inputs.bullet_mass * result.impact_velocity.powi(2);
7778 assert!((result.impact_energy - expected_energy).abs() < 1e-12);
7779 assert!(terminal.velocity_magnitude < muzzle.velocity_magnitude);
7780 assert!(terminal.kinetic_energy < muzzle.kinetic_energy);
7781
7782 let terminal_sample = result
7783 .sampled_points
7784 .as_ref()
7785 .and_then(|samples| samples.last())
7786 .expect("terminal trajectory sample is missing");
7787 assert_eq!(
7788 terminal_sample.distance_m.to_bits(),
7789 target_range.to_bits(),
7790 "{name} sampling did not include max_range"
7791 );
7792 assert_eq!(
7793 terminal_sample.time_s.to_bits(),
7794 result.time_of_flight.to_bits()
7795 );
7796 assert_eq!(
7797 terminal_sample.velocity_mps.to_bits(),
7798 result.impact_velocity.to_bits()
7799 );
7800 assert!((terminal_sample.energy_j - result.impact_energy).abs() < 1e-12);
7801 }
7802 }
7803}
7804
7805#[cfg(test)]
7806mod precession_inertia_wiring_tests {
7807 use super::*;
7808
7809 #[test]
7810 fn solver_uses_projectile_specific_moments_of_inertia() {
7811 let mass_kg = 55.0 * crate::constants::GRAINS_TO_KG;
7812 let caliber_m = 0.224 * 0.0254;
7813 let length_m = 0.75 * 0.0254;
7814 let inputs = BallisticInputs {
7815 bullet_mass: mass_kg,
7816 bullet_diameter: caliber_m,
7817 bullet_length: length_m,
7818 muzzle_velocity: 800.0,
7819 twist_rate: 7.0,
7820 enable_precession_nutation: true,
7821 use_rk4: false,
7822 use_adaptive_rk45: false,
7823 ..BallisticInputs::default()
7824 };
7825 let mut solver = TrajectorySolver::new(
7826 inputs,
7827 WindConditions::default(),
7828 AtmosphericConditions::default(),
7829 );
7830 solver.set_max_range(0.1);
7831
7832 let (air_density, speed_of_sound, _, _) = solver.resolved_atmosphere();
7833 let velocity_mps = solver.inputs.muzzle_velocity;
7834 let velocity_fps = velocity_mps * 3.28084;
7835 let twist_rate_ft = solver.inputs.twist_rate / 12.0;
7836 let spin_rate_rad_s = (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI;
7837 let initial_state = AngularState {
7838 pitch_angle: 0.001,
7839 yaw_angle: 0.001,
7840 pitch_rate: 0.0,
7841 yaw_rate: 0.0,
7842 precession_angle: 0.0,
7843 nutation_phase: 0.0,
7844 };
7845 let params = PrecessionNutationParams {
7846 mass_kg,
7847 caliber_m,
7848 length_m,
7849 spin_rate_rad_s,
7850 spin_inertia: crate::spin_decay::calculate_moment_of_inertia(
7851 mass_kg, caliber_m, length_m, "ogive",
7852 ),
7853 transverse_inertia: crate::pitch_damping::calculate_transverse_moment_of_inertia(
7854 mass_kg, caliber_m, length_m, "ogive",
7855 ),
7856 velocity_mps,
7857 air_density_kg_m3: air_density,
7858 mach: velocity_mps / speed_of_sound,
7859 pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
7860 nutation_damping_factor: 0.05,
7861 };
7862 let expected = calculate_combined_angular_motion(
7863 ¶ms,
7864 &initial_state,
7865 0.0,
7866 solver.time_step,
7867 0.001,
7868 );
7869 let actual = solver
7870 .solve()
7871 .expect("one-step solve should succeed")
7872 .angular_state
7873 .expect("precession/nutation was enabled");
7874
7875 assert!(
7876 (actual.precession_angle - expected.precession_angle).abs() < 1e-15,
7877 "precession phase used the wrong inertia: actual={}, expected={}",
7878 actual.precession_angle,
7879 expected.precession_angle
7880 );
7881 assert!(
7882 (actual.nutation_phase - expected.nutation_phase).abs() < 1e-15,
7883 "nutation phase used the wrong inertia: actual={}, expected={}",
7884 actual.nutation_phase,
7885 expected.nutation_phase
7886 );
7887 }
7888}
7889
7890#[cfg(test)]
7891mod form_factor_drag_tests {
7892 use super::*;
7893
7894 fn acceleration_with_form_factor_flag(enabled: bool) -> Vector3<f64> {
7895 let inputs = BallisticInputs {
7896 bc_value: 0.462,
7897 bc_type: DragModel::G1,
7898 bullet_model: Some("168gr SMK Match".to_string()),
7899 use_form_factor: enabled,
7900 ..BallisticInputs::default()
7901 };
7902 let solver = TrajectorySolver::new(
7903 inputs,
7904 WindConditions::default(),
7905 AtmosphericConditions::default(),
7906 );
7907 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7908 solver.calculate_acceleration(
7909 &Vector3::zeros(),
7910 &Vector3::new(600.0, 0.0, 0.0),
7911 &Vector3::zeros(),
7912 (temp_c, pressure_hpa, density / 1.225),
7913 )
7914 }
7915
7916 #[test]
7917 fn measured_bc_drag_does_not_apply_name_based_form_factor_again() {
7918 let baseline = acceleration_with_form_factor_flag(false);
7919 let flagged = acceleration_with_form_factor_flag(true);
7920
7921 assert!(
7922 (flagged - baseline).norm() < 1e-12,
7923 "published BC already encodes form factor: baseline={baseline:?} flagged={flagged:?}"
7924 );
7925 }
7926}
7927
7928#[cfg(test)]
7929mod rk45_adaptivity_tests {
7930 use super::*;
7931
7932 #[test]
7933 fn cli_rk45_error_norm_scales_components_independently() {
7934 let position = Vector3::new(1.0e9, 0.0, 0.0);
7935 let velocity = Vector3::new(800.0, 0.0, 0.0);
7936 let fifth_position = position;
7937 let fifth_velocity = velocity;
7938 let fourth_position = position;
7939 let fourth_velocity = Vector3::new(800.0, 1.0e-3, 0.0);
7940
7941 let error = cli_rk45_error_norm(
7942 &position,
7943 &velocity,
7944 &fifth_position,
7945 &fifth_velocity,
7946 &fourth_position,
7947 &fourth_velocity,
7948 );
7949 let expected = 1.0e-3 / 6.0_f64.sqrt();
7950
7951 assert!(
7952 (error - expected).abs() <= 1e-15,
7953 "large downrange position masked a velocity-component error: {error}"
7954 );
7955 }
7956
7957 fn discontinuous_wind_solver() -> TrajectorySolver {
7958 let inputs = BallisticInputs::default();
7959 let mut solver = TrajectorySolver::new(
7960 inputs,
7961 WindConditions::default(),
7962 AtmosphericConditions::default(),
7963 );
7964 solver.set_wind_segments(vec![
7965 crate::wind::WindSegment::new(0.0, 90.0, 4.0),
7966 crate::wind::WindSegment::new(1_000.0, 90.0, 10_000.0),
7967 ]);
7968 solver
7969 }
7970
7971 #[test]
7972 fn rk45_retries_discontinuous_trial_before_advancing() {
7973 let solver = discontinuous_wind_solver();
7974 let position = Vector3::new(0.0, solver.inputs.muzzle_height, 0.0);
7975 let velocity = Vector3::new(solver.inputs.muzzle_velocity, 0.0, 0.0);
7976 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7977 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
7978 let dt = 0.01;
7979
7980 let rejected_trial = solver.rk45_step(
7981 &position,
7982 &velocity,
7983 dt,
7984 &Vector3::zeros(),
7985 RK45_TOLERANCE,
7986 resolved_atmo,
7987 );
7988 assert!(
7989 rejected_trial.error > RK45_TOLERANCE,
7990 "discontinuous full step must exceed tolerance, got {}",
7991 rejected_trial.error
7992 );
7993
7994 let accepted =
7995 solver.adaptive_rk45_step(&position, &velocity, dt, &Vector3::zeros(), resolved_atmo);
7996 assert!(accepted.used_dt < dt, "oversized trial was not retried");
7997 assert!(
7998 accepted.error <= RK45_TOLERANCE || accepted.used_dt <= RK45_MIN_DT,
7999 "accepted error {} exceeds tolerance at dt {}",
8000 accepted.error,
8001 accepted.used_dt
8002 );
8003
8004 let accepted_trial = solver.rk45_step(
8005 &position,
8006 &velocity,
8007 accepted.used_dt,
8008 &Vector3::zeros(),
8009 RK45_TOLERANCE,
8010 resolved_atmo,
8011 );
8012 assert_eq!(accepted.position, accepted_trial.position);
8013 assert_eq!(accepted.velocity, accepted_trial.velocity);
8014 assert!((RK45_MIN_DT..=RK45_MAX_DT).contains(&accepted.next_dt));
8015 }
8016}
8017
8018#[cfg(test)]
8019mod ground_termination_tests {
8020 use super::*;
8021 use crate::trajectory_observation::TrajectoryObservationFlag;
8022
8023 #[test]
8024 fn every_solver_reports_one_exact_early_ground_endpoint() {
8025 for (name, use_rk4, use_adaptive_rk45) in [
8026 ("Euler", false, false),
8027 ("RK4", true, false),
8028 ("RK45", true, true),
8029 ] {
8030 let inputs = BallisticInputs {
8031 muzzle_height: 1.0,
8032 muzzle_angle: -0.2,
8033 ground_threshold: 0.0,
8034 use_rk4,
8035 use_adaptive_rk45,
8036 ..BallisticInputs::default()
8037 };
8038 let mut solver = TrajectorySolver::new(
8039 inputs,
8040 WindConditions::default(),
8041 AtmosphericConditions::default(),
8042 );
8043 solver.set_max_range(1_000.0);
8044
8045 let result = solver.solve().expect("early-ground solve should succeed");
8046 let terminal = result.points.last().expect("terminal point is missing");
8047
8048 assert_eq!(result.termination, TrajectoryTermination::GroundThreshold);
8049 assert_eq!(terminal.position.y.to_bits(), 0.0_f64.to_bits());
8050 assert!(
8051 terminal.position.x < 1_000.0,
8052 "{name} incorrectly reached max range"
8053 );
8054 assert_eq!(result.max_range.to_bits(), terminal.position.x.to_bits());
8055 assert_eq!(
8056 result
8057 .points
8058 .iter()
8059 .filter(|point| point.position.y == 0.0)
8060 .count(),
8061 1,
8062 "{name} did not retain exactly one ground endpoint"
8063 );
8064
8065 let observations = result
8066 .sample_observations(1.0, 100)
8067 .expect("checked early-ground sampling should succeed");
8068 assert!(observations[..observations.len() - 1]
8069 .iter()
8070 .all(|observation| observation.distance_m < terminal.position.x));
8071 let terminal_observation = observations.last().expect("terminal observation");
8072 assert_eq!(
8073 terminal_observation.distance_m.to_bits(),
8074 terminal.position.x.to_bits()
8075 );
8076 assert!(terminal_observation
8077 .flags
8078 .contains(&TrajectoryObservationFlag::Terminal));
8079 assert!(terminal_observation
8080 .flags
8081 .contains(&TrajectoryObservationFlag::GroundThreshold));
8082 assert_eq!(
8083 observations
8084 .iter()
8085 .filter(|observation| observation
8086 .flags
8087 .contains(&TrajectoryObservationFlag::Terminal))
8088 .count(),
8089 1,
8090 "{name} repeated the terminal observation"
8091 );
8092 }
8093 }
8094
8095 #[test]
8100 fn rk4_and_rk45_descend_to_ground_threshold() {
8101 for adaptive in [false, true] {
8102 let inputs = BallisticInputs {
8103 muzzle_angle: 0.1, use_rk4: true,
8105 use_adaptive_rk45: adaptive,
8106 ..BallisticInputs::default()
8107 };
8108 assert_eq!(
8109 inputs.ground_threshold, -100.0,
8110 "default ground_threshold is -100 m"
8111 );
8112
8113 let mut solver = TrajectorySolver::new(
8114 inputs,
8115 WindConditions::default(),
8116 AtmosphericConditions::default(),
8117 );
8118 solver.set_max_range(1.0e7);
8120
8121 let result = solver.solve().expect("solve should succeed");
8122 let final_y = result
8123 .points
8124 .last()
8125 .expect("trajectory has points")
8126 .position
8127 .y;
8128 assert!(
8129 final_y < -1.0,
8130 "adaptive_rk45={adaptive}: final y = {final_y} m; a lofted shot should descend \
8131 past launch level toward the ground_threshold floor, not stop at y = 0"
8132 );
8133 }
8134 }
8135}
8136
8137#[cfg(test)]
8138mod magnus_stability_tests {
8139 use super::*;
8140
8141 #[test]
8142 fn yaw_of_repose_magnus_force_is_vertical_and_twist_signed() {
8143 let acceleration = |enable_magnus, is_twist_right| {
8144 let inputs = BallisticInputs {
8145 muzzle_velocity: 822.96,
8146 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
8147 bullet_diameter: 0.308 * 0.0254,
8148 bullet_length: 1.215 * 0.0254,
8149 twist_rate: 10.0,
8150 is_twist_right,
8151 enable_magnus,
8152 ..BallisticInputs::default()
8153 };
8154 let solver = TrajectorySolver::new(
8155 inputs,
8156 WindConditions::default(),
8157 AtmosphericConditions::default(),
8158 );
8159 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
8160 solver.calculate_acceleration(
8161 &Vector3::zeros(),
8162 &Vector3::new(822.96, 0.0, 0.0),
8163 &Vector3::zeros(),
8164 (temp_c, pressure_hpa, density / 1.225),
8165 )
8166 };
8167
8168 let baseline = acceleration(false, true);
8169 let right_twist = acceleration(true, true) - baseline;
8170 let left_twist = acceleration(true, false) - baseline;
8171
8172 assert!(
8173 right_twist.y < 0.0,
8174 "right-hand Magnus must point down, got {right_twist:?}"
8175 );
8176 assert!(
8177 left_twist.y > 0.0,
8178 "left-hand Magnus must point up, got {left_twist:?}"
8179 );
8180 assert!((right_twist.y + left_twist.y).abs() < 1e-12);
8181 assert!(right_twist.x.abs() < 1e-12 && right_twist.z.abs() < 1e-12);
8182 assert!(left_twist.x.abs() < 1e-12 && left_twist.z.abs() < 1e-12);
8183 }
8184
8185 #[test]
8186 fn magnus_uses_velocity_corrected_muzzle_stability_gate() {
8187 let muzzle_velocity = 1_400.0 / 3.28084;
8188 let inputs = BallisticInputs {
8189 muzzle_velocity,
8190 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
8191 bullet_diameter: 0.308 * 0.0254,
8192 bullet_length: 1.215 * 0.0254,
8193 twist_rate: 15.0,
8194 enable_magnus: true,
8195 ..BallisticInputs::default()
8196 };
8197 let solver = TrajectorySolver::new(
8198 inputs.clone(),
8199 WindConditions::default(),
8200 AtmosphericConditions::default(),
8201 );
8202
8203 let bare_sg = crate::spin_drift::miller_stability(0.308, 168.0, 15.0, 1.215);
8204 let canonical_sg = solver.effective_spin_drift_sg();
8205 assert!(bare_sg > 1.0, "test requires bare Sg above the Magnus gate");
8206 assert!(
8207 canonical_sg < 1.0,
8208 "velocity-corrected Sg must be below the gate, got {canonical_sg}"
8209 );
8210
8211 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
8212 let acceleration = solver.calculate_acceleration(
8213 &Vector3::zeros(),
8214 &Vector3::new(muzzle_velocity, 0.0, 0.0),
8215 &Vector3::zeros(),
8216 (temp_c, pressure_hpa, density / 1.225),
8217 );
8218 let mut baseline_inputs = inputs;
8219 baseline_inputs.enable_magnus = false;
8220 let baseline_solver = TrajectorySolver::new(
8221 baseline_inputs,
8222 WindConditions::default(),
8223 AtmosphericConditions::default(),
8224 );
8225 let baseline = baseline_solver.calculate_acceleration(
8226 &Vector3::zeros(),
8227 &Vector3::new(muzzle_velocity, 0.0, 0.0),
8228 &Vector3::zeros(),
8229 (temp_c, pressure_hpa, density / 1.225),
8230 );
8231
8232 assert_eq!(
8233 acceleration, baseline,
8234 "canonical Sg below 1 must suppress every Magnus acceleration component"
8235 );
8236 }
8237
8238 #[test]
8239 fn magnus_force_grows_as_fixed_spin_projectile_slows() {
8240 let inputs = BallisticInputs {
8241 muzzle_velocity: 800.0,
8242 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
8243 bullet_diameter: 0.308 * 0.0254,
8244 bullet_length: 1.215 * 0.0254,
8245 twist_rate: 12.0,
8246 enable_magnus: true,
8247 ..BallisticInputs::default()
8248 };
8249
8250 let magnus_acceleration = |speed_mps| {
8251 let evaluate = |enable_magnus| {
8252 let mut run_inputs = inputs.clone();
8253 run_inputs.enable_magnus = enable_magnus;
8254 let solver = TrajectorySolver::new(
8255 run_inputs,
8256 WindConditions::default(),
8257 AtmosphericConditions::default(),
8258 );
8259 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
8260 solver
8261 .calculate_acceleration(
8262 &Vector3::zeros(),
8263 &Vector3::new(speed_mps, 0.0, 0.0),
8264 &Vector3::zeros(),
8265 (temp_c, pressure_hpa, density / 1.225),
8266 )
8267 .y
8268 };
8269 (evaluate(true) - evaluate(false)).abs()
8270 };
8271
8272 let fast = magnus_acceleration(200.0);
8273 let slow = magnus_acceleration(100.0);
8274 let ratio = slow / fast;
8275 let expected_ratio = 2.0_f64.powf(5.0 / 3.0);
8276
8277 assert!(fast > 0.0 && slow > 0.0, "fast={fast}, slow={slow}");
8278 assert!(
8279 (ratio - expected_ratio).abs() < 1e-3,
8280 "fixed-spin Magnus acceleration must grow downrange; slow/fast={ratio}, \
8281 expected={expected_ratio}"
8282 );
8283 }
8284}
8285
8286#[cfg(test)]
8287mod coriolis_direction_tests {
8288 use super::*;
8289 use std::f64::consts::FRAC_PI_2;
8290
8291 #[test]
8292 fn supersonic_crossing_flags_a_positive_range_sample() {
8293 use crate::trajectory_sampling::TrajectoryFlag;
8297
8298 for (solver_name, use_rk4, use_adaptive_rk45) in [
8299 ("Euler", false, false),
8300 ("RK4", true, false),
8301 ("RK45", true, true),
8302 ] {
8303 let inputs = BallisticInputs {
8304 muzzle_velocity: 850.0,
8305 bc_value: 0.2,
8306 bc_type: DragModel::G7,
8307 muzzle_angle: 0.03,
8308 enable_trajectory_sampling: true,
8309 sample_interval: 50.0,
8310 use_rk4,
8311 use_adaptive_rk45,
8312 ..BallisticInputs::default()
8313 };
8314 let mut solver = TrajectorySolver::new(
8315 inputs,
8316 WindConditions::default(),
8317 AtmosphericConditions::default(),
8318 );
8319 solver.set_max_range(2000.0);
8320 let samples = solver
8321 .solve()
8322 .expect("supersonic solve should succeed")
8323 .sampled_points
8324 .expect("sampling was enabled");
8325 let flagged_distances: Vec<_> = samples
8326 .iter()
8327 .filter(|sample| sample.flags.contains(&TrajectoryFlag::MachTransition))
8328 .map(|sample| sample.distance_m)
8329 .collect();
8330
8331 assert!(
8332 !flagged_distances.is_empty()
8333 && flagged_distances.iter().all(|distance| *distance > 0.0),
8334 "{solver_name} must flag genuine crossings only at positive range: {flagged_distances:?}"
8335 );
8336 }
8337 }
8338
8339 #[test]
8340 fn subsonic_launch_does_not_flag_a_muzzle_transition() {
8341 use crate::trajectory_sampling::TrajectoryFlag;
8342
8343 for (solver_name, use_rk4, use_adaptive_rk45) in [
8344 ("Euler", false, false),
8345 ("RK4", true, false),
8346 ("RK45", true, true),
8347 ] {
8348 let inputs = BallisticInputs {
8349 muzzle_velocity: 250.0,
8350 muzzle_angle: 0.02,
8351 enable_trajectory_sampling: true,
8352 sample_interval: 25.0,
8353 use_rk4,
8354 use_adaptive_rk45,
8355 ..BallisticInputs::default()
8356 };
8357 let mut solver = TrajectorySolver::new(
8358 inputs,
8359 WindConditions::default(),
8360 AtmosphericConditions::default(),
8361 );
8362 solver.set_max_range(300.0);
8363 let samples = solver
8364 .solve()
8365 .expect("subsonic solve should succeed")
8366 .sampled_points
8367 .expect("sampling was enabled");
8368
8369 assert!(
8370 samples
8371 .iter()
8372 .all(|sample| !sample.flags.contains(&TrajectoryFlag::MachTransition)),
8373 "{solver_name} marked a Mach transition for a launch already below Mach 1"
8374 );
8375 }
8376 }
8377
8378 #[test]
8379 fn mach_transition_tracker_requires_a_downward_crossing() {
8380 fn record(mach_values: &[f64]) -> Vec<f64> {
8381 let mut tracker = MachTransitionTracker::default();
8382 let mut distances = Vec::new();
8383 for (index, mach) in mach_values.iter().copied().enumerate() {
8384 tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
8385 }
8386 distances
8387 }
8388
8389 assert!(record(&[0.9, 0.8, 0.7]).is_empty());
8390 assert_eq!(record(&[1.1, 1.05, 0.99]), vec![20.0]);
8391 assert_eq!(record(&[1.2, 1.19, 1.0, 0.99]), vec![10.0, 30.0]);
8392 assert_eq!(record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]), vec![20.0, 30.0]);
8393 assert!(record(&[1.3, f64::NAN, 1.1]).is_empty());
8394 }
8395
8396 #[test]
8397 fn mach_transition_tracker_labels_0_9_without_touching_the_flat_vec() {
8398 fn record(mach_values: &[f64]) -> (Vec<f64>, MachTransitionTracker) {
8404 let mut tracker = MachTransitionTracker::default();
8405 let mut distances = Vec::new();
8406 for (index, mach) in mach_values.iter().copied().enumerate() {
8407 tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
8408 }
8409 (distances, tracker)
8410 }
8411
8412 let (distances, tracker) = record(&[0.9, 0.8, 0.7]);
8415 assert!(distances.is_empty()); assert_eq!(tracker.mach_1_2_distance_m, None);
8417 assert_eq!(tracker.mach_1_0_distance_m, None);
8418 assert_eq!(tracker.mach_0_9_distance_m, Some(10.0));
8419
8420 let (distances, tracker) = record(&[1.1, 1.05, 0.99]);
8422 assert_eq!(distances, vec![20.0]);
8423 assert_eq!(tracker.mach_1_2_distance_m, None);
8424 assert_eq!(tracker.mach_1_0_distance_m, Some(20.0));
8425 assert_eq!(tracker.mach_0_9_distance_m, None);
8426
8427 let (distances, tracker) = record(&[1.2, 1.19, 1.0, 0.99]);
8429 assert_eq!(distances, vec![10.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(10.0));
8431 assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
8432 assert_eq!(tracker.mach_0_9_distance_m, None);
8433
8434 let (distances, tracker) = record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]);
8437 assert_eq!(distances, vec![20.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(20.0));
8439 assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
8440 assert_eq!(tracker.mach_0_9_distance_m, Some(50.0));
8441 assert!(
8442 tracker.mach_1_2_distance_m < tracker.mach_1_0_distance_m
8443 && tracker.mach_1_0_distance_m < tracker.mach_0_9_distance_m,
8444 "labeled crossings must be strictly increasing downrange"
8445 );
8446
8447 let (distances, tracker) = record(&[1.3, f64::NAN, 1.1]);
8449 assert!(distances.is_empty());
8450 assert_eq!(tracker.mach_1_2_distance_m, None);
8451 assert_eq!(tracker.mach_1_0_distance_m, None);
8452 assert_eq!(tracker.mach_0_9_distance_m, None);
8453 }
8454
8455 #[test]
8456 fn humidity_percent_converts_and_clamps() {
8457 let mut i = BallisticInputs {
8459 humidity: 0.5,
8460 ..BallisticInputs::default()
8461 };
8462 assert!((i.humidity_percent() - 50.0).abs() < 1e-9, "0.5 -> 50%");
8463 i.humidity = 0.0;
8464 assert_eq!(i.humidity_percent(), 0.0);
8465 i.humidity = 1.0;
8466 assert_eq!(i.humidity_percent(), 100.0);
8467 i.humidity = 1.5; assert_eq!(i.humidity_percent(), 100.0);
8469 }
8470
8471 fn vertical_at(shot_azimuth: f64, range_m: f64) -> f64 {
8474 let inputs = BallisticInputs {
8475 muzzle_velocity: 800.0,
8476 bc_value: 0.5,
8477 bc_type: DragModel::G7,
8478 muzzle_angle: 0.02, enable_coriolis: true,
8480 latitude: Some(45.0),
8481 shot_azimuth,
8482 ground_threshold: f64::NEG_INFINITY, ..BallisticInputs::default()
8484 };
8485 let mut solver = TrajectorySolver::new(
8486 inputs,
8487 WindConditions::default(),
8488 AtmosphericConditions::default(),
8489 );
8490 solver.set_max_range(range_m + 50.0);
8491 let r = solver.solve().expect("solve");
8492 let pts = &r.points;
8493 for i in 1..pts.len() {
8494 if pts[i].position.x >= range_m {
8495 let p1 = &pts[i - 1];
8496 let p2 = &pts[i];
8497 let t = (range_m - p1.position.x) / (p2.position.x - p1.position.x);
8498 return p1.position.y + t * (p2.position.y - p1.position.y);
8499 }
8500 }
8501 panic!("range {range_m} not reached");
8502 }
8503
8504 #[test]
8509 fn eotvos_east_higher_than_west() {
8510 let range = 600.0;
8511 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!(
8515 east > west,
8516 "east ({east:.5}) must be higher than west ({west:.5}) at {range} m (Eötvös)"
8517 );
8518 assert!(
8519 east > north && north > west,
8520 "north ({north:.5}) must lie between east ({east:.5}) and west ({west:.5})"
8521 );
8522 assert!(
8523 (east - west) > 1e-3,
8524 "E-W vertical separation ({:.6} m) should be physically meaningful, not FP noise",
8525 east - west
8526 );
8527 }
8528
8529 #[test]
8537 fn labeled_mach_crossings_match_pinned_pre_change_flat_vec_across_solvers() {
8538 let cases = [
8540 (
8541 "Euler",
8542 false,
8543 false,
8544 670.9878683238721_f64,
8545 805.5274119916264_f64,
8546 ),
8547 (
8548 "RK4",
8549 true,
8550 false,
8551 671.7257336844475_f64,
8552 805.933409072171_f64,
8553 ),
8554 (
8555 "RK45",
8556 true,
8557 true,
8558 672.4905711917901_f64,
8559 806.5709746782849_f64,
8560 ),
8561 ];
8562
8563 for (solver_name, use_rk4, use_adaptive_rk45, expected_1_2, expected_1_0) in cases {
8564 let inputs = BallisticInputs {
8565 muzzle_velocity: 850.0,
8566 bc_value: 0.2,
8567 bc_type: DragModel::G7,
8568 muzzle_angle: 0.03,
8569 use_rk4,
8570 use_adaptive_rk45,
8571 ..BallisticInputs::default()
8572 };
8573 let mut solver = TrajectorySolver::new(
8574 inputs,
8575 WindConditions::default(),
8576 AtmosphericConditions::default(),
8577 );
8578 solver.set_max_range(2000.0);
8579 let result = solver.solve().expect("solve should succeed");
8580
8581 assert_eq!(
8582 result.mach_1_2_distance_m,
8583 Some(expected_1_2),
8584 "{solver_name}: mach_1_2_distance_m must match the pinned pre-change flat-Vec value"
8585 );
8586 assert_eq!(
8587 result.mach_1_0_distance_m,
8588 Some(expected_1_0),
8589 "{solver_name}: mach_1_0_distance_m must match the pinned pre-change flat-Vec value"
8590 );
8591
8592 let mach_1_2 = result.mach_1_2_distance_m.expect("crosses 1.2");
8593 let mach_1_0 = result.mach_1_0_distance_m.expect("crosses 1.0");
8594 let mach_0_9 = result
8595 .mach_0_9_distance_m
8596 .expect("this trajectory also goes past 0.9 within 2000 m");
8597 assert!(
8598 mach_1_2 < mach_1_0 && mach_1_0 < mach_0_9,
8599 "{solver_name}: labeled crossings must be strictly increasing downrange \
8600 (1.2={mach_1_2}, 1.0={mach_1_0}, 0.9={mach_0_9})"
8601 );
8602 }
8603 }
8604
8605 #[test]
8608 fn labeled_mach_crossings_are_none_for_a_fully_supersonic_trajectory() {
8609 for (solver_name, use_rk4, use_adaptive_rk45) in [
8610 ("Euler", false, false),
8611 ("RK4", true, false),
8612 ("RK45", true, true),
8613 ] {
8614 let inputs = BallisticInputs {
8615 muzzle_velocity: 850.0,
8616 bc_value: 0.2,
8617 bc_type: DragModel::G7,
8618 muzzle_angle: 0.03,
8619 use_rk4,
8620 use_adaptive_rk45,
8621 ..BallisticInputs::default()
8622 };
8623 let mut solver = TrajectorySolver::new(
8624 inputs,
8625 WindConditions::default(),
8626 AtmosphericConditions::default(),
8627 );
8628 solver.set_max_range(200.0);
8630 let result = solver.solve().expect("solve should succeed");
8631
8632 assert_eq!(
8633 result.mach_1_2_distance_m, None,
8634 "{solver_name}: must not report a 1.2 crossing that never happens"
8635 );
8636 assert_eq!(
8637 result.mach_1_0_distance_m, None,
8638 "{solver_name}: must not report a 1.0 crossing that never happens"
8639 );
8640 assert_eq!(
8641 result.mach_0_9_distance_m, None,
8642 "{solver_name}: must not report a 0.9 crossing that never happens"
8643 );
8644 }
8645 }
8646}
8647
8648#[cfg(test)]
8649mod cant_tests {
8650 use super::*;
8651
8652 fn base_inputs() -> BallisticInputs {
8653 BallisticInputs {
8654 muzzle_velocity: 800.0,
8655 bc_value: 0.5,
8656 bc_type: DragModel::G7,
8657 bullet_mass: 0.0109,
8658 bullet_diameter: 0.00782,
8659 bullet_length: 0.0309,
8660 sight_height: 0.05,
8661 twist_rate: 10.0,
8662 use_rk4: true,
8663 ..BallisticInputs::default()
8664 }
8665 }
8666
8667 fn solve_with(inputs: BallisticInputs, max_range: f64) -> TrajectoryResult {
8668 let mut s = TrajectorySolver::new(
8669 inputs,
8670 WindConditions::default(),
8671 AtmosphericConditions::default(),
8672 );
8673 s.set_max_range(max_range);
8674 s.solve().expect("solve")
8675 }
8676
8677 fn yz_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
8679 let pts = &result.points;
8680 for i in 1..pts.len() {
8681 if pts[i].position.x >= x {
8682 let (p1, p2) = (&pts[i - 1], &pts[i]);
8683 let dx = p2.position.x - p1.position.x;
8684 let t = if dx.abs() < 1e-12 {
8685 0.0
8686 } else {
8687 (x - p1.position.x) / dx
8688 };
8689 return (
8690 p1.position.y + t * (p2.position.y - p1.position.y),
8691 p1.position.z + t * (p2.position.z - p1.position.z),
8692 );
8693 }
8694 }
8695 panic!("trajectory never reached {x} m");
8696 }
8697
8698 #[test]
8699 fn cant_sign_clockwise_up_offset_goes_right_and_low() {
8700 let mut level = base_inputs();
8702 level.muzzle_angle = 0.003; let mut canted = level.clone();
8704 canted.cant_angle = 10f64.to_radians();
8705
8706 let (y0, z0) = yz_at(&solve_with(level, 400.0), 300.0);
8707 let (y1, z1) = yz_at(&solve_with(canted, 400.0), 300.0);
8708 assert!(
8709 z1 > z0 + 0.01,
8710 "clockwise cant must move POI right: z0={z0} z1={z1}"
8711 );
8712 assert!(
8713 y1 < y0 - 0.001,
8714 "clockwise cant must move POI low: y0={y0} y1={y1}"
8715 );
8716 }
8717
8718 #[test]
8719 fn pure_cant_shows_bore_offset_near_range() {
8720 let mut i = base_inputs();
8723 i.muzzle_angle = 0.0;
8724 i.cant_angle = 10f64.to_radians();
8725 let sh = i.sight_height;
8726 let r = solve_with(i, 60.0);
8727 let first = &r.points[1]; let expected = -sh * 10f64.to_radians().sin();
8729 assert!(
8730 (first.position.z - expected).abs() < 0.005,
8731 "near-muzzle lateral {} should be ~bore offset {expected}",
8732 first.position.z
8733 );
8734 }
8735
8736 #[test]
8737 fn zero_angle_is_independent_of_cant() {
8738 let a = base_inputs();
8739 let mut b = base_inputs();
8740 b.cant_angle = 15f64.to_radians();
8741 let za = calculate_zero_angle(a.clone(), 100.0, 0.0).expect("zero a");
8742 let zb = calculate_zero_angle(b.clone(), 100.0, 0.0).expect("zero b");
8743 assert_eq!(
8744 za.to_bits(),
8745 zb.to_bits(),
8746 "zeroing must ignore cant: {za} vs {zb}"
8747 );
8748 let _ = (a.cant_angle, b.cant_angle);
8750 }
8751
8752 #[test]
8753 fn nonfinite_cant_is_rejected() {
8754 let mut i = base_inputs();
8755 i.cant_angle = f64::NAN;
8756 let s = TrajectorySolver::new(
8757 i,
8758 WindConditions::default(),
8759 AtmosphericConditions::default(),
8760 );
8761 assert!(s.solve().is_err());
8762 }
8763
8764 #[test]
8765 fn incline_and_cant_compose_without_breaking() {
8766 let mut flat = base_inputs();
8768 flat.muzzle_angle = 0.003;
8769 flat.shooting_angle = 15f64.to_radians();
8770 let mut canted = flat.clone();
8771 canted.cant_angle = 10f64.to_radians();
8772 let (_, z_flat) = yz_at(&solve_with(flat, 400.0), 300.0);
8773 let (_, z_cant) = yz_at(&solve_with(canted, 400.0), 300.0);
8774 assert!(
8775 z_cant > z_flat,
8776 "cant must still deflect right on an incline"
8777 );
8778 }
8779}
8780
8781#[cfg(test)]
8782mod vertical_wind_tests {
8783 use super::*;
8784
8785 fn base_inputs() -> BallisticInputs {
8786 BallisticInputs {
8787 muzzle_velocity: 800.0,
8788 bc_value: 0.5,
8789 bc_type: DragModel::G7,
8790 bullet_mass: 0.0109,
8791 bullet_diameter: 0.00782,
8792 bullet_length: 0.0309,
8793 sight_height: 0.05,
8794 twist_rate: 10.0,
8795 use_rk4: true,
8796 ..BallisticInputs::default()
8797 }
8798 }
8799
8800 fn y_at(result: &TrajectoryResult, x: f64) -> f64 {
8802 let pts = &result.points;
8803 for i in 1..pts.len() {
8804 if pts[i].position.x >= x {
8805 let (p1, p2) = (&pts[i - 1], &pts[i]);
8806 let dx = p2.position.x - p1.position.x;
8807 let t = if dx.abs() < 1e-12 {
8808 0.0
8809 } else {
8810 (x - p1.position.x) / dx
8811 };
8812 return p1.position.y + t * (p2.position.y - p1.position.y);
8813 }
8814 }
8815 panic!("trajectory never reached {x} m");
8816 }
8817
8818 fn solve_with(
8819 inputs: BallisticInputs,
8820 wind: WindConditions,
8821 max_range: f64,
8822 ) -> TrajectoryResult {
8823 let mut s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
8824 s.set_max_range(max_range);
8825 s.solve().expect("solve")
8826 }
8827
8828 #[test]
8829 fn updraft_raises_poi_downrange() {
8830 let calm_inputs = base_inputs();
8833 let calm_wind = WindConditions::default();
8834 let updraft = WindConditions {
8835 vertical_speed: 5.0,
8836 ..Default::default()
8837 };
8838
8839 let calm = solve_with(calm_inputs.clone(), calm_wind, 500.0);
8840 let updraft_result = solve_with(calm_inputs, updraft, 500.0);
8841
8842 let y_calm = y_at(&calm, 400.0);
8843 let y_updraft = y_at(&updraft_result, 400.0);
8844 assert!(
8845 y_updraft > y_calm,
8846 "5 m/s updraft must raise POI at 400m: calm={y_calm}, updraft={y_updraft}"
8847 );
8848 }
8849
8850 #[test]
8851 fn zero_vertical_is_default_and_finite_required() {
8852 assert_eq!(WindConditions::default().vertical_speed, 0.0);
8853
8854 let inputs = base_inputs();
8855 let wind = WindConditions {
8856 vertical_speed: f64::NAN,
8857 ..Default::default()
8858 };
8859 let s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
8860 assert!(
8861 s.solve().is_err(),
8862 "NaN wind.vertical_speed must be rejected by validate_for_solve"
8863 );
8864 }
8865}
8866
8867#[cfg(test)]
8869mod bc_reference_standard_tests {
8870 use super::*;
8871
8872 fn base_inputs() -> BallisticInputs {
8873 BallisticInputs {
8874 muzzle_velocity: 800.0,
8875 bc_value: 0.5,
8876 bc_type: DragModel::G7,
8877 bullet_mass: 0.0109,
8878 bullet_diameter: 0.00782,
8879 bullet_length: 0.0309,
8880 sight_height: 0.05,
8881 twist_rate: 10.0,
8882 use_rk4: true,
8883 ..BallisticInputs::default()
8884 }
8885 }
8886
8887 fn y_and_speed_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
8889 let pts = &result.points;
8890 for i in 1..pts.len() {
8891 if pts[i].position.x >= x {
8892 let (p1, p2) = (&pts[i - 1], &pts[i]);
8893 let dx = p2.position.x - p1.position.x;
8894 let t = if dx.abs() < 1e-12 {
8895 0.0
8896 } else {
8897 (x - p1.position.x) / dx
8898 };
8899 return (
8900 p1.position.y + t * (p2.position.y - p1.position.y),
8901 p1.velocity_magnitude + t * (p2.velocity_magnitude - p1.velocity_magnitude),
8902 );
8903 }
8904 }
8905 panic!("trajectory never reached {x} m");
8906 }
8907
8908 #[test]
8911 fn asm_to_icao_ratio_matches_documented_value() {
8912 assert!(
8913 (crate::constants::ASM_TO_ICAO_BC - 0.98237).abs() < 1e-5,
8914 "ASM_TO_ICAO_BC = {} must equal 0.98237 to 5 decimal places",
8915 crate::constants::ASM_TO_ICAO_BC
8916 );
8917 assert_eq!(
8922 crate::constants::ASM_TO_ICAO_BC,
8923 crate::constants::ASM_DENSITY_LB_FT3 / crate::constants::ICAO_DENSITY_LB_FT3
8924 );
8925 }
8926
8927 #[test]
8930 fn default_bc_reference_standard_is_icao() {
8931 assert_eq!(
8932 BallisticInputs::default().bc_reference_standard,
8933 BcReferenceStandard::Icao
8934 );
8935 }
8936
8937 #[test]
8942 fn icao_reference_leaves_bc_value_bit_identical() {
8943 let raw_bc: f64 = 0.4372911; let inputs = BallisticInputs {
8945 bc_value: raw_bc,
8946 bc_reference_standard: BcReferenceStandard::Icao,
8947 ..base_inputs()
8948 };
8949 let solver = TrajectorySolver::new(
8950 inputs,
8951 WindConditions::default(),
8952 AtmosphericConditions::default(),
8953 );
8954 assert_eq!(solver.inputs.bc_value.to_bits(), raw_bc.to_bits());
8955 }
8956
8957 #[test]
8958 fn default_inputs_solve_is_unaffected_by_the_new_field_existing() {
8959 let a = TrajectorySolver::new(
8963 base_inputs(),
8964 WindConditions::default(),
8965 AtmosphericConditions::default(),
8966 )
8967 .solve()
8968 .expect("solve a");
8969 let b = TrajectorySolver::new(
8970 BallisticInputs { ..base_inputs() },
8971 WindConditions::default(),
8972 AtmosphericConditions::default(),
8973 )
8974 .solve()
8975 .expect("solve b");
8976 assert_eq!(a.impact_velocity.to_bits(), b.impact_velocity.to_bits());
8977 assert_eq!(a.max_range.to_bits(), b.max_range.to_bits());
8978 }
8979
8980 #[test]
8983 fn army_standard_metro_scales_bc_value_by_exactly_the_derived_ratio() {
8984 let raw_bc = 0.5;
8985 let inputs = BallisticInputs {
8986 bc_value: raw_bc,
8987 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8988 ..base_inputs()
8989 };
8990 let solver = TrajectorySolver::new(
8991 inputs,
8992 WindConditions::default(),
8993 AtmosphericConditions::default(),
8994 );
8995 assert_eq!(
8996 solver.inputs.bc_value,
8997 raw_bc * crate::constants::ASM_TO_ICAO_BC
8998 );
8999 }
9000
9001 #[test]
9002 fn army_standard_metro_scales_mach_keyed_bc_segments() {
9003 let inputs = BallisticInputs {
9004 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
9005 bc_segments: Some(vec![(0.5, 0.40), (1.5, 0.30)]),
9006 ..base_inputs()
9007 };
9008 let solver = TrajectorySolver::new(
9009 inputs,
9010 WindConditions::default(),
9011 AtmosphericConditions::default(),
9012 );
9013 let segments = solver.inputs.bc_segments.as_ref().expect("segments");
9014 assert_eq!(segments[0], (0.5, 0.40 * crate::constants::ASM_TO_ICAO_BC));
9015 assert_eq!(segments[1], (1.5, 0.30 * crate::constants::ASM_TO_ICAO_BC));
9016 }
9017
9018 #[test]
9019 fn army_standard_metro_scales_velocity_keyed_bc_segments_data() {
9020 let inputs = BallisticInputs {
9021 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
9022 bc_segments_data: Some(vec![
9023 crate::BCSegmentData {
9024 velocity_min: 0.0,
9025 velocity_max: 500.0,
9026 bc_value: 0.40,
9027 },
9028 crate::BCSegmentData {
9029 velocity_min: 500.0,
9030 velocity_max: 900.0,
9031 bc_value: 0.45,
9032 },
9033 ]),
9034 ..base_inputs()
9035 };
9036 let solver = TrajectorySolver::new(
9037 inputs,
9038 WindConditions::default(),
9039 AtmosphericConditions::default(),
9040 );
9041 let segments = solver.inputs.bc_segments_data.as_ref().expect("segments");
9042 assert_eq!(
9043 segments[0].bc_value,
9044 0.40 * crate::constants::ASM_TO_ICAO_BC
9045 );
9046 assert_eq!(
9047 segments[1].bc_value,
9048 0.45 * crate::constants::ASM_TO_ICAO_BC
9049 );
9050 assert_eq!(segments[0].velocity_min, 0.0);
9052 assert_eq!(segments[1].velocity_max, 900.0);
9053 }
9054
9055 #[test]
9060 fn army_standard_metro_moves_impact_in_the_more_drag_direction() {
9061 let solve_at = |standard: BcReferenceStandard| {
9062 let inputs = BallisticInputs {
9063 bc_value: 0.475,
9064 bc_reference_standard: standard,
9065 ..base_inputs()
9066 };
9067 let mut solver = TrajectorySolver::new(
9068 inputs,
9069 WindConditions::default(),
9070 AtmosphericConditions::default(),
9071 );
9072 solver.set_max_range(500.0);
9073 solver.solve().expect("solve")
9074 };
9075
9076 let icao = solve_at(BcReferenceStandard::Icao);
9077 let asm = solve_at(BcReferenceStandard::ArmyStandardMetro);
9078
9079 let (y_icao, v_icao) = y_and_speed_at(&icao, 400.0);
9080 let (y_asm, v_asm) = y_and_speed_at(&asm, 400.0);
9081
9082 assert!(
9083 y_asm < y_icao,
9084 "ArmyStandardMetro must drop MORE (lower y) at 400m than Icao for the same raw \
9085 bc_value: icao_y={y_icao}, asm_y={y_asm}"
9086 );
9087 assert!(
9088 v_asm < v_icao,
9089 "ArmyStandardMetro must retain LESS velocity at 400m than Icao for the same raw \
9090 bc_value: icao_v={v_icao}, asm_v={v_asm}"
9091 );
9092 }
9093
9094 #[test]
9101 fn monte_carlo_inherits_the_normalized_bc_reference() {
9102 let base_inputs_asm = BallisticInputs {
9103 bc_value: 0.475,
9104 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
9105 ..base_inputs()
9106 };
9107 let wind = WindConditions::default();
9108
9109 let mut direct_solver = TrajectorySolver::new(
9115 base_inputs_asm.clone(),
9116 wind.clone(),
9117 AtmosphericConditions::default(),
9118 );
9119 direct_solver.set_max_range(base_inputs_asm.target_distance.max(1000.0) * 2.0);
9120 let direct = direct_solver.solve().expect("direct solve");
9121
9122 let mc_params = MonteCarloParams {
9123 num_simulations: 1,
9124 velocity_std_dev: 0.0,
9125 angle_std_dev: 0.0,
9126 bc_std_dev: 0.0,
9127 wind_speed_std_dev: 0.0,
9128 target_distance: None,
9129 base_wind_speed: 0.0,
9130 base_wind_direction: 0.0,
9131 azimuth_std_dev: 0.0,
9132 };
9133 let mc = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
9134 base_inputs_asm,
9135 wind,
9136 mc_params,
9137 0.0,
9138 42,
9139 )
9140 .expect("monte carlo");
9141
9142 assert_eq!(mc.ranges.len(), 1);
9143 assert_eq!(
9144 mc.ranges[0].to_bits(),
9145 direct.max_range.to_bits(),
9146 "a zero-dispersion single MC sample must match a plain solve of the same \
9147 ASM-referenced inputs bit-for-bit"
9148 );
9149 assert_eq!(
9150 mc.impact_velocities[0].to_bits(),
9151 direct.impact_velocity.to_bits()
9152 );
9153 }
9154
9155 #[test]
9163 fn estimate_bc_fit_recovers_an_icao_referenced_bc() {
9164 let known_bc = 0.475;
9165 let velocity = 800.0;
9166 let mass = 0.0109;
9167 let diameter = 0.00782;
9168 let atmosphere = AtmosphericConditions::default();
9169
9170 let synth_inputs = BallisticInputs {
9171 muzzle_velocity: velocity,
9172 bc_value: known_bc,
9173 bc_type: DragModel::G7,
9174 bullet_mass: mass,
9175 bullet_diameter: diameter,
9176 bullet_length: 0.0309,
9177 sight_height: 0.05,
9178 twist_rate: 10.0,
9179 use_rk4: true,
9180 bc_reference_standard: BcReferenceStandard::Icao,
9181 ..BallisticInputs::default()
9182 };
9183 let mut solver =
9184 TrajectorySolver::new(synth_inputs, WindConditions::default(), atmosphere.clone());
9185 solver.set_max_range(500.0);
9186 let trajectory = solver.solve().expect("synthetic solve");
9187
9188 let points: Vec<(f64, f64)> = [100.0, 200.0, 300.0, 400.0]
9189 .iter()
9190 .map(|&d| {
9191 let (y, _) = {
9192 let pts = &trajectory.points;
9193 let mut found = None;
9194 for i in 1..pts.len() {
9195 if pts[i].position.x >= d {
9196 let (p1, p2) = (&pts[i - 1], &pts[i]);
9197 let dx = p2.position.x - p1.position.x;
9198 let t = if dx.abs() < 1e-12 {
9199 0.0
9200 } else {
9201 (d - p1.position.x) / dx
9202 };
9203 found =
9204 Some((p1.position.y + t * (p2.position.y - p1.position.y), 0.0));
9205 break;
9206 }
9207 }
9208 found.expect("trajectory reached observation distance")
9209 };
9210 (d, -y) })
9212 .collect();
9213
9214 let estimate = estimate_bc_fit(
9215 velocity,
9216 mass,
9217 diameter,
9218 &points,
9219 DragModel::G7,
9220 BcFitMode::Drop,
9221 atmosphere,
9222 None,
9223 0.05,
9224 )
9225 .expect("fit should converge");
9226
9227 assert!(
9228 (estimate.bc - known_bc).abs() < 0.02,
9229 "fit should recover the known ICAO-referenced bc={known_bc}, got {}",
9230 estimate.bc
9231 );
9232 }
9233
9234 #[test]
9237 fn custom_drag_table_makes_bc_reference_standard_numerically_inert() {
9238 let table =
9239 crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0, 3.0], vec![0.3, 0.4, 0.3, 0.2])
9240 .expect("valid table");
9241
9242 let solve_with = |standard: BcReferenceStandard| {
9243 let inputs = BallisticInputs {
9244 bc_value: 0.5, bc_reference_standard: standard,
9246 custom_drag_table: Some(table.clone()),
9247 ..base_inputs()
9248 };
9249 let mut solver = TrajectorySolver::new(
9250 inputs,
9251 WindConditions::default(),
9252 AtmosphericConditions::default(),
9253 );
9254 solver.set_max_range(500.0);
9255 solver.solve().expect("solve")
9256 };
9257
9258 let icao = solve_with(BcReferenceStandard::Icao);
9259 let asm = solve_with(BcReferenceStandard::ArmyStandardMetro);
9260
9261 assert_eq!(
9262 icao.impact_velocity.to_bits(),
9263 asm.impact_velocity.to_bits(),
9264 "a custom drag table must make bc_reference_standard fully inert"
9265 );
9266 assert_eq!(icao.max_range.to_bits(), asm.max_range.to_bits());
9267 }
9268
9269 #[test]
9270 fn custom_drag_table_inert_warning_fires_only_for_army_standard_metro_with_a_table() {
9271 let table = crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0], vec![0.3, 0.4, 0.3])
9272 .expect("valid table");
9273
9274 let no_table_icao = base_inputs();
9276 assert!(no_table_icao
9277 .bc_reference_standard_inert_warning()
9278 .is_none());
9279 let no_table_asm = BallisticInputs {
9280 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
9281 ..base_inputs()
9282 };
9283 assert!(no_table_asm.bc_reference_standard_inert_warning().is_none());
9284
9285 let table_icao = BallisticInputs {
9287 custom_drag_table: Some(table.clone()),
9288 ..base_inputs()
9289 };
9290 assert!(table_icao.bc_reference_standard_inert_warning().is_none());
9291
9292 let table_asm = BallisticInputs {
9294 custom_drag_table: Some(table),
9295 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
9296 ..base_inputs()
9297 };
9298 let warning = table_asm
9299 .bc_reference_standard_inert_warning()
9300 .expect("must warn");
9301 assert!(warning.contains("--bc-reference"));
9302 assert!(warning.contains("--drag-table"));
9303 }
9304}
9305
9306#[cfg(test)]
9312mod effective_drag_coefficient_tests {
9313 use super::*;
9314
9315 fn inputs_175gr_g7() -> BallisticInputs {
9316 let mut inputs = BallisticInputs {
9317 bc_value: 0.243,
9318 bc_type: DragModel::G7,
9319 muzzle_velocity: 823.0,
9320 ..Default::default()
9321 };
9322 inputs.bullet_mass = 175.0 * crate::constants::GRAINS_TO_KG;
9326 inputs.bullet_diameter = 0.308 * 0.0254;
9327 inputs.weight_grains = 175.0;
9328 inputs.caliber_inches = 0.308;
9329 inputs
9330 }
9331
9332 fn solver(inputs: BallisticInputs) -> TrajectorySolver {
9333 TrajectorySolver::new(
9334 inputs,
9335 WindConditions::default(),
9336 AtmosphericConditions::default(),
9337 )
9338 }
9339
9340 #[test]
9344 fn reports_the_projectiles_own_cd_not_the_reference_tables() {
9345 let inputs = inputs_175gr_g7();
9346 let sd = inputs.sectional_density_lb_in2().expect("SD");
9347 let solver = solver(inputs);
9348
9349 let sos = 340.0;
9350 let velocity = 800.0;
9351 let mach = velocity / sos;
9352
9353 let reference = crate::drag::get_drag_coefficient(mach, &DragModel::G7);
9354 let reported = solver
9355 .effective_drag_coefficient(velocity, sos)
9356 .expect("mass and diameter are set");
9357
9358 let expected = reference * sd / 0.243;
9359 assert!(
9360 (reported - expected).abs() < 1e-12,
9361 "reported {reported} != Cd_ref * SD / BC {expected}"
9362 );
9363 assert!(
9366 (reported - reference).abs() > 1e-6,
9367 "form factor collapsed to 1; this fixture no longer distinguishes the two values"
9368 );
9369 }
9370
9371 #[test]
9374 fn a_custom_drag_table_passes_through_unscaled() {
9375 let mut inputs = inputs_175gr_g7();
9376 inputs.custom_drag_table = Some(crate::drag::DragTable::new(
9377 vec![0.5, 3.0],
9378 vec![0.15, 0.40],
9379 ));
9380 let solver = solver(inputs);
9381
9382 let sos = 340.0;
9383 let velocity = 0.9 * sos;
9384 let table_value = solver
9385 .inputs
9386 .custom_drag_table
9387 .as_ref()
9388 .expect("table")
9389 .interpolate(0.9);
9390
9391 let reported = solver
9392 .effective_drag_coefficient(velocity, sos)
9393 .expect("mass and diameter are set");
9394 assert!(
9395 (reported - table_value).abs() < 1e-12,
9396 "custom table Cd {table_value} was rescaled to {reported}"
9397 );
9398 }
9399
9400 #[test]
9403 fn a_velocity_segmented_bc_steps_the_reported_cd() {
9404 let mut inputs = inputs_175gr_g7();
9405 inputs.use_bc_segments = true;
9406 inputs.bc_segments_data = Some(vec![
9407 crate::BCSegmentData {
9408 velocity_min: 2400.0,
9409 velocity_max: 4000.0,
9410 bc_value: 0.243,
9411 },
9412 crate::BCSegmentData {
9413 velocity_min: 0.0,
9414 velocity_max: 2400.0,
9415 bc_value: 0.200,
9416 },
9417 ]);
9418 let solver = solver(inputs);
9419
9420 let sos = 340.0;
9421 let above = solver
9423 .effective_drag_coefficient(2500.0 / 3.28084, sos)
9424 .expect("cd");
9425 let below = solver
9426 .effective_drag_coefficient(2300.0 / 3.28084, sos)
9427 .expect("cd");
9428
9429 assert!(
9431 below > above,
9432 "expected the 0.200 band to report a higher Cd than the 0.243 band; got {below} vs {above}"
9433 );
9434 }
9435
9436 #[test]
9439 fn is_absent_when_sectional_density_is_unknown() {
9440 let mut inputs = inputs_175gr_g7();
9441 inputs.weight_grains = 0.0;
9442 inputs.bullet_mass = 0.0;
9443 let solver = solver(inputs);
9444 assert!(solver.effective_drag_coefficient(800.0, 340.0).is_none());
9445 }
9446
9447 #[test]
9452 fn the_json_emit_rule_is_flag_gated_and_absent_when_cd_is_unknown() {
9453 let mut point = TrajectoryPoint {
9454 time: 0.0,
9455 position: nalgebra::Vector3::new(0.0, 0.0, 0.0),
9456 velocity_magnitude: 800.0,
9457 kinetic_energy: 3000.0,
9458 drag_coefficient: Some(0.31),
9459 };
9460 assert_eq!(point.drag_coefficient_json_value(true), Some(0.31));
9461 assert_eq!(
9462 point.drag_coefficient_json_value(false),
9463 None,
9464 "without the flag the key must not exist, so default JSON stays byte-identical"
9465 );
9466 point.drag_coefficient = None;
9467 assert_eq!(
9468 point.drag_coefficient_json_value(true),
9469 None,
9470 "unknown sectional density must yield an ABSENT key, not null"
9471 );
9472 }
9473
9474 #[test]
9476 fn every_point_of_a_solved_trajectory_carries_the_value() {
9477 let mut solver = solver(inputs_175gr_g7());
9478 solver.set_max_range(300.0);
9479 let result = solver.solve().expect("solve");
9480 assert!(!result.points.is_empty());
9481 assert!(
9482 result.points.iter().all(|p| p.drag_coefficient.is_some()),
9483 "the post-integration pass missed at least one point"
9484 );
9485 }
9486}