1use crate::cluster_bc::ClusterBCDegradation;
3use crate::mc_stats::{
4 wilson_interval, BernoulliConfidenceSequence, ConfidenceLevel, Welford,
5};
6use crate::pitch_damping::{calculate_pitch_damping_coefficient, PitchDampingCoefficients};
7use crate::precession_nutation::{
8 calculate_combined_angular_motion, projectile_moments_of_inertia, AngularState,
9 PrecessionNutationParams,
10};
11use crate::trajectory_sampling::{
12 projected_sample_count, sample_trajectory, TrajectoryData, TrajectoryOutputs,
13 TrajectorySample,
14};
15use crate::trajectory_observation::{bracket_param, Bracket, TrajectoryTermination};
16use crate::wind_shear::WindShearModel;
17use crate::DragModel;
18use nalgebra::{Vector3, Vector6};
19use std::error::Error;
20use std::fmt;
21
22#[derive(Debug, Clone, Copy, PartialEq)]
32#[cfg_attr(feature = "cli", derive(clap::ValueEnum))]
33pub enum UnitSystem {
34 Metric,
36 Imperial,
38}
39
40#[derive(Debug, Clone, Copy, PartialEq)]
42pub enum OutputFormat {
43 Table,
44 Json,
45 Csv,
46}
47
48#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
67#[cfg_attr(feature = "cli", derive(clap::ValueEnum))]
68pub enum BcReferenceStandard {
69 #[default]
72 Icao,
73 ArmyStandardMetro,
76}
77
78pub const BC_REFERENCE_STANDARD_INERT_WARNING: &str =
82 "warning: --bc-reference army-standard-metro has no effect together with a custom drag \
83 table (--drag-table): the deck's Cd is divided by sectional density, not a BC value, so \
84 no BC-reference conversion applies";
85
86#[derive(Debug)]
88pub struct BallisticsError {
89 message: String,
90}
91
92impl fmt::Display for BallisticsError {
93 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
94 write!(f, "{}", self.message)
95 }
96}
97
98impl Error for BallisticsError {}
99
100impl From<String> for BallisticsError {
101 fn from(msg: String) -> Self {
102 BallisticsError { message: msg }
103 }
104}
105
106impl From<&str> for BallisticsError {
107 fn from(msg: &str) -> Self {
108 BallisticsError {
109 message: msg.to_string(),
110 }
111 }
112}
113
114#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
122pub enum DropsReference {
123 #[default]
125 Los,
126 Target,
129}
130
131#[derive(Debug, Clone)]
135pub struct BallisticInputs {
136 pub bc_value: f64, pub bc_type: DragModel, pub bc_reference_standard: BcReferenceStandard,
144 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,
158 pub shooting_angle: f64, pub cant_angle: f64,
168 pub sight_height: f64, pub sight_offset_lateral_m: f64,
180 pub muzzle_height: f64, pub target_height: f64, pub zero_poi_vertical_m: f64,
191 pub zero_poi_horizontal_m: f64,
197 pub ground_threshold: f64, pub altitude: f64, pub temperature: f64, pub pressure: f64, pub humidity: f64,
208 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,
230 pub enable_magnus: bool, pub enable_coriolis: bool, pub use_powder_sensitivity: bool,
233 pub powder_temp_sensitivity: f64, pub powder_temp: f64, pub powder_temp_curve: Option<Vec<(f64, f64)>>,
242 pub powder_curve_temp_c: Option<f64>,
246 pub tipoff_yaw: f64, pub tipoff_decay_distance: f64, pub use_bc_segments: bool,
251 pub bc_segments: Option<Vec<(f64, f64)>>, pub bc_segments_data: Option<Vec<crate::BCSegmentData>>, pub use_enhanced_spin_drift: bool,
254 pub use_form_factor: bool,
257 pub enable_wind_shear: bool,
258 pub wind_shear_model: String,
259 pub enable_trajectory_sampling: bool,
260 pub sample_interval: f64, pub drops_reference: DropsReference,
270 pub enable_pitch_damping: bool,
271 pub enable_precession_nutation: bool,
272 pub enable_aerodynamic_jump: bool,
275 pub use_cluster_bc: bool, pub custom_drag_table: Option<crate::drag::DragTable>,
279 pub cd_scale: f64,
287
288 pub bc_type_str: Option<String>,
290}
291
292impl BallisticInputs {
293 pub fn humidity_percent(&self) -> f64 {
298 (self.humidity * 100.0).clamp(0.0, 100.0)
299 }
300
301 pub fn windage_zero_bias_rad(&self, zero_distance_m: f64) -> f64 {
316 if zero_distance_m > 0.0 {
317 (self.zero_poi_horizontal_m + self.sight_offset_lateral_m) / zero_distance_m
318 } else {
319 0.0
320 }
321 }
322
323 pub fn sectional_density_lb_in2(&self) -> Option<f64> {
329 let weight_gr = if self.weight_grains > 0.0 {
330 self.weight_grains
331 } else {
332 self.bullet_mass / crate::constants::GRAINS_TO_KG };
334 let diameter_in = if self.caliber_inches > 0.0 {
335 self.caliber_inches
336 } else {
337 self.bullet_diameter / 0.0254 };
339 if weight_gr > 0.0 && diameter_in > 0.0 {
340 Some(weight_gr / 7000.0 / (diameter_in * diameter_in))
341 } else {
342 None
343 }
344 }
345
346 pub fn custom_drag_denominator(&self, fallback_bc: f64) -> f64 {
358 match self.sectional_density_lb_in2() {
359 Some(sd) => sd,
360 None => {
361 static WARN_ONCE: std::sync::Once = std::sync::Once::new();
362 WARN_ONCE.call_once(|| {
363 eprintln!(
364 "Warning: custom drag table active but bullet mass/diameter are \
365 unavailable; falling back to bc_value for the retardation denominator"
366 );
367 });
368 fallback_bc
369 }
370 }
371 }
372
373 pub fn bc_reference_standard_inert_warning(&self) -> Option<&'static str> {
384 if self.custom_drag_table.is_some()
385 && matches!(self.bc_reference_standard, BcReferenceStandard::ArmyStandardMetro)
386 {
387 Some(BC_REFERENCE_STANDARD_INERT_WARNING)
388 } else {
389 None
390 }
391 }
392
393 pub fn normalize_for_solve(&mut self) {
412 if matches!(
423 self.bc_reference_standard,
424 BcReferenceStandard::ArmyStandardMetro
425 ) {
426 self.bc_value *= crate::constants::ASM_TO_ICAO_BC;
427 if let Some(segments) = self.bc_segments.as_mut() {
428 for (_mach, bc) in segments.iter_mut() {
429 *bc *= crate::constants::ASM_TO_ICAO_BC;
430 }
431 }
432 if let Some(segments) = self.bc_segments_data.as_mut() {
433 for segment in segments.iter_mut() {
434 segment.bc_value *= crate::constants::ASM_TO_ICAO_BC;
435 }
436 }
437 self.bc_reference_standard = BcReferenceStandard::Icao;
440 }
441
442 self.caliber_inches = self.bullet_diameter / 0.0254;
447 self.weight_grains = self.bullet_mass / crate::constants::GRAINS_TO_KG;
448
449 self.muzzle_velocity = resolve_powder_adjusted_velocity(
460 self.muzzle_velocity,
461 self.temperature,
462 self.use_powder_sensitivity,
463 self.powder_temp_sensitivity,
464 self.powder_temp,
465 self.powder_temp_curve.as_deref(),
466 self.powder_curve_temp_c,
467 );
468 }
469}
470
471impl Default for BallisticInputs {
472 fn default() -> Self {
473 let mass_kg = 0.01;
474 let diameter_m = 0.00762;
475 let bc = 0.5;
476 let muzzle_angle_rad = 0.0;
477 let bc_type = DragModel::G1;
478
479 Self {
480 bc_value: bc,
482 bc_type,
483 bc_reference_standard: BcReferenceStandard::Icao,
484 bullet_mass: mass_kg,
485 muzzle_velocity: 800.0,
486 bullet_diameter: diameter_m,
487 bullet_length: crate::stability::estimate_bullet_length_m(diameter_m, mass_kg),
491
492 muzzle_angle: muzzle_angle_rad,
494 target_distance: 100.0,
495 azimuth_angle: 0.0,
496 shot_azimuth: 0.0,
497 shooting_angle: 0.0,
498 cant_angle: 0.0,
499 sight_height: 0.05,
500 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,
505 ground_threshold: -100.0, altitude: 0.0,
509 temperature: 15.0,
510 pressure: 1013.25, humidity: 0.5, latitude: None,
513
514 wind_speed: 0.0,
516 wind_angle: 0.0,
517
518 twist_rate: 12.0, is_twist_right: true,
521 caliber_inches: diameter_m / 0.0254, weight_grains: mass_kg / crate::constants::GRAINS_TO_KG, manufacturer: None,
524 bullet_model: None,
525 bullet_id: None,
526 bullet_cluster: None,
527
528 use_rk4: true, use_adaptive_rk45: true, enable_advanced_effects: false,
534 enable_magnus: false,
535 enable_coriolis: false,
536 use_powder_sensitivity: false,
537 powder_temp_sensitivity: 0.0,
538 powder_temp: 15.0,
539 powder_temp_curve: None,
540 powder_curve_temp_c: None,
541 tipoff_yaw: 0.0,
542 tipoff_decay_distance: 50.0,
543 use_bc_segments: false,
544 bc_segments: None,
545 bc_segments_data: None,
546 use_enhanced_spin_drift: false,
547 use_form_factor: false,
548 enable_wind_shear: false,
549 wind_shear_model: "none".to_string(),
550 enable_trajectory_sampling: false,
551 sample_interval: 10.0, drops_reference: DropsReference::Los, enable_pitch_damping: false,
554 enable_precession_nutation: false,
555 enable_aerodynamic_jump: false,
556 use_cluster_bc: false, custom_drag_table: None,
560 cd_scale: 1.0,
561
562 bc_type_str: None,
564 }
565 }
566}
567
568pub fn interpolate_powder_temp_curve(curve: &[(f64, f64)], temp_c: f64) -> f64 {
574 debug_assert!(!curve.is_empty());
575 if curve.is_empty() {
576 return 0.0;
577 }
578 let mut sorted;
581 let pts: &[(f64, f64)] = if curve.windows(2).all(|w| w[0].0 <= w[1].0) {
582 curve
583 } else {
584 sorted = curve.to_vec();
585 sorted.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
586 &sorted
587 };
588 let n = pts.len();
589 if temp_c <= pts[0].0 {
590 return pts[0].1; }
592 if temp_c >= pts[n - 1].0 {
593 return pts[n - 1].1; }
595 for i in 1..n {
596 let (t0, v0) = pts[i - 1];
597 let (t1, v1) = pts[i];
598 if temp_c <= t1 {
599 let span = t1 - t0;
600 if span.abs() < f64::EPSILON {
601 return v1; }
603 let f = (temp_c - t0) / span;
604 return v0 + f * (v1 - v0);
605 }
606 }
607 pts[n - 1].1
608}
609
610pub fn parse_powder_sweep(s: &str) -> Result<Vec<f64>, String> {
616 const MAX_SWEEP_ROWS: usize = 500;
617 let parts: Vec<&str> = s.split(':').collect();
618 if parts.len() != 3 {
619 return Err(format!(
620 "Invalid --sweep '{}': expected START:END:STEP (e.g. \"20:110:10\")",
621 s
622 ));
623 }
624 let parse = |p: &str, name: &str| -> Result<f64, String> {
625 p.trim()
626 .parse::<f64>()
627 .map_err(|_| format!("Invalid --sweep {}: '{}' is not a number", name, p.trim()))
628 };
629 let start = parse(parts[0], "START")?;
630 let end = parse(parts[1], "END")?;
631 let step = parse(parts[2], "STEP")?;
632 if !step.is_finite() || step <= 0.0 {
633 return Err(format!("Invalid --sweep STEP {}: must be positive", step));
634 }
635 if !start.is_finite() || !end.is_finite() || end < start {
636 return Err(format!(
637 "Invalid --sweep range {}:{}: END must be >= START",
638 start, end
639 ));
640 }
641 let n_f = ((end - start) / step + 1e-9).floor();
647 if !n_f.is_finite() || n_f + 1.0 > MAX_SWEEP_ROWS as f64 {
648 return Err(format!(
649 "--sweep would produce more than {} rows; use a larger STEP",
650 MAX_SWEEP_ROWS
651 ));
652 }
653 let n = n_f as usize + 1;
654 Ok((0..n).map(|i| start + step * i as f64).collect())
656}
657
658pub fn resolve_powder_adjusted_velocity(
667 nominal_velocity_mps: f64,
668 ambient_temperature_c: f64,
669 use_powder_sensitivity: bool,
670 powder_temp_sensitivity_mps_per_c: f64,
671 powder_reference_temp_c: f64,
672 powder_temp_curve: Option<&[(f64, f64)]>,
673 powder_curve_temp_c: Option<f64>,
674) -> f64 {
675 if let Some(curve) = powder_temp_curve {
676 if !curve.is_empty() {
677 let lookup_c = powder_curve_temp_c.unwrap_or(ambient_temperature_c);
678 return interpolate_powder_temp_curve(curve, lookup_c);
679 }
680 return nominal_velocity_mps;
683 }
684 if use_powder_sensitivity {
685 let temp_delta_c = ambient_temperature_c - powder_reference_temp_c;
686 return nominal_velocity_mps + powder_temp_sensitivity_mps_per_c * temp_delta_c;
687 }
688 nominal_velocity_mps
689}
690
691#[derive(Debug, Clone)]
693pub struct WindConditions {
694 pub speed: f64, pub direction: f64,
698 pub vertical_speed: f64,
706}
707
708impl Default for WindConditions {
709 fn default() -> Self {
710 Self {
711 speed: 0.0,
712 direction: 0.0,
713 vertical_speed: 0.0,
714 }
715 }
716}
717
718#[derive(Debug, Clone)]
720pub struct AtmosphericConditions {
721 pub temperature: f64, pub pressure: f64, pub humidity: f64,
727 pub altitude: f64, }
729
730impl Default for AtmosphericConditions {
731 fn default() -> Self {
732 Self {
733 temperature: 15.0,
734 pressure: 1013.25,
735 humidity: 50.0,
736 altitude: 0.0,
737 }
738 }
739}
740
741#[derive(Debug, Clone)]
743pub struct TrajectoryPoint {
744 pub time: f64,
745 pub position: Vector3<f64>,
746 pub velocity_magnitude: f64,
747 pub kinetic_energy: f64,
748 pub drag_coefficient: Option<f64>,
755}
756
757impl TrajectoryPoint {
758 pub fn drag_coefficient_json_value(&self, with_drag_coefficient: bool) -> Option<f64> {
766 if with_drag_coefficient {
767 self.drag_coefficient
768 } else {
769 None
770 }
771 }
772}
773
774#[derive(Debug, Clone)]
776pub struct TrajectoryResult {
777 pub max_range: f64,
778 pub max_height: f64,
779 pub time_of_flight: f64,
780 pub impact_velocity: f64,
781 pub impact_energy: f64,
782 pub projectile_mass_kg: f64,
784 pub line_of_sight_height_m: f64,
786 pub station_speed_of_sound_mps: f64,
788 pub termination: TrajectoryTermination,
790 pub points: Vec<TrajectoryPoint>,
791 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>,
800 pub mach_1_2_distance_m: Option<f64>,
805 pub mach_1_0_distance_m: Option<f64>,
809 pub mach_0_9_distance_m: Option<f64>,
817}
818
819const RK45_TOLERANCE: f64 = 1e-6;
820const RK45_SAFETY_FACTOR: f64 = 0.9;
821const RK45_MAX_DT: f64 = 0.01;
822const RK45_MIN_DT: f64 = 1e-6;
823const TRAJECTORY_TIME_LIMIT_S: f64 = 100.0;
824
825pub const MAX_TRAJECTORY_POINTS: usize = 250_000;
831
832fn cli_rk45_error_norm(
834 position: &Vector3<f64>,
835 velocity: &Vector3<f64>,
836 fifth_position: &Vector3<f64>,
837 fifth_velocity: &Vector3<f64>,
838 fourth_position: &Vector3<f64>,
839 fourth_velocity: &Vector3<f64>,
840) -> f64 {
841 let pack_state = |position: &Vector3<f64>, velocity: &Vector3<f64>| {
842 Vector6::new(
843 position.x, position.y, position.z, velocity.x, velocity.y, velocity.z,
844 )
845 };
846 let state = pack_state(position, velocity);
847 let fifth_order = pack_state(fifth_position, fifth_velocity);
848 let fourth_order = pack_state(fourth_position, fourth_velocity);
849
850 crate::trajectory_integration::rk45_error_norm(&state, &fifth_order, &fourth_order)
851}
852
853struct Rk45Trial {
854 position: Vector3<f64>,
855 velocity: Vector3<f64>,
856 suggested_dt: f64,
857 error: f64,
858}
859
860struct Rk45AcceptedStep {
861 position: Vector3<f64>,
862 velocity: Vector3<f64>,
863 used_dt: f64,
864 next_dt: f64,
865 error: f64,
866}
867
868#[derive(Default)]
882struct MachTransitionTracker {
883 previous_mach: Option<f64>,
884 crossed_transonic: bool,
885 crossed_subsonic: bool,
886 crossed_narrow: bool,
887 mach_1_2_distance_m: Option<f64>,
890 mach_1_0_distance_m: Option<f64>,
893 mach_0_9_distance_m: Option<f64>,
896}
897
898impl MachTransitionTracker {
899 fn record_downward_crossings(&mut self, mach: f64, downrange_m: f64, distances: &mut Vec<f64>) {
900 if !mach.is_finite() {
901 self.previous_mach = None;
902 return;
903 }
904
905 if let Some(previous_mach) = self.previous_mach {
906 if !self.crossed_transonic && previous_mach >= 1.2 && mach < 1.2 {
907 self.crossed_transonic = true;
908 distances.push(downrange_m);
909 self.mach_1_2_distance_m = Some(downrange_m);
910 }
911 if !self.crossed_subsonic && previous_mach >= 1.0 && mach < 1.0 {
912 self.crossed_subsonic = true;
913 distances.push(downrange_m);
914 self.mach_1_0_distance_m = Some(downrange_m);
915 }
916 if !self.crossed_narrow && previous_mach >= 0.9 && mach < 0.9 {
917 self.crossed_narrow = true;
918 self.mach_0_9_distance_m = Some(downrange_m);
920 }
921 }
922 self.previous_mach = Some(mach);
923 }
924}
925
926impl TrajectoryResult {
927 pub fn position_at_range(&self, target_range: f64) -> Option<Vector3<f64>> {
931 if self.points.is_empty() {
932 return None;
933 }
934
935 for i in 0..self.points.len() - 1 {
937 let p1 = &self.points[i];
938 let p2 = &self.points[i + 1];
939
940 if p1.position.x <= target_range && p2.position.x >= target_range {
942 let dx = p2.position.x - p1.position.x;
944 if dx.abs() < 1e-10 {
945 return Some(p1.position);
946 }
947 let t = (target_range - p1.position.x) / dx;
948
949 return Some(Vector3::new(
951 target_range,
952 p1.position.y + t * (p2.position.y - p1.position.y),
953 p1.position.z + t * (p2.position.z - p1.position.z),
954 ));
955 }
956 }
957
958 self.points.last().map(|p| p.position)
960 }
961}
962
963#[derive(Debug, Clone, Copy, PartialEq, Eq)]
965enum StationAtmosphereResolution {
966 LegacyDefaultSentinels,
969 Authoritative,
972}
973
974#[derive(Clone)]
975pub struct TrajectorySolver {
976 inputs: BallisticInputs,
977 wind: WindConditions,
978 atmosphere: AtmosphericConditions,
979 station_atmosphere_resolution: StationAtmosphereResolution,
980 max_range: f64,
981 time_step: f64,
982 max_trajectory_points: usize,
983 cluster_bc: Option<ClusterBCDegradation>,
984 precession_nutation_inertias: (f64, f64),
986 wind_sock: Option<crate::wind::WindSock>,
991 atmo_sock: Option<crate::atmosphere::AtmoSock>,
998}
999
1000#[derive(Debug, Clone, Copy, PartialEq, Eq)]
1012pub(crate) enum ZeroTargetFrame {
1013 SightLine,
1014 WorldVertical,
1015}
1016
1017#[derive(Debug, Clone, Copy, PartialEq)]
1030pub struct ZeroCrossings {
1031 pub near_m: Option<f64>,
1034 pub far_m: Option<f64>,
1037}
1038
1039impl TrajectorySolver {
1040 pub fn new(
1041 inputs: BallisticInputs,
1042 wind: WindConditions,
1043 atmosphere: AtmosphericConditions,
1044 ) -> Self {
1045 Self::new_with_station_atmosphere_resolution(
1046 inputs,
1047 wind,
1048 atmosphere,
1049 StationAtmosphereResolution::LegacyDefaultSentinels,
1050 )
1051 }
1052
1053 pub fn new_with_resolved_station_atmosphere(
1064 inputs: BallisticInputs,
1065 wind: WindConditions,
1066 atmosphere: AtmosphericConditions,
1067 ) -> Self {
1068 Self::new_with_station_atmosphere_resolution(
1069 inputs,
1070 wind,
1071 atmosphere,
1072 StationAtmosphereResolution::Authoritative,
1073 )
1074 }
1075
1076 fn new_with_station_atmosphere_resolution(
1077 mut inputs: BallisticInputs,
1078 wind: WindConditions,
1079 atmosphere: AtmosphericConditions,
1080 station_atmosphere_resolution: StationAtmosphereResolution,
1081 ) -> Self {
1082 inputs.normalize_for_solve();
1087
1088 let cluster_bc = if inputs.use_cluster_bc {
1090 Some(ClusterBCDegradation::new())
1091 } else {
1092 None
1093 };
1094 let precession_nutation_inertias = projectile_moments_of_inertia(
1095 inputs.bullet_mass,
1096 inputs.bullet_diameter,
1097 inputs.bullet_length,
1098 );
1099
1100 Self {
1101 inputs,
1102 wind,
1103 atmosphere,
1104 station_atmosphere_resolution,
1105 max_range: 1000.0,
1106 time_step: 0.001,
1107 max_trajectory_points: MAX_TRAJECTORY_POINTS,
1108 cluster_bc,
1109 precession_nutation_inertias,
1110 wind_sock: None,
1111 atmo_sock: None,
1112 }
1113 }
1114
1115 pub fn set_max_range(&mut self, range: f64) {
1116 self.max_range = range;
1117 }
1118
1119 pub fn set_time_step(&mut self, step: f64) {
1120 self.time_step = step;
1121 }
1122
1123 pub(crate) fn calculate_and_set_zero_angle(
1127 &mut self,
1128 target_distance_m: f64,
1129 target_height_m: f64,
1130 frame: ZeroTargetFrame,
1131 ) -> Result<f64, BallisticsError> {
1132 let angle = self.find_zero_angle(target_distance_m, target_height_m, frame)?;
1133 let angle = if target_distance_m > 0.0 {
1142 angle + self.inputs.zero_poi_vertical_m / target_distance_m
1143 } else {
1144 angle
1145 };
1146 self.inputs.muzzle_angle = angle;
1147 self.apply_windage_zero_bias(target_distance_m);
1148 Ok(angle)
1149 }
1150
1151 pub(crate) fn apply_windage_zero_bias(&mut self, target_distance_m: f64) {
1165 self.inputs.azimuth_angle += self.inputs.windage_zero_bias_rad(target_distance_m);
1166 }
1167
1168 fn find_zero_angle(
1169 &self,
1170 target_distance_m: f64,
1171 target_height_m: f64,
1172 frame: ZeroTargetFrame,
1173 ) -> Result<f64, BallisticsError> {
1174 let mut low_angle = 0.0;
1177 let mut high_angle = 0.2; let tolerance = 1e-7;
1179 let max_iterations = 60;
1180
1181 let low_height = self.zero_trial_height_at(low_angle, target_distance_m, frame)?;
1183 let high_height = self.zero_trial_height_at(high_angle, target_distance_m, frame)?;
1184
1185 match (low_height, high_height) {
1186 (Some(low_height), Some(high_height)) => {
1187 let low_error = low_height - target_height_m;
1188 let high_error = high_height - target_height_m;
1189
1190 if low_error > 0.0 && high_error > 0.0 {
1191 } else if low_error < 0.0 && high_error < 0.0 {
1194 let mut expanded = false;
1196 for multiplier in [2.0, 3.0, 4.0] {
1197 let new_high = (high_angle * multiplier).min(0.785);
1198 if let Ok(Some(height)) =
1199 self.zero_trial_height_at(new_high, target_distance_m, frame)
1200 {
1201 if height - target_height_m > 0.0 {
1202 high_angle = new_high;
1203 expanded = true;
1204 break;
1205 }
1206 }
1207 if new_high >= 0.785 {
1208 break;
1209 }
1210 }
1211 if !expanded {
1212 return Err("Cannot find zero angle: target beyond effective range even at maximum angle".into());
1213 }
1214 }
1215 }
1216 (None, Some(_)) => {
1217 }
1220 (Some(_), None) => {
1221 return Err(
1222 "Cannot find zero angle: high angle trajectory doesn't reach target distance"
1223 .into(),
1224 );
1225 }
1226 (None, None) => {
1227 return Err(
1228 "Cannot find zero angle: trajectory cannot reach target distance at any angle"
1229 .into(),
1230 );
1231 }
1232 }
1233
1234 for _ in 0..max_iterations {
1235 let mid_angle = (low_angle + high_angle) / 2.0;
1236 match self.zero_trial_height_at(mid_angle, target_distance_m, frame)? {
1237 Some(height) => {
1238 let error = height - target_height_m;
1239 if error.abs() < 0.0001 {
1242 return Ok(mid_angle);
1243 }
1244
1245 if (high_angle - low_angle).abs() < tolerance {
1248 if error.abs() < 0.01 {
1249 return Ok(mid_angle);
1250 }
1251 return Err("Zero angle did not converge: residual height error too large (target not reachable / not bracketed)".into());
1252 }
1253
1254 if error > 0.0 {
1255 high_angle = mid_angle;
1256 } else {
1257 low_angle = mid_angle;
1258 }
1259 }
1260 None => {
1261 low_angle = mid_angle;
1262 if (high_angle - low_angle).abs() < tolerance {
1263 return Err("Trajectory cannot reach target distance - angle converged without valid solution".into());
1264 }
1265 }
1266 }
1267 }
1268
1269 Err("Failed to find zero angle".into())
1270 }
1271
1272 fn zero_trial_height_at(
1275 &self,
1276 angle_rad: f64,
1277 target_distance_m: f64,
1278 frame: ZeroTargetFrame,
1279 ) -> Result<Option<f64>, BallisticsError> {
1280 let mut trial = self.clone();
1281 trial.inputs.muzzle_angle = angle_rad;
1282 trial.inputs.enable_aerodynamic_jump = false;
1285 trial.inputs.cant_angle = 0.0;
1288 if frame == ZeroTargetFrame::SightLine {
1295 trial.inputs.shooting_angle = 0.0;
1296 }
1297 trial.set_max_range(target_distance_m * 2.0);
1298 let result = trial.solve()?;
1299
1300 for (index, point) in result.points.iter().enumerate() {
1301 if point.position.x >= target_distance_m {
1302 let shot_y_m = if index == 0 {
1303 point.position.y
1304 } else {
1305 let previous = &result.points[index - 1];
1306 let span = point.position.x - previous.position.x;
1307 let fraction = (target_distance_m - previous.position.x) / span;
1308 previous.position.y + fraction * (point.position.y - previous.position.y)
1309 };
1310 return Ok(Some(crate::atmosphere::shot_frame_altitude(
1311 0.0,
1312 target_distance_m,
1313 shot_y_m,
1314 trial.inputs.shooting_angle,
1315 )));
1316 }
1317 }
1318 Ok(None)
1319 }
1320
1321 fn find_zero_range(
1348 &self,
1349 angle_rad: f64,
1350 target_height_m: f64,
1351 frame: ZeroTargetFrame,
1352 ) -> Result<ZeroCrossings, BallisticsError> {
1353 let mut trial = self.clone();
1354 trial.inputs.muzzle_angle = angle_rad;
1355 trial.inputs.enable_aerodynamic_jump = false;
1358 trial.inputs.cant_angle = 0.0;
1359 if frame == ZeroTargetFrame::SightLine {
1360 trial.inputs.shooting_angle = 0.0;
1361 }
1362 let result = trial.solve()?;
1363
1364 let mut near_crossing: Option<f64> = None;
1372 let mut far_crossing: Option<f64> = None;
1373 let mut previous: Option<(f64, f64)> = None; for point in &result.points {
1375 let height = crate::atmosphere::shot_frame_altitude(
1376 0.0,
1377 point.position.x,
1378 point.position.y,
1379 trial.inputs.shooting_angle,
1380 );
1381 let error = height - target_height_m;
1382 if let Some((prev_x, prev_error)) = previous {
1383 if prev_error == 0.0 {
1384 if near_crossing.is_none() {
1387 near_crossing = Some(prev_x);
1388 } else {
1389 far_crossing = Some(prev_x);
1390 }
1391 }
1392 if prev_error * error < 0.0 {
1393 let fraction = prev_error / (prev_error - error);
1394 let crossing = prev_x + fraction * (point.position.x - prev_x);
1395 if prev_error < 0.0 && error > 0.0 {
1396 if near_crossing.is_none() {
1398 near_crossing = Some(crossing);
1399 }
1400 } else {
1401 far_crossing = Some(crossing);
1403 }
1404 }
1405 }
1406 previous = Some((point.position.x, error));
1407 }
1408 if let Some((last_x, last_error)) = previous {
1411 if last_error == 0.0 {
1412 if near_crossing.is_none() {
1413 near_crossing = Some(last_x);
1414 } else {
1415 far_crossing = Some(last_x);
1416 }
1417 }
1418 }
1419
1420 if near_crossing.is_none() && far_crossing.is_none() {
1421 return Err(BallisticsError::from(
1422 "Cannot find zero range: this angle never crosses the target height within the \
1423 solved range (angle too shallow to reach it, or both crossings lie beyond \
1424 the solver's max range)."
1425 .to_string(),
1426 ));
1427 }
1428
1429 Ok(ZeroCrossings {
1430 near_m: near_crossing,
1431 far_m: far_crossing,
1432 })
1433 }
1434
1435 pub fn equivalent_horizontal_range(
1459 &self,
1460 target_range_m: f64,
1461 zero_distance_m: f64,
1462 ) -> Option<f64> {
1463 if !target_range_m.is_finite() || !zero_distance_m.is_finite() {
1464 return None;
1465 }
1466 if target_range_m <= zero_distance_m || target_range_m <= 0.0 {
1467 return None;
1468 }
1469
1470 fn path_y_at(points: &[TrajectoryPoint], distance_m: f64) -> Option<f64> {
1475 match bracket_param(points.len(), |i| points[i].position.x, distance_m) {
1476 Bracket::Below => Some(points[0].position.y),
1481 Bracket::Above | Bracket::Degenerate => None,
1482 Bracket::Inside { lo, t } => {
1483 let hi = lo + 1;
1484 Some(
1485 points[lo].position.y
1486 + t * (points[hi].position.y - points[lo].position.y),
1487 )
1488 }
1489 }
1490 }
1491
1492 let mut inclined = self.clone();
1496 inclined.inputs.enable_trajectory_sampling = false;
1497 let inclined_result = inclined.solve().ok()?;
1498 let los_height = inclined_result.line_of_sight_height_m;
1499 let inclined_drop = los_height - path_y_at(&inclined_result.points, target_range_m)?;
1500 let correction = inclined_drop / target_range_m;
1501 if correction <= 0.0 {
1502 return None;
1503 }
1504
1505 let mut flat = self.clone();
1507 flat.inputs.enable_trajectory_sampling = false;
1508 flat.inputs.shooting_angle = 0.0;
1509 let flat_result = flat.solve().ok()?;
1510 let flat_correction_at = |range_m: f64| -> Option<f64> {
1511 Some((los_height - path_y_at(&flat_result.points, range_m)?) / range_m)
1512 };
1513
1514 let flat_terminal_m = flat_result.points.last().map(|p| p.position.x)?;
1522 let mut low = zero_distance_m.max(1.0);
1523 let mut high = target_range_m.min(flat_terminal_m);
1524 if high <= low {
1525 return None;
1526 }
1527 if flat_correction_at(low)? - correction > 0.0 {
1528 return None; }
1530 if flat_correction_at(high)? - correction < 0.0 {
1531 return None; }
1533 for _ in 0..60 {
1534 let mid = 0.5 * (low + high);
1535 let error = flat_correction_at(mid)? - correction;
1536 if error.abs() == 0.0 {
1537 return Some(mid);
1538 }
1539 if error < 0.0 {
1540 low = mid;
1541 } else {
1542 high = mid;
1543 }
1544 if high - low < 0.01 {
1545 break;
1546 }
1547 }
1548 Some(0.5 * (low + high))
1549 }
1550
1551 fn validate_for_solve(&self) -> Result<(), BallisticsError> {
1557 let require_finite = |name: &str, value: f64| {
1558 if value.is_finite() {
1559 Ok(())
1560 } else {
1561 Err(BallisticsError::from(format!("{name} must be finite")))
1562 }
1563 };
1564 let require_positive = |name: &str, value: f64| {
1565 if value.is_finite() && value > 0.0 {
1566 Ok(())
1567 } else {
1568 Err(BallisticsError::from(format!(
1569 "{name} must be finite and greater than zero"
1570 )))
1571 }
1572 };
1573
1574 if self.inputs.custom_drag_table.is_none() {
1580 require_positive("bc_value", self.inputs.bc_value)?;
1581 }
1582 require_positive("bullet_mass", self.inputs.bullet_mass)?;
1583 require_positive("bullet_diameter", self.inputs.bullet_diameter)?;
1584 require_positive("muzzle_velocity", self.inputs.muzzle_velocity)?;
1585 require_positive("cd_scale", self.inputs.cd_scale)?;
1591
1592 require_finite("muzzle_angle", self.inputs.muzzle_angle)?;
1593 require_finite("azimuth_angle", self.inputs.azimuth_angle)?;
1594 require_finite("shooting_angle", self.inputs.shooting_angle)?;
1595 require_finite("cant_angle", self.inputs.cant_angle)?;
1596 require_finite("muzzle_height", self.inputs.muzzle_height)?;
1597
1598 for (name, value) in [
1602 ("zero_poi_vertical_m", self.inputs.zero_poi_vertical_m),
1603 ("zero_poi_horizontal_m", self.inputs.zero_poi_horizontal_m),
1604 ] {
1605 require_finite(name, value)?;
1606 if value.abs() >= 1.0 {
1607 return Err(BallisticsError::from(format!(
1608 "{name} must be smaller than 1.0 m in magnitude (it is a linear POI \
1609 offset at the zero range, in meters)"
1610 )));
1611 }
1612 }
1613
1614 require_finite(
1618 "sight_offset_lateral_m",
1619 self.inputs.sight_offset_lateral_m,
1620 )?;
1621 if self.inputs.sight_offset_lateral_m.abs() >= 0.5 {
1622 return Err(BallisticsError::from(
1623 "sight_offset_lateral_m must be smaller than 0.5 m in magnitude (it is \
1624 the lateral sight-to-bore mount offset, in meters)",
1625 ));
1626 }
1627
1628 if !(self.inputs.ground_threshold.is_finite()
1631 || self.inputs.ground_threshold == f64::NEG_INFINITY)
1632 {
1633 return Err(BallisticsError::from(
1634 "ground_threshold must be finite or negative infinity",
1635 ));
1636 }
1637
1638 match &self.wind_sock {
1639 Some(wind_sock) => wind_sock
1640 .validate_segments()
1641 .map_err(BallisticsError::from)?,
1642 None => {
1643 require_finite("wind.speed", self.wind.speed)?;
1644 require_finite("wind.direction", self.wind.direction)?;
1645 require_finite("wind.vertical_speed", self.wind.vertical_speed)?;
1646 }
1647 }
1648
1649 require_finite("atmosphere.temperature", self.atmosphere.temperature)?;
1650 require_finite("atmosphere.pressure", self.atmosphere.pressure)?;
1651 require_finite("atmosphere.humidity", self.atmosphere.humidity)?;
1652 require_finite("atmosphere.altitude", self.atmosphere.altitude)?;
1653
1654 require_positive("max_range", self.max_range)?;
1655 if !self.inputs.use_rk4 || !self.inputs.use_adaptive_rk45 {
1658 require_positive("time_step", self.time_step)?;
1659 }
1660
1661 if self.inputs.enable_trajectory_sampling {
1662 require_finite("sight_height", self.inputs.sight_height)?;
1663 require_positive("sample_interval", self.inputs.sample_interval)?;
1664 projected_sample_count(self.max_range, self.inputs.sample_interval)?;
1665 }
1666
1667 if self.inputs.drops_reference == DropsReference::Target {
1671 require_finite("target_height", self.inputs.target_height)?;
1672 if self.inputs.shooting_angle.cos() <= 1e-9 {
1673 return Err(BallisticsError::from(
1674 "drops reference 'target' is undefined for shooting angles at or beyond 90 degrees",
1675 ));
1676 }
1677 }
1678
1679 if self.inputs.enable_coriolis {
1680 require_finite("shot_azimuth", self.inputs.shot_azimuth)?;
1681 if let Some(latitude) = self.inputs.latitude {
1682 require_finite("latitude", latitude)?;
1683 }
1684 }
1685
1686 Ok(())
1687 }
1688
1689 fn validate_result_sanity(&self, result: &TrajectoryResult) -> Result<(), BallisticsError> {
1696 let require_finite = |name: &str, value: f64| {
1697 if value.is_finite() {
1698 Ok(())
1699 } else {
1700 Err(BallisticsError::from(format!(
1701 "trajectory result contains non-finite {name}"
1702 )))
1703 }
1704 };
1705 let require_non_negative = |name: &str, value: f64| {
1706 if value >= 0.0 {
1707 Ok(())
1708 } else {
1709 Err(BallisticsError::from(format!(
1710 "trajectory result contains non-physical negative {name} ({value})"
1711 )))
1712 }
1713 };
1714 let require_indexed_finite = |collection: &str, index: usize, field: &str, value: f64| {
1715 if value.is_finite() {
1716 Ok(())
1717 } else {
1718 Err(BallisticsError::from(format!(
1719 "trajectory result contains non-finite {collection}[{index}].{field}"
1720 )))
1721 }
1722 };
1723 let require_indexed_non_negative =
1724 |collection: &str, index: usize, field: &str, value: f64| {
1725 if value >= 0.0 {
1726 Ok(())
1727 } else {
1728 Err(BallisticsError::from(format!(
1729 "trajectory result contains non-physical negative {collection}[{index}].{field} ({value})"
1730 )))
1731 }
1732 };
1733
1734 require_finite("max_range", result.max_range)?;
1735 require_finite("max_height", result.max_height)?;
1736 require_finite("time_of_flight", result.time_of_flight)?;
1737 require_finite("impact_velocity", result.impact_velocity)?;
1738 require_finite("impact_energy", result.impact_energy)?;
1739 require_finite("projectile_mass_kg", result.projectile_mass_kg)?;
1740 require_finite(
1741 "line_of_sight_height_m",
1742 result.line_of_sight_height_m,
1743 )?;
1744 require_finite(
1745 "station_speed_of_sound_mps",
1746 result.station_speed_of_sound_mps,
1747 )?;
1748
1749 require_non_negative("max_range", result.max_range)?;
1753 require_non_negative("time_of_flight", result.time_of_flight)?;
1754 require_non_negative("impact_velocity", result.impact_velocity)?;
1755 require_non_negative("impact_energy", result.impact_energy)?;
1756 require_non_negative("projectile_mass_kg", result.projectile_mass_kg)?;
1757 require_non_negative(
1758 "station_speed_of_sound_mps",
1759 result.station_speed_of_sound_mps,
1760 )?;
1761
1762 for (index, point) in result.points.iter().enumerate() {
1763 require_indexed_finite("points", index, "time", point.time)?;
1764 require_indexed_finite("points", index, "position.x", point.position.x)?;
1765 require_indexed_finite("points", index, "position.y", point.position.y)?;
1766 require_indexed_finite("points", index, "position.z", point.position.z)?;
1767 require_indexed_finite(
1768 "points",
1769 index,
1770 "velocity_magnitude",
1771 point.velocity_magnitude,
1772 )?;
1773 require_indexed_finite("points", index, "kinetic_energy", point.kinetic_energy)?;
1774 require_indexed_non_negative("points", index, "time", point.time)?;
1775 require_indexed_non_negative(
1776 "points",
1777 index,
1778 "velocity_magnitude",
1779 point.velocity_magnitude,
1780 )?;
1781 require_indexed_non_negative("points", index, "kinetic_energy", point.kinetic_energy)?;
1782 }
1783
1784 if let Some(samples) = &result.sampled_points {
1785 for (index, sample) in samples.iter().enumerate() {
1786 require_indexed_finite("sampled_points", index, "distance_m", sample.distance_m)?;
1787 require_indexed_finite("sampled_points", index, "drop_m", sample.drop_m)?;
1788 require_indexed_finite(
1789 "sampled_points",
1790 index,
1791 "wind_drift_m",
1792 sample.wind_drift_m,
1793 )?;
1794 require_indexed_finite(
1795 "sampled_points",
1796 index,
1797 "velocity_mps",
1798 sample.velocity_mps,
1799 )?;
1800 require_indexed_finite("sampled_points", index, "energy_j", sample.energy_j)?;
1801 require_indexed_finite("sampled_points", index, "time_s", sample.time_s)?;
1802 }
1803 }
1804
1805 for (name, value) in [
1806 ("min_pitch_damping", result.min_pitch_damping),
1807 ("transonic_mach", result.transonic_mach),
1808 ("max_yaw_angle", result.max_yaw_angle),
1809 ("max_precession_angle", result.max_precession_angle),
1810 ] {
1811 if let Some(value) = value {
1812 require_finite(name, value)?;
1813 }
1814 }
1815
1816 if let Some(state) = result.angular_state {
1817 for (name, value) in [
1818 ("angular_state.pitch_angle", state.pitch_angle),
1819 ("angular_state.yaw_angle", state.yaw_angle),
1820 ("angular_state.pitch_rate", state.pitch_rate),
1821 ("angular_state.yaw_rate", state.yaw_rate),
1822 ("angular_state.precession_angle", state.precession_angle),
1823 ("angular_state.nutation_phase", state.nutation_phase),
1824 ] {
1825 require_finite(name, value)?;
1826 }
1827 }
1828
1829 if let Some(jump) = result.aerodynamic_jump {
1830 for (name, value) in [
1831 ("aerodynamic_jump.vertical_jump_moa", jump.vertical_jump_moa),
1832 (
1833 "aerodynamic_jump.horizontal_jump_moa",
1834 jump.horizontal_jump_moa,
1835 ),
1836 ("aerodynamic_jump.jump_angle_rad", jump.jump_angle_rad),
1837 (
1838 "aerodynamic_jump.magnus_component_moa",
1839 jump.magnus_component_moa,
1840 ),
1841 ("aerodynamic_jump.yaw_component_moa", jump.yaw_component_moa),
1842 (
1843 "aerodynamic_jump.stabilization_factor",
1844 jump.stabilization_factor,
1845 ),
1846 ] {
1847 require_finite(name, value)?;
1848 }
1849 }
1850
1851 Ok(())
1852 }
1853
1854 fn validate_integration_state(
1867 &self,
1868 position: &Vector3<f64>,
1869 velocity: &Vector3<f64>,
1870 time: f64,
1871 ) -> Result<(), BallisticsError> {
1872 if !(position.iter().all(|value| value.is_finite())
1873 && velocity.iter().all(|value| value.is_finite())
1874 && time.is_finite())
1875 {
1876 return Err(BallisticsError::from(
1877 "trajectory integration produced a non-finite state (often from physically \
1878 extreme inputs — e.g. an absurd bore/muzzle height placing the launch far \
1879 from sea level, or a degenerate atmosphere; check those inputs, or set \
1880 --altitude explicitly)",
1881 ));
1882 }
1883
1884 let speed = velocity.magnitude();
1885 let budget = self.speed_budget(time);
1886 if speed > budget {
1887 return Err(BallisticsError::from(format!(
1888 "trajectory integration diverged: speed {speed:.3e} m/s at t={time:.6}s exceeds \
1889 the physical budget of {budget:.3e} m/s"
1890 )));
1891 }
1892 Ok(())
1893 }
1894
1895 fn speed_budget(&self, time: f64) -> f64 {
1900 let scalar_wind = self.wind.speed.abs() + self.wind.vertical_speed.abs();
1901 let wind_bound = match &self.wind_sock {
1902 Some(sock) => scalar_wind.max(sock.max_speed_mps()),
1903 None => scalar_wind,
1904 };
1905 2.0 * (self.inputs.muzzle_velocity + wind_bound + 10.0)
1906 + crate::constants::G_ACCEL_MPS2 * time
1907 }
1908
1909 fn push_trajectory_point(
1911 &self,
1912 points: &mut Vec<TrajectoryPoint>,
1913 point: TrajectoryPoint,
1914 ) -> Result<(), BallisticsError> {
1915 if points.len() >= self.max_trajectory_points {
1916 return Err(BallisticsError::from(format!(
1917 "trajectory point limit of {} exceeded",
1918 self.max_trajectory_points
1919 )));
1920 }
1921 points.push(point);
1922 Ok(())
1923 }
1924
1925 pub fn set_wind_segments(&mut self, segments: Vec<crate::wind::WindSegment>) {
1932 self.wind_sock = if segments.is_empty() {
1933 None
1934 } else {
1935 Some(crate::wind::WindSock::new(segments))
1936 };
1937 }
1938
1939 pub fn set_atmo_segments(&mut self, segments: Vec<crate::atmosphere::AtmoSegment>) {
1948 self.atmo_sock = if segments.is_empty() {
1949 None
1950 } else {
1951 Some(crate::atmosphere::AtmoSock::new(segments))
1952 };
1953 }
1954
1955 fn launch_angles_from(
1964 &self,
1965 aj: Option<&crate::aerodynamic_jump::AerodynamicJumpComponents>,
1966 ) -> (f64, f64) {
1967 let (mut elev, mut azim) = (self.inputs.muzzle_angle, self.inputs.azimuth_angle);
1968 if self.inputs.cant_angle != 0.0 {
1975 let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
1976 let (e0, a0) = (elev, azim);
1977 elev = e0 * cos_c - a0 * sin_c;
1978 azim = a0 * cos_c + e0 * sin_c;
1979 }
1980 match aj {
1981 Some(c) => {
1982 const MOA_PER_RAD: f64 = 3437.7467707849;
1984 (
1985 elev + c.vertical_jump_moa / MOA_PER_RAD,
1986 azim + c.horizontal_jump_moa / MOA_PER_RAD,
1987 )
1988 }
1989 None => (elev, azim),
1990 }
1991 }
1992
1993 fn aerodynamic_jump_components(
2001 &self,
2002 ) -> Option<crate::aerodynamic_jump::AerodynamicJumpComponents> {
2003 if !self.inputs.enable_aerodynamic_jump {
2004 return None;
2005 }
2006 let diameter_m = self.inputs.bullet_diameter;
2010 if !(self.inputs.twist_rate.is_finite()
2011 && self.inputs.twist_rate != 0.0
2012 && diameter_m.is_finite()
2013 && diameter_m > 0.0
2014 && self.inputs.bullet_length.is_finite()
2015 && self.inputs.bullet_length > 0.0
2016 && self.inputs.muzzle_velocity.is_finite())
2017 {
2018 return None;
2019 }
2020
2021 let (_, _, temp_c, pressure_hpa) = self.resolved_atmosphere();
2023 let sg = crate::stability::compute_stability_coefficient(
2024 &self.inputs,
2025 (self.atmosphere.altitude, temp_c, pressure_hpa, 0.0),
2026 );
2027 if !(sg.is_finite() && sg > 0.0) {
2028 return None;
2029 }
2030 let length_calibers = self.inputs.bullet_length / diameter_m;
2031
2032 const MS_TO_MPH: f64 = 2.236_936_292_054_4;
2038 let crosswind_from_right_mps = if let Some(sock) = &self.wind_sock {
2039 -sock.vector_for_range_stateless(0.0)[2]
2040 } else {
2041 self.wind.speed * self.wind.direction.sin()
2042 };
2043 let crosswind_from_right_mph = crosswind_from_right_mps * MS_TO_MPH;
2044
2045 let vertical_jump_moa = crate::aerodynamic_jump::litz_crosswind_jump_moa(
2046 sg,
2047 length_calibers,
2048 crosswind_from_right_mph,
2049 self.inputs.is_twist_right,
2050 );
2051 if !vertical_jump_moa.is_finite() {
2052 return None;
2053 }
2054
2055 const MOA_PER_RAD: f64 = 3437.7467707849;
2056 Some(crate::aerodynamic_jump::AerodynamicJumpComponents {
2057 vertical_jump_moa,
2058 horizontal_jump_moa: 0.0,
2060 jump_angle_rad: vertical_jump_moa.abs() / MOA_PER_RAD,
2061 magnus_component_moa: 0.0,
2062 yaw_component_moa: 0.0,
2063 stabilization_factor: (sg / 1.5).clamp(0.0, 1.0),
2064 })
2065 }
2066
2067 fn resolved_atmosphere(&self) -> (f64, f64, f64, f64) {
2068 let (temp_c, pressure_hpa) = match self.station_atmosphere_resolution {
2069 StationAtmosphereResolution::LegacyDefaultSentinels => {
2070 crate::atmosphere::resolve_station_conditions(
2071 self.atmosphere.temperature,
2072 self.atmosphere.pressure,
2073 self.atmosphere.altitude,
2074 )
2075 }
2076 StationAtmosphereResolution::Authoritative => {
2077 (self.atmosphere.temperature, self.atmosphere.pressure)
2078 }
2079 };
2080 let (density, speed_of_sound) = crate::atmosphere::calculate_atmosphere(
2081 self.atmosphere.altitude,
2082 Some(temp_c),
2083 Some(pressure_hpa),
2084 self.atmosphere.humidity,
2085 );
2086 (density, speed_of_sound, temp_c, pressure_hpa)
2087 }
2088
2089 fn precession_nutation_params(
2090 &self,
2091 velocity_mps: f64,
2092 air_density_kg_m3: f64,
2093 speed_of_sound_mps: f64,
2094 ) -> PrecessionNutationParams {
2095 let (spin_inertia, transverse_inertia) = self.precession_nutation_inertias;
2096 let spin_rate_rad_s = if self.inputs.twist_rate > 0.0 {
2097 let velocity_fps = velocity_mps * 3.28084;
2098 let twist_rate_ft = self.inputs.twist_rate / 12.0;
2099 (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI
2100 } else {
2101 0.0
2102 };
2103
2104 PrecessionNutationParams {
2105 mass_kg: self.inputs.bullet_mass,
2106 caliber_m: self.inputs.bullet_diameter,
2107 length_m: self.inputs.bullet_length,
2108 spin_rate_rad_s,
2109 spin_inertia,
2110 transverse_inertia,
2111 velocity_mps,
2112 air_density_kg_m3,
2113 mach: velocity_mps / speed_of_sound_mps,
2114 pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
2115 nutation_damping_factor: 0.05,
2116 }
2117 }
2118
2119 fn append_terminal_endpoint(
2126 &self,
2127 points: &mut Vec<TrajectoryPoint>,
2128 post_position: Vector3<f64>,
2129 post_velocity: Vector3<f64>,
2130 post_time: f64,
2131 max_height: &mut f64,
2132 ) -> Result<TrajectoryTermination, BallisticsError> {
2133 let previous = points
2134 .last()
2135 .cloned()
2136 .ok_or_else(|| BallisticsError::from("No trajectory points generated"))?;
2137
2138 let mut crossings = Vec::with_capacity(3);
2139 if previous.position.x < self.max_range && post_position.x >= self.max_range {
2140 let span = post_position.x - previous.position.x;
2141 if span.is_finite() && span > 0.0 {
2142 crossings.push((
2143 (self.max_range - previous.position.x) / span,
2144 TrajectoryTermination::MaxRange,
2145 ));
2146 }
2147 }
2148 if self.inputs.ground_threshold.is_finite()
2149 && previous.position.y > self.inputs.ground_threshold
2150 && post_position.y <= self.inputs.ground_threshold
2151 {
2152 let span = post_position.y - previous.position.y;
2153 if span.is_finite() && span < 0.0 {
2154 crossings.push((
2155 (self.inputs.ground_threshold - previous.position.y) / span,
2156 TrajectoryTermination::GroundThreshold,
2157 ));
2158 }
2159 }
2160 if previous.time < TRAJECTORY_TIME_LIMIT_S && post_time >= TRAJECTORY_TIME_LIMIT_S {
2161 let span = post_time - previous.time;
2162 if span.is_finite() && span > 0.0 {
2163 crossings.push((
2164 (TRAJECTORY_TIME_LIMIT_S - previous.time) / span,
2165 TrajectoryTermination::TimeLimit,
2166 ));
2167 }
2168 }
2169
2170 let (fraction, termination) = crossings
2171 .into_iter()
2172 .filter(|(fraction, _)| fraction.is_finite() && (0.0..=1.0).contains(fraction))
2173 .min_by(|left, right| {
2174 let priority = |termination: TrajectoryTermination| match termination {
2175 TrajectoryTermination::GroundThreshold => 0,
2176 TrajectoryTermination::MaxRange => 1,
2177 TrajectoryTermination::TimeLimit => 2,
2178 TrajectoryTermination::VelocityFloor => 3,
2179 };
2180 left.0
2181 .total_cmp(&right.0)
2182 .then_with(|| priority(left.1).cmp(&priority(right.1)))
2183 })
2184 .ok_or_else(|| {
2185 BallisticsError::from(
2186 "trajectory integration stopped without crossing a supported boundary",
2187 )
2188 })?;
2189
2190 let mut position = previous.position + (post_position - previous.position) * fraction;
2191 match termination {
2192 TrajectoryTermination::MaxRange => position.x = self.max_range,
2193 TrajectoryTermination::GroundThreshold => {
2194 position.y = self.inputs.ground_threshold;
2195 }
2196 TrajectoryTermination::TimeLimit | TrajectoryTermination::VelocityFloor => {}
2197 }
2198 let velocity_magnitude = previous.velocity_magnitude
2199 + (post_velocity.magnitude() - previous.velocity_magnitude) * fraction;
2200 let mut time = previous.time + (post_time - previous.time) * fraction;
2201 if termination == TrajectoryTermination::TimeLimit {
2202 time = TRAJECTORY_TIME_LIMIT_S;
2203 }
2204 let kinetic_energy =
2205 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2206
2207 if position.y > *max_height {
2208 *max_height = position.y;
2209 }
2210 let terminal_point = TrajectoryPoint {
2211 time,
2212 position,
2213 velocity_magnitude,
2214 kinetic_energy,
2215 drag_coefficient: None,
2216 };
2217 if terminal_point.position.x < previous.position.x {
2218 return Err(BallisticsError::from(
2219 "trajectory terminal state reversed downrange before the crossed boundary",
2220 ));
2221 }
2222 if terminal_point.position.x == previous.position.x {
2223 let last = points.last_mut().ok_or_else(|| {
2228 BallisticsError::from("trajectory points disappeared during terminal finalization")
2229 })?;
2230 *last = terminal_point;
2231 } else {
2232 self.push_trajectory_point(points, terminal_point)?;
2233 }
2234 Ok(termination)
2235 }
2236
2237 fn gravity_acceleration(&self) -> Vector3<f64> {
2238 let theta = self.inputs.shooting_angle;
2239 Vector3::new(
2240 -crate::constants::G_ACCEL_MPS2 * theta.sin(),
2241 -crate::constants::G_ACCEL_MPS2 * theta.cos(),
2242 0.0,
2243 )
2244 }
2245
2246 fn get_wind_at_altitude(&self, altitude_m: f64) -> Vector3<f64> {
2247 let model = match self.inputs.wind_shear_model.as_str() {
2262 "logarithmic" => WindShearModel::Logarithmic,
2263 "power_law" | "powerlaw" | "exponential" => WindShearModel::PowerLaw,
2264 "ekman_spiral" | "ekman" => WindShearModel::EkmanSpiral,
2265 "custom_layers" | "custom" => WindShearModel::CustomLayers,
2266 _ => WindShearModel::PowerLaw,
2267 };
2268 let speed_ratio = crate::wind_shear::boundary_layer_speed_ratio(altitude_m, model);
2269
2270 crate::wind::wind_vector(self.wind.speed, self.wind.direction, 0.0) * speed_ratio
2277 + Vector3::new(0.0, self.wind.vertical_speed, 0.0)
2278 }
2279
2280 pub fn solve(&self) -> Result<TrajectoryResult, BallisticsError> {
2281 self.validate_for_solve()?;
2282 let mut result = if self.inputs.use_rk4 {
2283 if self.inputs.use_adaptive_rk45 {
2284 self.solve_rk45()?
2285 } else {
2286 self.solve_rk4()?
2287 }
2288 } else {
2289 self.solve_euler()?
2290 };
2291 self.apply_spin_drift(&mut result);
2292 self.validate_result_sanity(&result)?;
2293 Ok(result)
2294 }
2295
2296 fn apply_spin_drift(&self, result: &mut TrajectoryResult) {
2302 if !self.inputs.use_enhanced_spin_drift {
2303 return;
2304 }
2305 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 {
2309 return;
2310 }
2311
2312 let sg = self.effective_spin_drift_sg();
2319
2320 for p in result.points.iter_mut() {
2321 if p.time <= 0.0 {
2322 continue;
2323 }
2324 p.position.z +=
2326 crate::spin_drift::litz_drift_meters(sg, p.time, self.inputs.is_twist_right);
2327 }
2328
2329 if let Some(samples) = result.sampled_points.as_mut() {
2333 for s in samples.iter_mut() {
2334 if s.time_s <= 0.0 {
2335 continue;
2336 }
2337 s.wind_drift_m +=
2338 crate::spin_drift::litz_drift_meters(sg, s.time_s, self.inputs.is_twist_right);
2339 }
2340 }
2341 }
2342
2343 fn effective_spin_drift_sg(&self) -> f64 {
2348 let (_, _, temp_c, press_hpa) = self.resolved_atmosphere();
2349 crate::spin_drift::effective_sg_from_inputs(&self.inputs, temp_c, press_hpa)
2350 }
2351
2352 fn initial_position(&self) -> Vector3<f64> {
2370 if self.inputs.cant_angle == 0.0 && self.inputs.sight_offset_lateral_m == 0.0 {
2371 return Vector3::new(0.0, self.inputs.muzzle_height, 0.0);
2372 }
2373 let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
2374 let sh = self.inputs.sight_height;
2375 let off = self.inputs.sight_offset_lateral_m;
2376 Vector3::new(
2377 0.0,
2378 self.inputs.muzzle_height + sh * (1.0 - cos_c) + off * sin_c,
2379 -sh * sin_c - off * cos_c,
2380 )
2381 }
2382
2383 fn build_sampled_points(
2390 &self,
2391 points: &[TrajectoryPoint],
2392 max_height: f64,
2393 transonic_distances: Vec<f64>,
2394 mach_transitions: &MachTransitionTracker,
2395 ) -> Result<Option<Vec<TrajectorySample>>, BallisticsError> {
2396 if !self.inputs.enable_trajectory_sampling {
2397 return Ok(None);
2398 }
2399
2400 let last_point = points.last().ok_or("No trajectory points generated")?;
2401 let trajectory_data = TrajectoryData {
2402 times: points.iter().map(|p| p.time).collect(),
2403 positions: points.iter().map(|p| p.position).collect(),
2404 velocities: points
2405 .iter()
2406 .map(|p| {
2407 Vector3::new(0.0, 0.0, p.velocity_magnitude)
2409 })
2410 .collect(),
2411 transonic_distances, mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2413 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2414 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2415 };
2416
2417 let sight_position_m = self.inputs.muzzle_height + self.inputs.sight_height;
2422 let target_reference = self.inputs.drops_reference == DropsReference::Target;
2423 let target_vertical_height_m = if target_reference && self.inputs.target_height != 0.0 {
2430 self.inputs.target_height
2431 } else {
2432 sight_position_m
2433 };
2434 let outputs = TrajectoryOutputs {
2435 target_distance_horiz_m: last_point.position.x, target_vertical_height_m,
2437 time_of_flight_s: last_point.time,
2438 max_ord_dist_horiz_m: max_height,
2439 sight_height_m: sight_position_m,
2440 };
2441
2442 let mut samples = sample_trajectory(
2444 &trajectory_data,
2445 &outputs,
2446 self.inputs.sample_interval,
2447 self.inputs.bullet_mass,
2448 )?;
2449
2450 if target_reference {
2457 let cos_theta = self.inputs.shooting_angle.cos();
2458 for sample in &mut samples {
2459 sample.drop_m /= cos_theta;
2460 }
2461 }
2462 Ok(Some(samples))
2463 }
2464
2465 fn solve_euler(&self) -> Result<TrajectoryResult, BallisticsError> {
2466 let mut time = 0.0;
2468 let mut position = self.initial_position();
2472 let aj_components = self.aerodynamic_jump_components();
2478 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2479 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2480 let mut velocity = Vector3::new(
2481 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2485
2486 let mut points = Vec::new();
2487 let mut max_height = position.y;
2488 let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
2494 let mut mach_transitions = MachTransitionTracker::default();
2495
2496 let mut angular_state = if self.inputs.enable_precession_nutation {
2498 Some(AngularState {
2499 pitch_angle: 0.001, yaw_angle: 0.001,
2501 pitch_rate: 0.0,
2502 yaw_rate: 0.0,
2503 precession_angle: 0.0,
2504 nutation_phase: 0.0,
2505 })
2506 } else {
2507 None
2508 };
2509 let mut max_yaw_angle = 0.0;
2510 let mut max_precession_angle = 0.0;
2511
2512 let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2514 self.resolved_atmosphere();
2515 let base_ratio = air_density / 1.225;
2520
2521 let wind_vector =
2527 crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
2528
2529 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2532 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2533 );
2534
2535 while position.x < self.max_range
2537 && position.y > self.inputs.ground_threshold
2538 && time < TRAJECTORY_TIME_LIMIT_S
2539 {
2540 let velocity_magnitude = velocity.magnitude();
2542 let kinetic_energy =
2543 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2544
2545 self.push_trajectory_point(
2546 &mut points,
2547 TrajectoryPoint {
2548 time,
2549 position,
2550 velocity_magnitude,
2551 kinetic_energy,
2552 drag_coefficient: None,
2553 },
2554 )?;
2555
2556 {
2559 let mach_here = if speed_of_sound > 0.0 {
2560 velocity_magnitude / speed_of_sound
2561 } else {
2562 0.0
2563 };
2564 mach_transitions.record_downward_crossings(
2565 mach_here,
2566 position.x,
2567 &mut transonic_distances,
2568 );
2569 }
2570
2571 if position.y > max_height {
2573 max_height = position.y;
2574 }
2575
2576 if self.inputs.enable_pitch_damping {
2578 let mach = velocity_magnitude / speed_of_sound;
2579
2580 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
2582 transonic_mach = Some(mach);
2583 }
2584
2585 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
2587
2588 if pitch_damping < min_pitch_damping {
2590 min_pitch_damping = pitch_damping;
2591 }
2592 }
2593
2594 if self.inputs.enable_precession_nutation {
2596 if let Some(ref mut state) = angular_state {
2597 let velocity_magnitude = velocity.magnitude();
2598 let params = self.precession_nutation_params(
2599 velocity_magnitude,
2600 air_density,
2601 speed_of_sound,
2602 );
2603
2604 *state = calculate_combined_angular_motion(
2606 ¶ms,
2607 state,
2608 time,
2609 self.time_step,
2610 0.001, );
2612
2613 if state.yaw_angle.abs() > max_yaw_angle {
2615 max_yaw_angle = state.yaw_angle.abs();
2616 }
2617 if state.precession_angle.abs() > max_precession_angle {
2618 max_precession_angle = state.precession_angle.abs();
2619 }
2620 }
2621 }
2622
2623 let acceleration = self.calculate_acceleration(
2630 &position,
2631 &velocity,
2632 &wind_vector,
2633 (resolved_temp_c, resolved_press_hpa, base_ratio),
2634 );
2635
2636 velocity += acceleration * self.time_step;
2638 position += velocity * self.time_step;
2639 time += self.time_step;
2640 self.validate_integration_state(&position, &velocity, time)?;
2641 }
2642
2643 let termination =
2644 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
2645
2646 self.annotate_drag_coefficients(&mut points, speed_of_sound);
2651
2652 let last_point = points.last().ok_or("No trajectory points generated")?;
2653
2654 let sampled_points = self.build_sampled_points(
2656 &points,
2657 max_height,
2658 transonic_distances,
2659 &mach_transitions,
2660 )?;
2661
2662 Ok(TrajectoryResult {
2663 max_range: last_point.position.x, max_height,
2665 time_of_flight: last_point.time,
2666 impact_velocity: last_point.velocity_magnitude,
2667 impact_energy: last_point.kinetic_energy,
2668 projectile_mass_kg: self.inputs.bullet_mass,
2669 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
2670 station_speed_of_sound_mps: speed_of_sound,
2671 termination,
2672 points,
2673 sampled_points,
2674 min_pitch_damping: if self.inputs.enable_pitch_damping {
2675 Some(min_pitch_damping)
2676 } else {
2677 None
2678 },
2679 transonic_mach,
2680 angular_state,
2681 max_yaw_angle: if self.inputs.enable_precession_nutation {
2682 Some(max_yaw_angle)
2683 } else {
2684 None
2685 },
2686 max_precession_angle: if self.inputs.enable_precession_nutation {
2687 Some(max_precession_angle)
2688 } else {
2689 None
2690 },
2691 aerodynamic_jump: aj_components,
2692 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2693 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2694 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2695 })
2696 }
2697
2698 fn solve_rk4(&self) -> Result<TrajectoryResult, BallisticsError> {
2699 let mut time = 0.0;
2701 let mut position = self.initial_position();
2706
2707 let aj_components = self.aerodynamic_jump_components();
2713 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2714 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2715 let mut velocity = Vector3::new(
2716 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2720
2721 let mut points = Vec::new();
2722 let mut max_height = position.y;
2723 let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
2729 let mut mach_transitions = MachTransitionTracker::default();
2730
2731 let mut angular_state = if self.inputs.enable_precession_nutation {
2733 Some(AngularState {
2734 pitch_angle: 0.001, yaw_angle: 0.001,
2736 pitch_rate: 0.0,
2737 yaw_rate: 0.0,
2738 precession_angle: 0.0,
2739 nutation_phase: 0.0,
2740 })
2741 } else {
2742 None
2743 };
2744 let mut max_yaw_angle = 0.0;
2745 let mut max_precession_angle = 0.0;
2746
2747 let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2749 self.resolved_atmosphere();
2750 let base_ratio = air_density / 1.225;
2755
2756 let wind_vector =
2762 crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
2763
2764 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2767 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2768 );
2769
2770 while position.x < self.max_range
2772 && position.y > self.inputs.ground_threshold
2773 && time < TRAJECTORY_TIME_LIMIT_S
2774 {
2775 let velocity_magnitude = velocity.magnitude();
2777 let kinetic_energy =
2778 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2779
2780 self.push_trajectory_point(
2781 &mut points,
2782 TrajectoryPoint {
2783 time,
2784 position,
2785 velocity_magnitude,
2786 kinetic_energy,
2787 drag_coefficient: None,
2788 },
2789 )?;
2790
2791 {
2794 let mach_here = if speed_of_sound > 0.0 {
2795 velocity_magnitude / speed_of_sound
2796 } else {
2797 0.0
2798 };
2799 mach_transitions.record_downward_crossings(
2800 mach_here,
2801 position.x,
2802 &mut transonic_distances,
2803 );
2804 }
2805
2806 if position.y > max_height {
2807 max_height = position.y;
2808 }
2809
2810 if self.inputs.enable_pitch_damping {
2812 let mach = velocity_magnitude / speed_of_sound;
2813
2814 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
2816 transonic_mach = Some(mach);
2817 }
2818
2819 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
2821
2822 if pitch_damping < min_pitch_damping {
2824 min_pitch_damping = pitch_damping;
2825 }
2826 }
2827
2828 if self.inputs.enable_precession_nutation {
2830 if let Some(ref mut state) = angular_state {
2831 let velocity_magnitude = velocity.magnitude();
2832 let params = self.precession_nutation_params(
2833 velocity_magnitude,
2834 air_density,
2835 speed_of_sound,
2836 );
2837
2838 *state = calculate_combined_angular_motion(
2840 ¶ms,
2841 state,
2842 time,
2843 self.time_step,
2844 0.001, );
2846
2847 if state.yaw_angle.abs() > max_yaw_angle {
2849 max_yaw_angle = state.yaw_angle.abs();
2850 }
2851 if state.precession_angle.abs() > max_precession_angle {
2852 max_precession_angle = state.precession_angle.abs();
2853 }
2854 }
2855 }
2856
2857 let dt = self.time_step;
2859
2860 let acc1 = self.calculate_acceleration(
2862 &position,
2863 &velocity,
2864 &wind_vector,
2865 (resolved_temp_c, resolved_press_hpa, base_ratio),
2866 );
2867
2868 let pos2 = position + velocity * (dt * 0.5);
2870 let vel2 = velocity + acc1 * (dt * 0.5);
2871 let acc2 = self.calculate_acceleration(
2872 &pos2,
2873 &vel2,
2874 &wind_vector,
2875 (resolved_temp_c, resolved_press_hpa, base_ratio),
2876 );
2877
2878 let pos3 = position + vel2 * (dt * 0.5);
2880 let vel3 = velocity + acc2 * (dt * 0.5);
2881 let acc3 = self.calculate_acceleration(
2882 &pos3,
2883 &vel3,
2884 &wind_vector,
2885 (resolved_temp_c, resolved_press_hpa, base_ratio),
2886 );
2887
2888 let pos4 = position + vel3 * dt;
2890 let vel4 = velocity + acc3 * dt;
2891 let acc4 = self.calculate_acceleration(
2892 &pos4,
2893 &vel4,
2894 &wind_vector,
2895 (resolved_temp_c, resolved_press_hpa, base_ratio),
2896 );
2897
2898 position += (velocity + vel2 * 2.0 + vel3 * 2.0 + vel4) * (dt / 6.0);
2900 velocity += (acc1 + acc2 * 2.0 + acc3 * 2.0 + acc4) * (dt / 6.0);
2901 time += dt;
2902 self.validate_integration_state(&position, &velocity, time)?;
2903 }
2904
2905 let termination =
2906 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
2907
2908 self.annotate_drag_coefficients(&mut points, speed_of_sound);
2913
2914 let last_point = points.last().ok_or("No trajectory points generated")?;
2915
2916 let sampled_points = self.build_sampled_points(
2918 &points,
2919 max_height,
2920 transonic_distances,
2921 &mach_transitions,
2922 )?;
2923
2924 Ok(TrajectoryResult {
2925 max_range: last_point.position.x, max_height,
2927 time_of_flight: last_point.time,
2928 impact_velocity: last_point.velocity_magnitude,
2929 impact_energy: last_point.kinetic_energy,
2930 projectile_mass_kg: self.inputs.bullet_mass,
2931 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
2932 station_speed_of_sound_mps: speed_of_sound,
2933 termination,
2934 points,
2935 sampled_points,
2936 min_pitch_damping: if self.inputs.enable_pitch_damping {
2937 Some(min_pitch_damping)
2938 } else {
2939 None
2940 },
2941 transonic_mach,
2942 angular_state,
2943 max_yaw_angle: if self.inputs.enable_precession_nutation {
2944 Some(max_yaw_angle)
2945 } else {
2946 None
2947 },
2948 max_precession_angle: if self.inputs.enable_precession_nutation {
2949 Some(max_precession_angle)
2950 } else {
2951 None
2952 },
2953 aerodynamic_jump: aj_components,
2954 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2955 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2956 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2957 })
2958 }
2959
2960 fn solve_rk45(&self) -> Result<TrajectoryResult, BallisticsError> {
2961 let mut time = 0.0;
2963 let mut position = self.initial_position();
2967
2968 let aj_components = self.aerodynamic_jump_components();
2974 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2975 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2976 let mut velocity = Vector3::new(
2977 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2981
2982 let mut points = Vec::new();
2983 let mut max_height = position.y;
2984 let mut dt = 0.001; let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2989 self.resolved_atmosphere();
2990 let base_ratio = air_density / 1.225;
2995 let wind_vector =
3000 crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
3001
3002 let mut transonic_distances: Vec<f64> = Vec::new();
3004 let mut mach_transitions = MachTransitionTracker::default();
3005
3006 let mut min_pitch_damping = f64::INFINITY;
3011 let mut transonic_mach: Option<f64> = None;
3012 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
3013 self.inputs.bullet_model.as_deref().unwrap_or("default"),
3014 );
3015 let mut angular_state = if self.inputs.enable_precession_nutation {
3016 Some(AngularState {
3017 pitch_angle: 0.001,
3018 yaw_angle: 0.001,
3019 pitch_rate: 0.0,
3020 yaw_rate: 0.0,
3021 precession_angle: 0.0,
3022 nutation_phase: 0.0,
3023 })
3024 } else {
3025 None
3026 };
3027 let mut max_yaw_angle = 0.0;
3028 let mut max_precession_angle = 0.0;
3029
3030 while position.x < self.max_range
3031 && position.y > self.inputs.ground_threshold
3032 && time < TRAJECTORY_TIME_LIMIT_S
3033 {
3034 let velocity_magnitude = velocity.magnitude();
3036 let kinetic_energy = 0.5 * self.inputs.bullet_mass * velocity_magnitude.powi(2);
3037
3038 self.push_trajectory_point(
3039 &mut points,
3040 TrajectoryPoint {
3041 time,
3042 position,
3043 velocity_magnitude,
3044 kinetic_energy,
3045 drag_coefficient: None,
3046 },
3047 )?;
3048
3049 {
3052 let mach_here = if speed_of_sound > 0.0 {
3053 velocity_magnitude / speed_of_sound
3054 } else {
3055 0.0
3056 };
3057 mach_transitions.record_downward_crossings(
3058 mach_here,
3059 position.x,
3060 &mut transonic_distances,
3061 );
3062 }
3063
3064 if position.y > max_height {
3065 max_height = position.y;
3066 }
3067
3068 if self.inputs.enable_pitch_damping {
3071 let mach = velocity_magnitude / speed_of_sound;
3072 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
3073 transonic_mach = Some(mach);
3074 }
3075 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
3076 if pitch_damping < min_pitch_damping {
3077 min_pitch_damping = pitch_damping;
3078 }
3079 }
3080
3081 let accepted_step = self.adaptive_rk45_step(
3084 &position,
3085 &velocity,
3086 dt,
3087 &wind_vector,
3088 (resolved_temp_c, resolved_press_hpa, base_ratio),
3089 );
3090 debug_assert!(
3091 accepted_step.error <= RK45_TOLERANCE || accepted_step.used_dt <= RK45_MIN_DT
3092 );
3093
3094 if self.inputs.enable_precession_nutation {
3098 if let Some(ref mut state) = angular_state {
3099 let params = self.precession_nutation_params(
3100 velocity_magnitude,
3101 air_density,
3102 speed_of_sound,
3103 );
3104
3105 *state = calculate_combined_angular_motion(
3106 ¶ms,
3107 state,
3108 time,
3109 accepted_step.used_dt,
3110 0.001,
3111 );
3112
3113 if state.yaw_angle.abs() > max_yaw_angle {
3114 max_yaw_angle = state.yaw_angle.abs();
3115 }
3116 if state.precession_angle.abs() > max_precession_angle {
3117 max_precession_angle = state.precession_angle.abs();
3118 }
3119 }
3120 }
3121
3122 position = accepted_step.position;
3123 velocity = accepted_step.velocity;
3124 time += accepted_step.used_dt;
3125 self.validate_integration_state(&position, &velocity, time)?;
3126
3127 dt = accepted_step.next_dt;
3129 }
3130
3131 if points.is_empty() {
3133 return Err(BallisticsError::from("No trajectory points calculated"));
3134 }
3135
3136 let termination =
3138 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
3139
3140 self.annotate_drag_coefficients(&mut points, speed_of_sound);
3142
3143 let last_point = points.last().unwrap();
3144
3145 let sampled_points = self.build_sampled_points(
3147 &points,
3148 max_height,
3149 transonic_distances,
3150 &mach_transitions,
3151 )?;
3152
3153 Ok(TrajectoryResult {
3154 max_range: last_point.position.x, max_height,
3156 time_of_flight: last_point.time,
3157 impact_velocity: last_point.velocity_magnitude,
3158 impact_energy: last_point.kinetic_energy,
3159 projectile_mass_kg: self.inputs.bullet_mass,
3160 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
3161 station_speed_of_sound_mps: speed_of_sound,
3162 termination,
3163 points,
3164 sampled_points,
3165 min_pitch_damping: if self.inputs.enable_pitch_damping {
3166 Some(min_pitch_damping)
3167 } else {
3168 None
3169 },
3170 transonic_mach,
3171 angular_state,
3172 max_yaw_angle: if self.inputs.enable_precession_nutation {
3173 Some(max_yaw_angle)
3174 } else {
3175 None
3176 },
3177 max_precession_angle: if self.inputs.enable_precession_nutation {
3178 Some(max_precession_angle)
3179 } else {
3180 None
3181 },
3182 aerodynamic_jump: aj_components,
3183 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
3184 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
3185 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
3186 })
3187 }
3188
3189 fn adaptive_rk45_step(
3190 &self,
3191 position: &Vector3<f64>,
3192 velocity: &Vector3<f64>,
3193 initial_dt: f64,
3194 wind_vector: &Vector3<f64>,
3195 resolved_atmo: (f64, f64, f64),
3196 ) -> Rk45AcceptedStep {
3197 let mut trial_dt = initial_dt;
3198
3199 loop {
3200 let trial = self.rk45_step(
3201 position,
3202 velocity,
3203 trial_dt,
3204 wind_vector,
3205 RK45_TOLERANCE,
3206 resolved_atmo,
3207 );
3208 let next_dt = if trial.suggested_dt.is_finite() {
3213 (RK45_SAFETY_FACTOR * trial.suggested_dt).clamp(RK45_MIN_DT, RK45_MAX_DT)
3214 } else {
3215 RK45_MIN_DT
3216 };
3217
3218 if trial.error <= RK45_TOLERANCE || trial_dt <= RK45_MIN_DT {
3219 return Rk45AcceptedStep {
3220 position: trial.position,
3221 velocity: trial.velocity,
3222 used_dt: trial_dt,
3223 next_dt,
3224 error: trial.error,
3225 };
3226 }
3227
3228 trial_dt = next_dt;
3229 }
3230 }
3231
3232 fn rk45_step(
3233 &self,
3234 position: &Vector3<f64>,
3235 velocity: &Vector3<f64>,
3236 dt: f64,
3237 wind_vector: &Vector3<f64>,
3238 tolerance: f64,
3239 resolved_atmo: (f64, f64, f64), ) -> Rk45Trial {
3241 const A21: f64 = 1.0 / 5.0;
3243 const A31: f64 = 3.0 / 40.0;
3244 const A32: f64 = 9.0 / 40.0;
3245 const A41: f64 = 44.0 / 45.0;
3246 const A42: f64 = -56.0 / 15.0;
3247 const A43: f64 = 32.0 / 9.0;
3248 const A51: f64 = 19372.0 / 6561.0;
3249 const A52: f64 = -25360.0 / 2187.0;
3250 const A53: f64 = 64448.0 / 6561.0;
3251 const A54: f64 = -212.0 / 729.0;
3252 const A61: f64 = 9017.0 / 3168.0;
3253 const A62: f64 = -355.0 / 33.0;
3254 const A63: f64 = 46732.0 / 5247.0;
3255 const A64: f64 = 49.0 / 176.0;
3256 const A65: f64 = -5103.0 / 18656.0;
3257 const A71: f64 = 35.0 / 384.0;
3258 const A73: f64 = 500.0 / 1113.0;
3259 const A74: f64 = 125.0 / 192.0;
3260 const A75: f64 = -2187.0 / 6784.0;
3261 const A76: f64 = 11.0 / 84.0;
3262
3263 const B1: f64 = 35.0 / 384.0;
3265 const B3: f64 = 500.0 / 1113.0;
3266 const B4: f64 = 125.0 / 192.0;
3267 const B5: f64 = -2187.0 / 6784.0;
3268 const B6: f64 = 11.0 / 84.0;
3269
3270 const B1_ERR: f64 = 5179.0 / 57600.0;
3272 const B3_ERR: f64 = 7571.0 / 16695.0;
3273 const B4_ERR: f64 = 393.0 / 640.0;
3274 const B5_ERR: f64 = -92097.0 / 339200.0;
3275 const B6_ERR: f64 = 187.0 / 2100.0;
3276 const B7_ERR: f64 = 1.0 / 40.0;
3277
3278 let k1_v = self.calculate_acceleration(position, velocity, wind_vector, resolved_atmo);
3280 let k1_p = *velocity;
3281
3282 let p2 = position + dt * A21 * k1_p;
3283 let v2 = velocity + dt * A21 * k1_v;
3284 let k2_v = self.calculate_acceleration(&p2, &v2, wind_vector, resolved_atmo);
3285 let k2_p = v2;
3286
3287 let p3 = position + dt * (A31 * k1_p + A32 * k2_p);
3288 let v3 = velocity + dt * (A31 * k1_v + A32 * k2_v);
3289 let k3_v = self.calculate_acceleration(&p3, &v3, wind_vector, resolved_atmo);
3290 let k3_p = v3;
3291
3292 let p4 = position + dt * (A41 * k1_p + A42 * k2_p + A43 * k3_p);
3293 let v4 = velocity + dt * (A41 * k1_v + A42 * k2_v + A43 * k3_v);
3294 let k4_v = self.calculate_acceleration(&p4, &v4, wind_vector, resolved_atmo);
3295 let k4_p = v4;
3296
3297 let p5 = position + dt * (A51 * k1_p + A52 * k2_p + A53 * k3_p + A54 * k4_p);
3298 let v5 = velocity + dt * (A51 * k1_v + A52 * k2_v + A53 * k3_v + A54 * k4_v);
3299 let k5_v = self.calculate_acceleration(&p5, &v5, wind_vector, resolved_atmo);
3300 let k5_p = v5;
3301
3302 let p6 = position + dt * (A61 * k1_p + A62 * k2_p + A63 * k3_p + A64 * k4_p + A65 * k5_p);
3303 let v6 = velocity + dt * (A61 * k1_v + A62 * k2_v + A63 * k3_v + A64 * k4_v + A65 * k5_v);
3304 let k6_v = self.calculate_acceleration(&p6, &v6, wind_vector, resolved_atmo);
3305 let k6_p = v6;
3306
3307 let p7 = position + dt * (A71 * k1_p + A73 * k3_p + A74 * k4_p + A75 * k5_p + A76 * k6_p);
3308 let v7 = velocity + dt * (A71 * k1_v + A73 * k3_v + A74 * k4_v + A75 * k5_v + A76 * k6_v);
3309 let k7_v = self.calculate_acceleration(&p7, &v7, wind_vector, resolved_atmo);
3310 let k7_p = v7;
3311
3312 let new_pos = position + dt * (B1 * k1_p + B3 * k3_p + B4 * k4_p + B5 * k5_p + B6 * k6_p);
3314 let new_vel = velocity + dt * (B1 * k1_v + B3 * k3_v + B4 * k4_v + B5 * k5_v + B6 * k6_v);
3315
3316 let pos_err = position
3318 + dt * (B1_ERR * k1_p
3319 + B3_ERR * k3_p
3320 + B4_ERR * k4_p
3321 + B5_ERR * k5_p
3322 + B6_ERR * k6_p
3323 + B7_ERR * k7_p);
3324 let vel_err = velocity
3325 + dt * (B1_ERR * k1_v
3326 + B3_ERR * k3_v
3327 + B4_ERR * k4_v
3328 + B5_ERR * k5_v
3329 + B6_ERR * k6_v
3330 + B7_ERR * k7_v);
3331
3332 let error = cli_rk45_error_norm(position, velocity, &new_pos, &new_vel, &pos_err, &vel_err);
3334
3335 let dt_new = if error < tolerance {
3337 dt * (tolerance / error).powf(0.2).min(2.0)
3338 } else {
3339 dt * (tolerance / error).powf(0.25).max(0.1)
3340 };
3341
3342 Rk45Trial {
3343 position: new_pos,
3344 velocity: new_vel,
3345 suggested_dt: dt_new,
3346 error,
3347 }
3348 }
3349
3350 fn apply_cluster_bc_correction(&self, base_bc: f64, velocity_fps: f64) -> f64 {
3351 if let Some(ref cluster_bc) = self.cluster_bc {
3352 cluster_bc.apply_correction_for_drag_model(
3353 base_bc,
3354 self.inputs.caliber_inches,
3355 self.inputs.weight_grains,
3356 velocity_fps,
3357 self.inputs.bc_type,
3358 )
3359 } else {
3360 base_bc
3361 }
3362 }
3363
3364 fn calculate_acceleration(
3365 &self,
3366 position: &Vector3<f64>,
3367 velocity: &Vector3<f64>,
3368 wind_vector: &Vector3<f64>,
3369 resolved_atmo: (f64, f64, f64), ) -> Vector3<f64> {
3371 let actual_wind = if let Some(ref sock) = self.wind_sock {
3377 sock.vector_for_range_stateless(position.x)
3378 } else if self.inputs.enable_wind_shear {
3379 self.get_wind_at_altitude(position.y)
3380 } else {
3381 *wind_vector
3382 };
3383 let actual_wind =
3384 crate::derivatives::level_vector_to_shot_frame(actual_wind, self.inputs.shooting_angle);
3385
3386 let relative_velocity = velocity - actual_wind;
3387 let velocity_magnitude = relative_velocity.magnitude();
3388
3389 if velocity_magnitude < 0.001 {
3390 return self.gravity_acceleration();
3391 }
3392
3393 let (base_temp_c, base_press_hpa, station_ratio) = resolved_atmo;
3404
3405 let (drag_base_temp_c, drag_base_press_hpa, drag_base_ratio, drag_humidity_percent) =
3413 if let Some(ref sock) = self.atmo_sock {
3414 let (zone_temp_c, zone_press_hpa, zone_humidity) = sock.atmo_for_range(position.x);
3415 let zone_base_ratio = crate::atmosphere::calculate_air_density_cimp(
3416 zone_temp_c,
3417 zone_press_hpa,
3418 zone_humidity,
3419 ) / 1.225;
3420 (zone_temp_c, zone_press_hpa, zone_base_ratio, zone_humidity)
3421 } else {
3422 (
3423 base_temp_c,
3424 base_press_hpa,
3425 station_ratio,
3426 self.atmosphere.humidity,
3427 )
3428 };
3429 let local_alt = crate::atmosphere::shot_frame_altitude(
3430 self.atmosphere.altitude,
3431 position.x,
3432 position.y,
3433 self.inputs.shooting_angle,
3434 );
3435 let (air_density, speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
3436 local_alt,
3437 self.atmosphere.altitude,
3438 drag_base_temp_c,
3439 drag_base_press_hpa,
3440 drag_base_ratio,
3441 drag_humidity_percent,
3442 );
3443
3444 let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
3448
3449 let velocity_fps = velocity_magnitude * 3.28084;
3451
3452 let cd_to_retard = crate::constants::CD_TO_RETARD;
3457 let standard_factor = cd * cd_to_retard;
3458 let density_scale = air_density / 1.225; let a_drag_ft_s2 =
3462 (velocity_fps * velocity_fps) * standard_factor * density_scale / retard_denom;
3463 let a_drag_m_s2 = a_drag_ft_s2 * 0.3048; let drag_acceleration = -a_drag_m_s2 * (relative_velocity / velocity_magnitude);
3467
3468 let mut accel = drag_acceleration + self.gravity_acceleration();
3471
3472 if self.inputs.enable_coriolis {
3475 if let Some(lat_deg) = self.inputs.latitude {
3476 let omega_earth = 7.2921159e-5_f64; let lat = lat_deg.to_radians();
3478 let az = self.inputs.shot_azimuth; let omega = Vector3::new(
3485 omega_earth * lat.cos() * az.cos(), omega_earth * lat.sin(), -omega_earth * lat.cos() * az.sin(), );
3489 let omega = crate::derivatives::level_vector_to_shot_frame(
3490 omega,
3491 self.inputs.shooting_angle,
3492 );
3493 accel += -2.0 * omega.cross(velocity);
3498 }
3499 }
3500
3501 if self.inputs.enable_magnus
3508 && !self.inputs.use_enhanced_spin_drift
3509 && self.inputs.bullet_diameter > 0.0
3510 && self.inputs.twist_rate > 0.0
3511 {
3512 let diameter_m = self.inputs.bullet_diameter;
3513 let (spin_rad_s, spin_param) = crate::spin_drift::calculate_magnus_spin_state(
3514 self.inputs.muzzle_velocity,
3515 velocity_magnitude,
3516 self.inputs.twist_rate,
3517 diameter_m,
3518 );
3519 let mach = velocity_magnitude / speed_of_sound;
3521
3522 let d_in = self.inputs.bullet_diameter / 0.0254;
3524 let m_gr = self.inputs.bullet_mass / crate::constants::GRAINS_TO_KG;
3525 let l_in = if self.inputs.bullet_length > 0.0 {
3526 self.inputs.bullet_length / 0.0254
3527 } else {
3528 let est_m = crate::stability::estimate_bullet_length_m(
3530 self.inputs.bullet_diameter,
3531 self.inputs.bullet_mass,
3532 );
3533 if est_m > 0.0 {
3534 est_m / 0.0254
3535 } else {
3536 4.5 * d_in
3537 }
3538 };
3539 let sg = crate::spin_drift::calculate_dynamic_stability(
3543 m_gr,
3544 velocity_magnitude,
3545 spin_rad_s,
3546 d_in,
3547 l_in,
3548 air_density,
3549 );
3550
3551 let (yaw_rad, _) = crate::spin_drift::calculate_yaw_of_repose(
3553 sg,
3554 velocity_magnitude,
3555 spin_rad_s,
3556 0.0, 0.0, air_density,
3559 d_in,
3560 l_in,
3561 m_gr,
3562 mach,
3563 "match",
3564 false,
3565 );
3566
3567 let c_np = crate::derivatives::calculate_magnus_moment_coefficient(mach);
3569 let area = std::f64::consts::PI * (diameter_m / 2.0).powi(2);
3570 let magnus_force = 0.5
3571 * air_density
3572 * velocity_magnitude.powi(2)
3573 * area
3574 * c_np
3575 * spin_param
3576 * yaw_rad.sin();
3577
3578 if magnus_force.abs() > 1e-12 {
3582 if let Some(dir) = crate::derivatives::yaw_of_repose_magnus_direction(
3583 relative_velocity,
3584 self.gravity_acceleration(),
3585 self.inputs.is_twist_right,
3586 ) {
3587 accel += (magnus_force / self.inputs.bullet_mass) * dir;
3588 }
3589 }
3590 }
3591
3592 accel
3593 }
3594
3595 fn drag_terms(&self, velocity_magnitude: f64, speed_of_sound: f64) -> (f64, f64) {
3607 let cd = self.calculate_drag_coefficient(velocity_magnitude, speed_of_sound);
3608
3609 let velocity_fps = velocity_magnitude * 3.28084;
3610
3611 let (base_bc, bc_from_segments) = if let Some(segments) = self
3616 .inputs
3617 .bc_segments_data
3618 .as_ref()
3619 .filter(|segments| self.inputs.use_bc_segments && !segments.is_empty())
3620 {
3621 (
3623 crate::bc_estimation::velocity_segment_bc(
3624 velocity_fps,
3625 segments,
3626 self.inputs.bc_value,
3627 ),
3628 true,
3629 )
3630 } else if let Some(segments) = self
3631 .inputs
3632 .bc_segments
3633 .as_ref()
3634 .filter(|segments| !segments.is_empty())
3635 {
3636 (
3637 crate::derivatives::interpolated_bc(
3638 velocity_magnitude / speed_of_sound,
3639 segments,
3640 Some(&self.inputs),
3641 ),
3642 true,
3643 )
3644 } else {
3645 (self.inputs.bc_value, false)
3646 };
3647
3648 let effective_bc = if bc_from_segments {
3653 base_bc
3654 } else {
3655 self.apply_cluster_bc_correction(base_bc, velocity_fps)
3656 };
3657 let effective_bc = effective_bc.max(1e-6);
3660
3661 let retard_denom = if self.inputs.custom_drag_table.is_some() {
3666 self.inputs.custom_drag_denominator(effective_bc)
3667 } else {
3668 effective_bc
3669 };
3670
3671 (cd, retard_denom)
3672 }
3673
3674 pub fn effective_drag_coefficient(
3695 &self,
3696 velocity_magnitude: f64,
3697 speed_of_sound: f64,
3698 ) -> Option<f64> {
3699 if !velocity_magnitude.is_finite() || speed_of_sound <= 1e-9 {
3700 return None;
3701 }
3702 let sectional_density = self.inputs.sectional_density_lb_in2()?;
3703 let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
3704 if retard_denom <= 0.0 {
3705 return None;
3706 }
3707 let effective = cd * sectional_density / retard_denom;
3708 effective.is_finite().then_some(effective)
3709 }
3710
3711 fn annotate_drag_coefficients(&self, points: &mut [TrajectoryPoint], speed_of_sound: f64) {
3721 for point in points.iter_mut() {
3722 point.drag_coefficient =
3723 self.effective_drag_coefficient(point.velocity_magnitude, speed_of_sound);
3724 }
3725 }
3726
3727 fn calculate_drag_coefficient(&self, velocity: f64, speed_of_sound: f64) -> f64 {
3728 let mach = velocity / speed_of_sound;
3729
3730 if let Some(ref table) = self.inputs.custom_drag_table {
3734 return table.interpolate(mach) * self.inputs.cd_scale;
3739 }
3740
3741 crate::drag::get_drag_coefficient(mach, &self.inputs.bc_type)
3744 }
3745}
3746
3747#[derive(Debug, Clone)]
3749pub struct MonteCarloParams {
3750 pub num_simulations: usize,
3751 pub velocity_std_dev: f64,
3752 pub angle_std_dev: f64,
3753 pub bc_std_dev: f64,
3754 pub wind_speed_std_dev: f64,
3755 pub target_distance: Option<f64>,
3756 pub base_wind_speed: f64,
3757 pub base_wind_direction: f64,
3758 pub azimuth_std_dev: f64, }
3760
3761impl Default for MonteCarloParams {
3762 fn default() -> Self {
3763 Self {
3764 num_simulations: 1000,
3765 velocity_std_dev: 1.0,
3766 angle_std_dev: 0.001,
3767 bc_std_dev: 0.01,
3768 wind_speed_std_dev: 1.0,
3769 target_distance: None,
3770 base_wind_speed: 0.0,
3771 base_wind_direction: 0.0,
3772 azimuth_std_dev: 0.001, }
3774 }
3775}
3776
3777#[derive(Debug, Clone)]
3779pub struct MonteCarloResults {
3780 pub ranges: Vec<f64>,
3781 pub impact_velocities: Vec<f64>,
3782 pub impact_positions: Vec<Vector3<f64>>,
3788}
3789
3790pub const DEFAULT_HIT_RADIUS_M: f64 = 0.3;
3793
3794pub const TARGET_NOT_REACHED_SENTINEL_M: f64 = -1.0e9;
3800
3801impl MonteCarloResults {
3802 pub fn position_reached_target(position: &Vector3<f64>) -> bool {
3804 position.iter().all(|component| component.is_finite())
3805 && position.y != TARGET_NOT_REACHED_SENTINEL_M
3806 }
3807
3808 pub fn target_arrival_count(&self) -> usize {
3810 self.impact_positions
3811 .iter()
3812 .filter(|position| Self::position_reached_target(position))
3813 .count()
3814 }
3815
3816 pub fn target_shortfall_fraction(&self) -> f64 {
3819 if self.impact_positions.is_empty() {
3820 return 0.0;
3821 }
3822 (self.impact_positions.len() - self.target_arrival_count()) as f64
3823 / self.impact_positions.len() as f64
3824 }
3825
3826 pub fn target_plane_cep(&self) -> Option<f64> {
3832 let mut radial_misses: Vec<f64> = self
3833 .impact_positions
3834 .iter()
3835 .filter(|position| Self::position_reached_target(position))
3836 .map(Vector3::norm)
3837 .filter(|miss| miss.is_finite())
3838 .collect();
3839 radial_misses.sort_by(f64::total_cmp);
3840 if radial_misses.is_empty() {
3841 None
3842 } else {
3843 Some(radial_misses[radial_misses.len() / 2])
3844 }
3845 }
3846
3847 pub fn hit_probability(&self, hit_radius_m: f64) -> f64 {
3856 if self.impact_positions.is_empty() {
3857 return 0.0;
3858 }
3859 let hits = self
3860 .impact_positions
3861 .iter()
3862 .filter(|position| Self::position_is_hit(position, hit_radius_m))
3863 .count();
3864 hits as f64 / self.impact_positions.len() as f64
3865 }
3866
3867 pub fn position_is_hit(position: &Vector3<f64>, hit_radius_m: f64) -> bool {
3879 Self::position_reached_target(position) && position.norm() < hit_radius_m
3880 }
3881
3882 pub fn hit_probability_wilson(
3904 &self,
3905 hit_radius_m: f64,
3906 level: ConfidenceLevel,
3907 ) -> (f64, (f64, f64), u64) {
3908 let trials = self.impact_positions.len() as u64;
3909 let hits = self
3910 .impact_positions
3911 .iter()
3912 .filter(|position| Self::position_is_hit(position, hit_radius_m))
3913 .count() as u64;
3914 (
3915 self.hit_probability(hit_radius_m),
3916 wilson_interval(hits, trials, level),
3917 trials,
3918 )
3919 }
3920
3921 pub fn rect_hit_probability(&self, width_m: f64, height_m: f64) -> f64 {
3933 let dimensions_invalid = width_m.is_nan()
3934 || width_m <= 0.0
3935 || height_m.is_nan()
3936 || height_m <= 0.0;
3937 if self.impact_positions.is_empty() || dimensions_invalid {
3938 return 0.0;
3939 }
3940 let half_width = width_m / 2.0;
3941 let half_height = height_m / 2.0;
3942 let hits = self
3943 .impact_positions
3944 .iter()
3945 .filter(|position| {
3946 Self::position_reached_target(position)
3947 && position.z.abs() <= half_width
3948 && position.y.abs() <= half_height
3949 })
3950 .count();
3951 hits as f64 / self.impact_positions.len() as f64
3952 }
3953}
3954
3955fn wind_from_signed_speed_sample(
3956 signed_speed: f64,
3957 sampled_direction: f64,
3958 vertical_speed: f64,
3959) -> WindConditions {
3960 if signed_speed < 0.0 {
3965 WindConditions {
3966 speed: -signed_speed,
3967 direction: sampled_direction + std::f64::consts::PI,
3968 vertical_speed,
3969 }
3970 } else {
3971 WindConditions {
3972 speed: signed_speed,
3973 direction: sampled_direction,
3974 vertical_speed,
3975 }
3976 }
3977}
3978
3979struct MonteCarloWindSampler {
3980 speed: rand_distr::Normal<f64>,
3981 direction: rand_distr::Normal<f64>,
3982 vertical_speed: f64,
3984}
3985
3986impl MonteCarloWindSampler {
3987 fn new(
3988 base_wind: &WindConditions,
3989 wind_speed_std_dev: f64,
3990 wind_direction_std_dev: f64,
3991 ) -> Result<Self, BallisticsError> {
3992 use rand_distr::Normal;
3993
3994 if !wind_direction_std_dev.is_finite() || wind_direction_std_dev < 0.0 {
3995 return Err("Wind direction standard deviation must be finite and non-negative".into());
3996 }
3997
3998 let speed = Normal::new(base_wind.speed, wind_speed_std_dev)
3999 .map_err(|e| format!("Invalid wind speed distribution: {e}"))?;
4000 let direction = Normal::new(base_wind.direction, wind_direction_std_dev)
4001 .map_err(|e| format!("Invalid wind direction distribution: {e}"))?;
4002 Ok(Self { speed, direction, vertical_speed: base_wind.vertical_speed })
4003 }
4004
4005 fn sample<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> WindConditions {
4006 use rand_distr::Distribution;
4007
4008 wind_from_signed_speed_sample(
4009 self.speed.sample(rng),
4010 self.direction.sample(rng),
4011 self.vertical_speed,
4012 )
4013 }
4014}
4015
4016#[derive(Debug, Clone, Copy)]
4019struct TrialOutcome {
4020 range: f64,
4023 impact_velocity: f64,
4025 impact_position: Vector3<f64>,
4028}
4029
4030struct MonteCarloTrialSampler {
4050 base_inputs: BallisticInputs,
4051 atmosphere: AtmosphericConditions,
4052 solver_max_range: f64,
4053 target_distance: f64,
4056 baseline_at_target: Vector3<f64>,
4059 velocity_delta_dist: rand_distr::Normal<f64>,
4060 angle_dist: rand_distr::Normal<f64>,
4061 bc_dist: rand_distr::Normal<f64>,
4062 wind_sampler: MonteCarloWindSampler,
4063 azimuth_dist: rand_distr::Normal<f64>,
4064}
4065
4066impl MonteCarloTrialSampler {
4067 fn new(
4071 base_inputs: BallisticInputs,
4072 base_wind: &WindConditions,
4073 params: &MonteCarloParams,
4074 wind_direction_std_dev: f64,
4075 ) -> Result<Self, BallisticsError> {
4076 use rand_distr::Normal;
4077
4078 let atmosphere = AtmosphericConditions {
4079 temperature: base_inputs.temperature,
4080 pressure: base_inputs.pressure,
4081 humidity: base_inputs.humidity_percent(),
4082 altitude: base_inputs.altitude,
4083 };
4084 let target_hint = params
4085 .target_distance
4086 .unwrap_or(base_inputs.target_distance);
4087 let solver_max_range = target_hint.max(1000.0) * 2.0;
4088
4089 let mut baseline_solver =
4091 TrajectorySolver::new(base_inputs.clone(), base_wind.clone(), atmosphere.clone());
4092 baseline_solver.set_max_range(solver_max_range);
4093 let baseline_result = baseline_solver.solve()?;
4094
4095 let target_distance = params.target_distance.unwrap_or(baseline_result.max_range);
4097
4098 let baseline_at_target = baseline_result
4100 .position_at_range(target_distance)
4101 .ok_or("Could not interpolate baseline at target distance")?;
4102
4103 let velocity_delta_dist = Normal::new(0.0, params.velocity_std_dev)
4108 .map_err(|e| format!("Invalid velocity distribution: {}", e))?;
4109 let angle_dist = Normal::new(base_inputs.muzzle_angle, params.angle_std_dev)
4110 .map_err(|e| format!("Invalid angle distribution: {}", e))?;
4111 let bc_dist = Normal::new(base_inputs.bc_value, params.bc_std_dev)
4112 .map_err(|e| format!("Invalid BC distribution: {}", e))?;
4113 let wind_sampler = MonteCarloWindSampler::new(
4116 base_wind,
4117 params.wind_speed_std_dev,
4118 wind_direction_std_dev,
4119 )?;
4120 let azimuth_dist = Normal::new(base_inputs.azimuth_angle, params.azimuth_std_dev)
4121 .map_err(|e| format!("Invalid azimuth distribution: {}", e))?;
4122
4123 Ok(Self {
4124 base_inputs,
4125 atmosphere,
4126 solver_max_range,
4127 target_distance,
4128 baseline_at_target,
4129 velocity_delta_dist,
4130 angle_dist,
4131 bc_dist,
4132 wind_sampler,
4133 azimuth_dist,
4134 })
4135 }
4136
4137 fn sample_one_trial<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> Option<TrialOutcome> {
4147 use rand_distr::Distribution;
4148
4149 let mut inputs = self.base_inputs.clone();
4151 let muzzle_velocity_delta = self.velocity_delta_dist.sample(&mut *rng);
4152 inputs.muzzle_angle = self.angle_dist.sample(&mut *rng);
4153 inputs.bc_value = self.bc_dist.sample(&mut *rng).max(0.01);
4154 inputs.azimuth_angle = self.azimuth_dist.sample(&mut *rng); let wind = self.wind_sampler.sample(&mut *rng);
4158
4159 let mut solver = TrajectorySolver::new(inputs, wind, self.atmosphere.clone());
4163 solver.inputs.muzzle_velocity =
4164 (solver.inputs.muzzle_velocity + muzzle_velocity_delta).max(0.0);
4165 solver.set_max_range(self.solver_max_range);
4166 let result = solver.solve().ok()?;
4168
4169 let impact_position = if result.max_range < self.target_distance {
4175 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)
4178 } else {
4179 let pos_at_target = result.position_at_range(self.target_distance)?;
4181 Vector3::new(
4186 0.0,
4187 pos_at_target.y - self.baseline_at_target.y,
4188 pos_at_target.z - self.baseline_at_target.z,
4189 )
4190 };
4191
4192 Some(TrialOutcome {
4193 range: result.max_range,
4194 impact_velocity: result.impact_velocity,
4195 impact_position,
4196 })
4197 }
4198}
4199
4200pub fn run_monte_carlo(
4202 base_inputs: BallisticInputs,
4203 params: MonteCarloParams,
4204) -> Result<MonteCarloResults, BallisticsError> {
4205 run_monte_carlo_with_direction_std_dev(base_inputs, params, 0.0)
4206}
4207
4208pub fn run_monte_carlo_with_direction_std_dev(
4213 base_inputs: BallisticInputs,
4214 params: MonteCarloParams,
4215 wind_direction_std_dev: f64,
4216) -> Result<MonteCarloResults, BallisticsError> {
4217 let base_wind = WindConditions {
4218 speed: params.base_wind_speed,
4219 direction: params.base_wind_direction,
4220 vertical_speed: 0.0,
4221 };
4222 run_monte_carlo_with_wind_and_direction_std_dev(
4223 base_inputs,
4224 base_wind,
4225 params,
4226 wind_direction_std_dev,
4227 )
4228}
4229
4230pub fn run_monte_carlo_with_wind(
4232 base_inputs: BallisticInputs,
4233 base_wind: WindConditions,
4234 params: MonteCarloParams,
4235) -> Result<MonteCarloResults, BallisticsError> {
4236 run_monte_carlo_with_wind_and_direction_std_dev(base_inputs, base_wind, params, 0.0)
4237}
4238
4239pub fn run_monte_carlo_with_wind_and_direction_std_dev(
4244 base_inputs: BallisticInputs,
4245 base_wind: WindConditions,
4246 params: MonteCarloParams,
4247 wind_direction_std_dev: f64,
4248) -> Result<MonteCarloResults, BallisticsError> {
4249 let mut rng = rand::rng();
4250 run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
4251 base_inputs,
4252 base_wind,
4253 params,
4254 wind_direction_std_dev,
4255 &mut rng,
4256 )
4257}
4258
4259pub fn run_monte_carlo_with_wind_and_direction_std_dev_seeded(
4266 base_inputs: BallisticInputs,
4267 base_wind: WindConditions,
4268 params: MonteCarloParams,
4269 wind_direction_std_dev: f64,
4270 seed: u64,
4271) -> Result<MonteCarloResults, BallisticsError> {
4272 use rand::{rngs::StdRng, SeedableRng};
4273 let mut rng = StdRng::seed_from_u64(seed);
4274 run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
4275 base_inputs,
4276 base_wind,
4277 params,
4278 wind_direction_std_dev,
4279 &mut rng,
4280 )
4281}
4282
4283fn run_monte_carlo_with_wind_and_direction_std_dev_using_rng<R: rand::Rng + ?Sized>(
4284 base_inputs: BallisticInputs,
4285 base_wind: WindConditions,
4286 params: MonteCarloParams,
4287 wind_direction_std_dev: f64,
4288 rng: &mut R,
4289) -> Result<MonteCarloResults, BallisticsError> {
4290 let mut ranges = Vec::new();
4291 let mut impact_velocities = Vec::new();
4292 let mut impact_positions = Vec::new();
4293
4294 let sampler = MonteCarloTrialSampler::new(
4295 base_inputs,
4296 &base_wind,
4297 ¶ms,
4298 wind_direction_std_dev,
4299 )?;
4300
4301 for _ in 0..params.num_simulations {
4302 if let Some(outcome) = sampler.sample_one_trial(rng) {
4306 ranges.push(outcome.range);
4307 impact_velocities.push(outcome.impact_velocity);
4308 impact_positions.push(outcome.impact_position);
4309 }
4310 }
4311
4312 if ranges.is_empty() {
4313 return Err("No successful simulations".into());
4314 }
4315
4316 Ok(MonteCarloResults {
4317 ranges,
4318 impact_velocities,
4319 impact_positions,
4320 })
4321}
4322
4323pub const MC_ADAPTIVE_SCHEMA_VERSION_V1: u32 = 1;
4325
4326pub const MC_ADAPTIVE_METHOD_V1: &str = "anytime_beta_binomial_mixture_cs_v1";
4333
4334pub const MC_ADAPTIVE_ASSUMPTIONS_V1: [&str; 4] = [
4348 "Sampling uncertainty only: intervals cover Monte Carlo sampling error, not model error in the trajectory solver or its inputs.",
4349 "Anytime-valid stopping: the beta-binomial mixture confidence sequence keeps its coverage guarantee despite stopping the moment the target half-width is met.",
4350 "Input dispersions are the independent normal distributions declared in MonteCarloParams; correlations between inputs are not modeled.",
4351 "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.",
4352];
4353
4354#[derive(Debug, Clone)]
4361pub struct McConvergence {
4362 pub level: ConfidenceLevel,
4364 pub target_half_width: f64,
4367 pub min_samples: u64,
4371 pub max_samples: u64,
4374 pub batch_size: u64,
4377}
4378
4379impl Default for McConvergence {
4380 fn default() -> Self {
4381 Self {
4382 level: ConfidenceLevel::P95,
4383 target_half_width: 0.02,
4384 min_samples: 1_000,
4385 max_samples: 100_000,
4386 batch_size: 500,
4387 }
4388 }
4389}
4390
4391impl McConvergence {
4392 pub fn validate(&self) -> Result<(), String> {
4401 if !self.target_half_width.is_finite() || self.target_half_width <= 0.0 {
4402 return Err(format!(
4403 "McConvergence.target_half_width must be a finite value greater than zero (got {})",
4404 self.target_half_width
4405 ));
4406 }
4407 if self.batch_size == 0 {
4408 return Err("McConvergence.batch_size must be greater than zero".to_string());
4409 }
4410 if self.max_samples == 0 {
4411 return Err("McConvergence.max_samples must be greater than zero".to_string());
4412 }
4413 if self.max_samples < self.min_samples {
4414 return Err(format!(
4415 "McConvergence.max_samples ({}) must be at least McConvergence.min_samples ({})",
4416 self.max_samples, self.min_samples
4417 ));
4418 }
4419 Ok(())
4420 }
4421}
4422
4423#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize)]
4425#[serde(rename_all = "snake_case")]
4426pub enum McStopReason {
4427 TargetHalfWidthMet,
4430 MaxSamplesReached,
4433}
4434
4435#[derive(Debug, Clone, serde::Serialize)]
4442pub struct AdaptiveMcReportV1 {
4443 pub schema_version: u32,
4445 pub method: String,
4447 pub assumptions: Vec<String>,
4449 pub confidence_percent: u32,
4451 pub hit_probability: f64,
4454 pub ci_low: f64,
4456 pub ci_high: f64,
4458 pub samples: u64,
4460 pub attempts: u64,
4471 pub arrivals: u64,
4486 pub stop_reason: McStopReason,
4488 pub hit_radius_m: f64,
4490 pub target_distance_m: f64,
4492 pub mean_impact_velocity_mps: f64,
4497 pub std_impact_velocity_mps: f64,
4499 pub mean_drop_at_target_m: f64,
4507 pub std_drop_at_target_m: f64,
4510 pub mean_wind_drift_at_target_m: f64,
4513 pub std_wind_drift_at_target_m: f64,
4516}
4517
4518pub fn run_monte_carlo_adaptive_seeded(
4593 base_inputs: &BallisticInputs,
4594 base_wind: &WindConditions,
4595 params: &MonteCarloParams,
4596 convergence: &McConvergence,
4597 hit_radius_m: f64,
4598 seed: u64,
4599) -> Result<AdaptiveMcReportV1, String> {
4600 use rand::{rngs::StdRng, SeedableRng};
4601
4602 convergence.validate()?;
4603
4604 let sampler = MonteCarloTrialSampler::new(base_inputs.clone(), base_wind, params, 0.0)
4607 .map_err(|e| e.to_string())?;
4608
4609 let mut rng = StdRng::seed_from_u64(seed);
4610 let mut hits_cs = BernoulliConfidenceSequence::new(convergence.level);
4611 let mut impact_velocity = Welford::new();
4612 let mut drop_at_target = Welford::new();
4613 let mut drift_at_target = Welford::new();
4614
4615 let mut attempts: u64 = 0;
4616 let mut stop_reason = McStopReason::MaxSamplesReached;
4617
4618 while attempts < convergence.max_samples {
4619 let batch = convergence.batch_size.min(convergence.max_samples - attempts);
4622 let mut batch_hits: u64 = 0;
4623 let mut batch_trials: u64 = 0;
4624
4625 for _ in 0..batch {
4626 attempts += 1;
4627 let Some(outcome) = sampler.sample_one_trial(&mut rng) else {
4628 continue; };
4630 batch_trials += 1;
4631 if MonteCarloResults::position_is_hit(&outcome.impact_position, hit_radius_m) {
4632 batch_hits += 1;
4633 }
4634 if MonteCarloResults::position_reached_target(&outcome.impact_position) {
4638 impact_velocity.push(outcome.impact_velocity);
4639 drop_at_target.push(outcome.impact_position.y);
4640 drift_at_target.push(outcome.impact_position.z);
4641 }
4642 }
4643
4644 hits_cs.update_batch(batch_hits, batch_trials);
4645
4646 if hits_cs.trials() >= convergence.min_samples
4647 && hits_cs.half_width() <= convergence.target_half_width
4648 {
4649 stop_reason = McStopReason::TargetHalfWidthMet;
4650 break;
4651 }
4652 }
4653
4654 let samples = hits_cs.trials();
4655 if samples == 0 {
4656 return Err("No successful simulations".to_string());
4657 }
4658 let (ci_low, ci_high) = hits_cs.bounds();
4659
4660 Ok(AdaptiveMcReportV1 {
4661 schema_version: MC_ADAPTIVE_SCHEMA_VERSION_V1,
4662 method: MC_ADAPTIVE_METHOD_V1.to_string(),
4663 assumptions: MC_ADAPTIVE_ASSUMPTIONS_V1
4664 .iter()
4665 .map(|s| s.to_string())
4666 .collect(),
4667 confidence_percent: convergence.level.as_percent(),
4668 hit_probability: hits_cs.successes() as f64 / samples as f64,
4669 ci_low,
4670 ci_high,
4671 samples,
4672 attempts,
4673 arrivals: drop_at_target.count(),
4677 stop_reason,
4678 hit_radius_m,
4679 target_distance_m: sampler.target_distance,
4680 mean_impact_velocity_mps: impact_velocity.mean(),
4681 std_impact_velocity_mps: impact_velocity.sample_std(),
4682 mean_drop_at_target_m: drop_at_target.mean(),
4683 std_drop_at_target_m: drop_at_target.sample_std(),
4684 mean_wind_drift_at_target_m: drift_at_target.mean(),
4685 std_wind_drift_at_target_m: drift_at_target.sample_std(),
4686 })
4687}
4688
4689pub fn calculate_zero_angle(
4691 inputs: BallisticInputs,
4692 target_distance: f64,
4693 target_height: f64,
4694) -> Result<f64, BallisticsError> {
4695 calculate_zero_angle_with_conditions(
4696 inputs,
4697 target_distance,
4698 target_height,
4699 WindConditions::default(),
4700 AtmosphericConditions::default(),
4701 )
4702}
4703
4704pub fn calculate_zero_angle_with_conditions(
4705 inputs: BallisticInputs,
4706 target_distance: f64,
4707 target_height: f64,
4708 wind: WindConditions,
4709 atmosphere: AtmosphericConditions,
4710) -> Result<f64, BallisticsError> {
4711 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
4712 solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
4713}
4714
4715pub fn calculate_zero_angle_with_resolved_conditions(
4721 inputs: BallisticInputs,
4722 target_distance: f64,
4723 target_height: f64,
4724 wind: WindConditions,
4725 atmosphere: AtmosphericConditions,
4726) -> Result<f64, BallisticsError> {
4727 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
4728 solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
4729}
4730
4731pub const ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M: f64 = 2000.0;
4739
4740pub fn calculate_zero_range_from_angle_with_conditions(
4755 inputs: BallisticInputs,
4756 zero_angle_rad: f64,
4757 target_height: f64,
4758 wind: WindConditions,
4759 atmosphere: AtmosphericConditions,
4760) -> Result<ZeroCrossings, BallisticsError> {
4761 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
4762 solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
4763 solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
4764}
4765
4766pub fn calculate_zero_range_from_angle_with_resolved_conditions(
4772 inputs: BallisticInputs,
4773 zero_angle_rad: f64,
4774 target_height: f64,
4775 wind: WindConditions,
4776 atmosphere: AtmosphericConditions,
4777) -> Result<ZeroCrossings, BallisticsError> {
4778 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
4779 solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
4780 solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
4781}
4782
4783#[derive(Debug, Clone, Copy, PartialEq, Eq)]
4785pub enum BcFitMode {
4786 Drop,
4788 Velocity,
4791}
4792
4793#[derive(Debug, Clone, Copy)]
4795pub struct BcEstimate {
4796 pub bc: f64,
4798 pub rms_error: f64,
4800 pub drag_model: DragModel,
4802 pub mode: BcFitMode,
4804 pub at_bound: bool,
4808}
4809
4810fn fit_value_at(
4818 points: &[TrajectoryPoint],
4819 target_dist: f64,
4820 mode: BcFitMode,
4821 drop_offset: f64,
4822) -> Option<f64> {
4823 let val = |p: &TrajectoryPoint| match mode {
4824 BcFitMode::Drop => drop_offset - p.position.y,
4825 BcFitMode::Velocity => p.velocity_magnitude,
4826 };
4827 for i in 0..points.len() {
4828 if points[i].position.x >= target_dist {
4829 if i == 0 {
4830 return Some(val(&points[0]));
4831 }
4832 let p1 = &points[i - 1];
4833 let p2 = &points[i];
4834 let dx = p2.position.x - p1.position.x;
4835 if dx.abs() < 1e-9 {
4836 return Some(val(p2));
4837 }
4838 let t = (target_dist - p1.position.x) / dx;
4839 return Some(val(p1) + t * (val(p2) - val(p1)));
4840 }
4841 }
4842 None
4843}
4844
4845fn fit_residual_sse(
4846 trajectory: &[TrajectoryPoint],
4847 observations: &[(f64, f64)],
4848 mode: BcFitMode,
4849 drop_offset: f64,
4850) -> Option<f64> {
4851 if observations.is_empty() {
4852 return None;
4853 }
4854 let mut total = 0.0;
4855 for (target_dist, target_val) in observations {
4856 let value = fit_value_at(trajectory, *target_dist, mode, drop_offset)?;
4859 let error = value - target_val;
4860 total += error * error;
4861 }
4862 Some(total)
4863}
4864
4865#[allow(clippy::too_many_arguments)] pub fn estimate_bc_fit(
4883 velocity: f64,
4884 mass: f64,
4885 diameter: f64,
4886 points: &[(f64, f64)],
4887 drag_model: DragModel,
4888 mode: BcFitMode,
4889 atmosphere: AtmosphericConditions,
4890 zero_range: Option<f64>,
4891 sight_height: f64,
4892) -> Result<BcEstimate, BallisticsError> {
4893 if points.is_empty() {
4894 return Err(BallisticsError::from(
4895 "No data points provided for BC estimation.".to_string(),
4896 ));
4897 }
4898 let max_dist = points.iter().map(|(d, _)| *d).fold(0.0_f64, f64::max);
4899 let drop_offset = if zero_range.is_some() { sight_height } else { 0.0 };
4902
4903 let sse = |bc_value: f64| -> Option<f64> {
4905 let mut inputs = BallisticInputs {
4906 muzzle_velocity: velocity,
4907 bc_value,
4908 bc_type: drag_model,
4909 bullet_mass: mass,
4910 bullet_diameter: diameter,
4911 sight_height,
4912 ..Default::default()
4913 };
4914 if let Some(zr) = zero_range {
4917 let za = calculate_zero_angle_with_conditions(
4923 inputs.clone(),
4924 zr,
4925 sight_height,
4926 WindConditions::default(),
4927 atmosphere.clone(),
4928 )
4929 .ok()?;
4930 inputs.muzzle_angle = za;
4931 }
4932 let mut solver =
4933 TrajectorySolver::new(inputs, WindConditions::default(), atmosphere.clone());
4934 solver.set_max_range(max_dist * 1.5);
4935 let result = solver.solve().ok()?;
4936 fit_residual_sse(&result.points, points, mode, drop_offset)
4937 };
4938
4939 let (bc_min, bc_max) = match drag_model {
4943 DragModel::G7 => (0.05, 0.70),
4944 _ => (0.10, 1.20),
4945 };
4946
4947 let mut best_bc = f64::NAN;
4949 let mut best_sse = f64::MAX;
4950 let mut bc = bc_min;
4951 while bc <= bc_max + 1e-9 {
4952 if let Some(s) = sse(bc) {
4953 if s < best_sse {
4954 best_sse = s;
4955 best_bc = bc;
4956 }
4957 }
4958 bc += 0.01;
4959 }
4960 if !best_bc.is_finite() {
4961 return Err(BallisticsError::from(
4962 "Unable to estimate BC from provided data. Check that the values and units are correct."
4963 .to_string(),
4964 ));
4965 }
4966
4967 let lo = (best_bc - 0.01).max(bc_min);
4969 let hi = (best_bc + 0.01).min(bc_max);
4970 let mut bc = lo;
4971 while bc <= hi + 1e-9 {
4972 if let Some(s) = sse(bc) {
4973 if s < best_sse {
4974 best_sse = s;
4975 best_bc = bc;
4976 }
4977 }
4978 bc += 0.001;
4979 }
4980
4981 let at_bound = best_bc <= bc_min + 0.011 || best_bc >= bc_max - 0.011;
4984 let rms_error = (best_sse / points.len() as f64).sqrt();
4987 Ok(BcEstimate {
4988 bc: best_bc,
4989 rms_error,
4990 drag_model,
4991 mode,
4992 at_bound,
4993 })
4994}
4995
4996pub fn estimate_bc_from_trajectory(
4999 velocity: f64,
5000 mass: f64,
5001 diameter: f64,
5002 points: &[(f64, f64)], ) -> Result<f64, BallisticsError> {
5004 estimate_bc_fit(
5005 velocity,
5006 mass,
5007 diameter,
5008 points,
5009 DragModel::G1,
5010 BcFitMode::Drop,
5011 AtmosphericConditions::default(),
5012 None,
5013 0.05,
5014 )
5015 .map(|e| e.bc)
5016}
5017
5018use rand;
5020use rand_distr;
5021
5022#[cfg(test)]
5023mod mba737_powder_resolution_tests {
5024 use super::*;
5025
5026 #[test]
5027 fn linear_model_cold_powder_subtracts() {
5028 let v = resolve_powder_adjusted_velocity(823.0, 11.1, true, 0.5486, 21.1, None, None);
5030 assert!((v - (823.0 + 0.5486 * (11.1 - 21.1))).abs() < 1e-12);
5031 assert!(v < 823.0);
5032 }
5033
5034 #[test]
5035 fn linear_model_hot_powder_adds() {
5036 let v = resolve_powder_adjusted_velocity(823.0, 31.1, true, 0.5486, 21.1, None, None);
5037 assert!((v - (823.0 + 0.5486 * 10.0)).abs() < 1e-12);
5038 }
5039
5040 #[test]
5041 fn disabled_flag_is_passthrough() {
5042 let v = resolve_powder_adjusted_velocity(823.0, 40.0, false, 0.5486, 21.1, None, None);
5043 assert_eq!(v, 823.0);
5044 }
5045
5046 #[test]
5047 fn curve_overrides_linear_and_interpolates_at_powder_temp() {
5048 let curve = [(4.4, 798.6), (21.1, 823.0), (37.8, 841.2)];
5049 let v = resolve_powder_adjusted_velocity(823.0, 30.0, true, 99.0, 21.1, Some(&curve), Some(4.4));
5051 assert!((v - 798.6).abs() < 1e-9);
5052 }
5053
5054 #[test]
5055 fn curve_falls_back_to_ambient_and_clamps() {
5056 let curve = [(4.4, 798.6), (37.8, 841.2)];
5057 let v = resolve_powder_adjusted_velocity(823.0, -40.0, true, 1.0, 21.1, Some(&curve), None);
5059 assert!((v - 798.6).abs() < 1e-9);
5060 let v_hot = resolve_powder_adjusted_velocity(823.0, 60.0, true, 1.0, 21.1, Some(&curve), None);
5061 assert!((v_hot - 841.2).abs() < 1e-9);
5062 }
5063
5064 #[test]
5065 fn empty_curve_suppresses_linear_fallback() {
5066 let v = resolve_powder_adjusted_velocity(823.0, 40.0, true, 0.5486, 21.1, Some(&[]), None);
5068 assert_eq!(v, 823.0);
5069 }
5070
5071 #[test]
5072 fn sweep_huge_range_errors_instead_of_overflowing() {
5073 assert!(parse_powder_sweep("0:1e20:1").is_err());
5076 assert!(parse_powder_sweep("0:1e308:1e-3").is_err());
5077 }
5078
5079 #[test]
5080 fn sweep_fractional_step_keeps_end_row() {
5081 let rows = parse_powder_sweep("0:0.3:0.1").unwrap();
5083 assert_eq!(rows.len(), 4);
5084 assert!((rows[3] - 0.3).abs() < 1e-9);
5085 }
5086
5087 #[test]
5088 fn solver_and_helper_agree_on_linear_model() {
5089 let inputs = BallisticInputs {
5091 use_powder_sensitivity: true,
5092 powder_temp_sensitivity: 0.5486,
5093 powder_temp: 21.1,
5094 temperature: 4.4,
5095 ..Default::default()
5096 };
5097 let expected = resolve_powder_adjusted_velocity(
5098 inputs.muzzle_velocity,
5099 inputs.temperature,
5100 true,
5101 0.5486,
5102 21.1,
5103 None,
5104 None,
5105 );
5106 let solver = TrajectorySolver::new(
5107 inputs,
5108 WindConditions::default(),
5109 AtmosphericConditions::default(),
5110 );
5111 assert!((solver.inputs.muzzle_velocity - expected).abs() < 1e-12);
5112 }
5113}
5114
5115#[cfg(test)]
5116mod mba1302_solver_seam_tests {
5117 use super::*;
5118 use crate::wind::WindSegment;
5119
5120 #[test]
5121 fn authoritative_station_atmosphere_preserves_explicit_standard_values_at_altitude() {
5122 let atmosphere = AtmosphericConditions {
5123 temperature: 15.0,
5124 pressure: 1013.25,
5125 humidity: 50.0,
5126 altitude: 2_000.0,
5127 };
5128 let legacy = TrajectorySolver::new(
5129 BallisticInputs::default(),
5130 WindConditions::default(),
5131 atmosphere.clone(),
5132 );
5133 let authoritative = TrajectorySolver::new_with_resolved_station_atmosphere(
5134 BallisticInputs::default(),
5135 WindConditions::default(),
5136 atmosphere,
5137 );
5138
5139 let (legacy_density, _, legacy_temp_c, legacy_pressure_hpa) = legacy.resolved_atmosphere();
5140 let (authoritative_density, _, authoritative_temp_c, authoritative_pressure_hpa) =
5141 authoritative.resolved_atmosphere();
5142 let (icao_temp_k, icao_pressure_pa) =
5143 crate::atmosphere::calculate_icao_standard_atmosphere(2_000.0);
5144 let (expected_authoritative_density, _) =
5145 crate::atmosphere::calculate_atmosphere(2_000.0, Some(15.0), Some(1013.25), 50.0);
5146
5147 assert!((legacy_temp_c - (icao_temp_k - 273.15)).abs() < 1e-12);
5148 assert!((legacy_pressure_hpa - icao_pressure_pa / 100.0).abs() < 1e-12);
5149 assert_eq!(authoritative_temp_c.to_bits(), 15.0_f64.to_bits());
5150 assert_eq!(authoritative_pressure_hpa.to_bits(), 1013.25_f64.to_bits());
5151 assert_eq!(
5152 authoritative_density.to_bits(),
5153 expected_authoritative_density.to_bits()
5154 );
5155 assert!(
5156 (authoritative_density - legacy_density).abs() > 0.1,
5157 "explicit standard values at altitude must differ from ICAO-at-altitude: explicit={authoritative_density}, ICAO={legacy_density}"
5158 );
5159 }
5160
5161 #[test]
5170 fn precomputed_absolute_resolution_via_authoritative_matches_legacy_new() {
5171 for (temperature, pressure, altitude) in [
5172 (15.0, 1013.25, 0.0), (15.0, 1013.25, 2000.0), (-5.0, 850.0, 2000.0), (22.0, 950.0, 500.0),
5176 ] {
5177 let atmosphere = AtmosphericConditions {
5178 temperature,
5179 pressure,
5180 humidity: 50.0,
5181 altitude,
5182 };
5183 let legacy = TrajectorySolver::new(
5184 BallisticInputs::default(),
5185 WindConditions::default(),
5186 atmosphere.clone(),
5187 );
5188
5189 let (resolved_temp_c, resolved_pressure_hpa) =
5190 crate::atmosphere::resolve_station_conditions_with_pressure_mode(
5191 temperature,
5192 pressure,
5193 altitude,
5194 crate::atmosphere::PressureReferenceMode::Absolute,
5195 );
5196 let precomputed_atmosphere = AtmosphericConditions {
5197 temperature: resolved_temp_c,
5198 pressure: resolved_pressure_hpa,
5199 humidity: 50.0,
5200 altitude,
5201 };
5202 let precomputed = TrajectorySolver::new_with_resolved_station_atmosphere(
5203 BallisticInputs::default(),
5204 WindConditions::default(),
5205 precomputed_atmosphere,
5206 );
5207
5208 let (legacy_density, legacy_sos, legacy_temp_c, legacy_pressure_hpa) =
5209 legacy.resolved_atmosphere();
5210 let (pre_density, pre_sos, pre_temp_c, pre_pressure_hpa) =
5211 precomputed.resolved_atmosphere();
5212
5213 assert_eq!(
5214 legacy_temp_c.to_bits(),
5215 pre_temp_c.to_bits(),
5216 "temperature=({temperature}, {pressure}, {altitude})"
5217 );
5218 assert_eq!(
5219 legacy_pressure_hpa.to_bits(),
5220 pre_pressure_hpa.to_bits(),
5221 "pressure=({temperature}, {pressure}, {altitude})"
5222 );
5223 assert_eq!(legacy_density.to_bits(), pre_density.to_bits());
5224 assert_eq!(legacy_sos.to_bits(), pre_sos.to_bits());
5225 }
5226 }
5227
5228 fn configured_euler_zero(vertical_wind_mps: f64, time_step_s: f64) -> TrajectorySolver {
5229 let inputs = BallisticInputs {
5230 muzzle_velocity: 800.0,
5231 bc_value: 0.5,
5232 bc_type: DragModel::G7,
5233 bullet_mass: 0.0109,
5234 bullet_diameter: 0.00782,
5235 bullet_length: 0.0309,
5236 sight_height: 0.05,
5237 ground_threshold: -100.0,
5238 use_rk4: false,
5239 use_adaptive_rk45: false,
5240 ..BallisticInputs::default()
5241 };
5242 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(
5243 inputs,
5244 WindConditions::default(),
5245 AtmosphericConditions::default(),
5246 );
5247 solver.set_max_range(300.0);
5248 solver.set_time_step(time_step_s);
5249 if vertical_wind_mps != 0.0 {
5250 solver.set_wind_segments(vec![WindSegment {
5251 speed_kmh: 0.0,
5252 angle_deg: 0.0,
5253 until_m: 400.0,
5254 vertical_mps: vertical_wind_mps,
5255 }]);
5256 }
5257 solver
5258 }
5259
5260 #[test]
5261 fn inclined_shot_zeroes_like_a_level_rifle() {
5262 const ZERO_DISTANCE_M: f64 = 91.44; const SIGHT_HEIGHT_M: f64 = 0.0381; let inputs = BallisticInputs {
5271 bc_value: 0.5,
5272 bullet_mass: 150.0 * 0.06479891 / 1000.0,
5273 muzzle_velocity: 2700.0 * 0.3048,
5274 sight_height: SIGHT_HEIGHT_M,
5275 ..Default::default()
5276 };
5277
5278 let mut level = inputs.clone();
5279 level.shooting_angle = 0.0;
5280 let level_angle = TrajectorySolver::new(level, Default::default(), Default::default())
5281 .find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
5282 .expect("level zero must solve");
5283
5284 let mut inclined = inputs;
5285 inclined.shooting_angle = 5.71_f64.to_radians();
5286 let inclined_angle =
5287 TrajectorySolver::new(inclined, Default::default(), Default::default())
5288 .find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
5289 .expect("MBA-1412: a 5.71 deg incline at a 100 yd zero must be solvable");
5290
5291 assert!(
5292 (inclined_angle - level_angle).abs() < 1e-9,
5293 "zeroing is level-rifle sight geometry; incline must not move the solved zero: \
5294 level={level_angle}, inclined={inclined_angle}"
5295 );
5296 }
5297
5298 #[test]
5299 fn configured_zero_keeps_segments_method_and_time_step_then_sets_base_angle() {
5300 const TARGET_DISTANCE_M: f64 = 150.0;
5301 const TARGET_HEIGHT_M: f64 = 0.05;
5302
5303 let mut segmented = configured_euler_zero(-10.0, 0.02);
5306 let coarse_height = segmented
5307 .zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
5308 .expect("coarse configured trial")
5309 .expect("coarse trial reaches target");
5310 let mut fine = segmented.clone();
5311 fine.set_time_step(0.001);
5312 let fine_height = fine
5313 .zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
5314 .expect("fine configured trial")
5315 .expect("fine trial reaches target");
5316 assert!(
5317 (coarse_height - fine_height).abs() > 1e-5,
5318 "configured Euler step must affect zero trials: coarse={coarse_height}, fine={fine_height}"
5319 );
5320
5321 let segmented_angle = segmented
5322 .calculate_and_set_zero_angle(TARGET_DISTANCE_M, TARGET_HEIGHT_M, ZeroTargetFrame::SightLine)
5323 .expect("segmented zero");
5324 assert_eq!(
5325 segmented.inputs.muzzle_angle.to_bits(),
5326 segmented_angle.to_bits(),
5327 "successful zero must install its angle on the configured solver"
5328 );
5329 assert_eq!(segmented.time_step.to_bits(), 0.02_f64.to_bits());
5330 assert_eq!(segmented.max_range.to_bits(), 300.0_f64.to_bits());
5331 assert!(segmented.wind_sock.is_some());
5332 assert_eq!(
5333 segmented.station_atmosphere_resolution,
5334 StationAtmosphereResolution::Authoritative
5335 );
5336 let zero_height = segmented
5337 .zero_trial_height_at(segmented_angle, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
5338 .expect("verify segmented zero")
5339 .expect("zeroed trial reaches target");
5340 assert!(
5341 (zero_height - TARGET_HEIGHT_M).abs() < 0.0001,
5342 "configured zero missed target: height={zero_height}"
5343 );
5344
5345 let mut calm = configured_euler_zero(0.0, 0.02);
5346 let calm_angle = calm
5347 .calculate_and_set_zero_angle(TARGET_DISTANCE_M, TARGET_HEIGHT_M, ZeroTargetFrame::SightLine)
5348 .expect("calm zero");
5349 assert!(
5350 (segmented_angle - calm_angle).abs() > 1e-5,
5351 "segmented vertical wind must participate in zero trials: segmented={segmented_angle}, calm={calm_angle}"
5352 );
5353 }
5354}
5355
5356#[cfg(test)]
5357mod result_sanity_tests {
5358 use super::*;
5359
5360 fn default_solver() -> TrajectorySolver {
5361 TrajectorySolver::new(
5362 BallisticInputs::default(),
5363 WindConditions::default(),
5364 AtmosphericConditions::default(),
5365 )
5366 }
5367
5368 fn minimal_result() -> TrajectoryResult {
5369 TrajectoryResult {
5370 max_range: 100.0,
5371 max_height: 1.0,
5372 time_of_flight: 0.5,
5373 impact_velocity: 700.0,
5374 impact_energy: 2450.0,
5375 projectile_mass_kg: 0.01,
5376 line_of_sight_height_m: 1.5,
5377 station_speed_of_sound_mps: 340.0,
5378 termination: TrajectoryTermination::MaxRange,
5379 points: vec![],
5380 sampled_points: None,
5381 min_pitch_damping: None,
5382 transonic_mach: None,
5383 angular_state: None,
5384 max_yaw_angle: None,
5385 max_precession_angle: None,
5386 aerodynamic_jump: None,
5387 mach_1_2_distance_m: None,
5388 mach_1_0_distance_m: None,
5389 mach_0_9_distance_m: None,
5390 }
5391 }
5392
5393 #[test]
5394 fn mba1293_negative_scalars_fail_the_result_postcondition() {
5395 let solver = default_solver();
5396 solver
5397 .validate_result_sanity(&minimal_result())
5398 .expect("a sane result must pass");
5399
5400 for (name, mutate) in [
5401 ("max_range", (|r| r.max_range = -50.588) as fn(&mut TrajectoryResult)),
5402 ("time_of_flight", |r| r.time_of_flight = -1.0),
5403 ("impact_velocity", |r| r.impact_velocity = -700.0),
5404 ("impact_energy", |r| r.impact_energy = -1.0),
5405 ] {
5406 let mut result = minimal_result();
5407 mutate(&mut result);
5408 let error = solver
5409 .validate_result_sanity(&result)
5410 .expect_err("negative scalar must fail");
5411 assert!(
5412 error.to_string().contains(name),
5413 "error for {name} did not name the field: {error}"
5414 );
5415 }
5416 }
5417
5418 #[test]
5419 fn mba1293_speed_budget_bounds_legitimate_states_and_rejects_divergence() {
5420 let solver = default_solver();
5421 let mv = solver.inputs.muzzle_velocity;
5422
5423 let position = Vector3::new(10.0, 0.0, 0.0);
5425 solver
5426 .validate_integration_state(&position, &Vector3::new(mv, 0.0, 0.0), 0.01)
5427 .expect("muzzle-speed state must pass");
5428
5429 let error = solver
5431 .validate_integration_state(&position, &Vector3::new(-13.0 * mv, 0.0, 0.0), 0.01)
5432 .expect_err("13x muzzle speed must fail the budget");
5433 assert!(error.to_string().contains("diverged"), "{error}");
5434
5435 let after_fall = mv + crate::constants::G_ACCEL_MPS2 * 60.0;
5437 solver
5438 .validate_integration_state(&position, &Vector3::new(0.0, -after_fall, 0.0), 60.0)
5439 .expect("gravity-accelerated speed within g*t must pass");
5440 }
5441}
5442
5443#[cfg(test)]
5444mod trajectory_point_budget_tests {
5445 use super::*;
5446 use crate::MAX_TRAJECTORY_SAMPLES;
5447
5448 fn solver_with_budget(
5449 use_rk4: bool,
5450 use_adaptive_rk45: bool,
5451 point_budget: usize,
5452 max_range: f64,
5453 ) -> TrajectorySolver {
5454 let inputs = BallisticInputs {
5455 use_rk4,
5456 use_adaptive_rk45,
5457 ground_threshold: f64::NEG_INFINITY,
5458 ..BallisticInputs::default()
5459 };
5460 let mut solver = TrajectorySolver::new(
5461 inputs,
5462 WindConditions::default(),
5463 AtmosphericConditions::default(),
5464 );
5465 solver.max_trajectory_points = point_budget;
5466 solver.set_max_range(max_range);
5467 solver.set_time_step(0.001);
5468 solver
5469 }
5470
5471 #[test]
5472 fn mba1283_every_solver_errors_instead_of_exceeding_point_budget() {
5473 for (mode, use_rk4, use_adaptive_rk45) in [
5474 ("Euler", false, false),
5475 ("RK4", true, false),
5476 ("RK45", true, true),
5477 ] {
5478 let error = solver_with_budget(use_rk4, use_adaptive_rk45, 3, 10.0)
5479 .solve()
5480 .expect_err("a solve requiring more than three points must fail");
5481 assert!(
5482 error.to_string().contains("point limit of 3"),
5483 "unexpected {mode} point-budget error: {error}"
5484 );
5485 }
5486 }
5487
5488 #[test]
5489 fn mba1283_interpolated_endpoint_counts_toward_point_budget() {
5490 for (mode, use_rk4, use_adaptive_rk45) in [
5491 ("Euler", false, false),
5492 ("RK4", true, false),
5493 ("RK45", true, true),
5494 ] {
5495 let result = solver_with_budget(use_rk4, use_adaptive_rk45, 2, 0.1)
5496 .solve()
5497 .expect("the initial point plus exact endpoint fit a two-point budget");
5498 assert_eq!(result.points.len(), 2, "unexpected {mode} point count");
5499
5500 let error = solver_with_budget(use_rk4, use_adaptive_rk45, 1, 0.1)
5501 .solve()
5502 .expect_err("the exact endpoint must not exceed a one-point budget");
5503 assert!(
5504 error.to_string().contains("point limit of 1"),
5505 "unexpected {mode} endpoint-budget error: {error}"
5506 );
5507 }
5508 }
5509
5510 #[test]
5511 fn mba1299_every_solver_preflights_the_sample_budget() {
5512 for (mode, use_rk4, use_adaptive_rk45) in [
5513 ("Euler", false, false),
5514 ("RK4", true, false),
5515 ("RK45", true, true),
5516 ] {
5517 let inputs = BallisticInputs {
5518 use_rk4,
5519 use_adaptive_rk45,
5520 enable_trajectory_sampling: true,
5521 sample_interval: 1.0,
5522 ground_threshold: f64::NEG_INFINITY,
5523 ..BallisticInputs::default()
5524 };
5525 let mut solver = TrajectorySolver::new(
5526 inputs,
5527 WindConditions::default(),
5528 AtmosphericConditions::default(),
5529 );
5530 solver.set_max_range(MAX_TRAJECTORY_SAMPLES as f64);
5531 solver.max_trajectory_points = 0;
5534
5535 let error = solver
5536 .solve()
5537 .expect_err("an over-limit sample grid must fail before integration");
5538 assert!(
5539 error
5540 .to_string()
5541 .contains("trajectory sample limit of 250000 exceeded"),
5542 "unexpected {mode} sample-budget error: {error}"
5543 );
5544 }
5545 }
5546
5547 #[test]
5548 fn mba1299_normal_sampling_does_not_change_solver_results() {
5549 for (mode, use_rk4, use_adaptive_rk45) in [
5550 ("Euler", false, false),
5551 ("RK4", true, false),
5552 ("RK45", true, true),
5553 ] {
5554 let solve = |enable_trajectory_sampling| {
5555 let inputs = BallisticInputs {
5556 use_rk4,
5557 use_adaptive_rk45,
5558 enable_trajectory_sampling,
5559 sample_interval: 0.5,
5560 ground_threshold: f64::NEG_INFINITY,
5561 ..BallisticInputs::default()
5562 };
5563 let mut solver = TrajectorySolver::new(
5564 inputs,
5565 WindConditions::default(),
5566 AtmosphericConditions::default(),
5567 );
5568 solver.set_max_range(2.0);
5569 solver.solve().expect("normal short-range solve")
5570 };
5571
5572 let baseline = solve(false);
5573 let sampled = solve(true);
5574 for (field, left, right) in [
5575 ("max_range", baseline.max_range, sampled.max_range),
5576 ("max_height", baseline.max_height, sampled.max_height),
5577 (
5578 "time_of_flight",
5579 baseline.time_of_flight,
5580 sampled.time_of_flight,
5581 ),
5582 (
5583 "impact_velocity",
5584 baseline.impact_velocity,
5585 sampled.impact_velocity,
5586 ),
5587 (
5588 "impact_energy",
5589 baseline.impact_energy,
5590 sampled.impact_energy,
5591 ),
5592 ] {
5593 assert_eq!(
5594 left.to_bits(),
5595 right.to_bits(),
5596 "{mode} sampling changed {field}"
5597 );
5598 }
5599 assert_eq!(baseline.points.len(), sampled.points.len());
5600 for (index, (left, right)) in baseline
5601 .points
5602 .iter()
5603 .zip(&sampled.points)
5604 .enumerate()
5605 {
5606 assert_eq!(left.time.to_bits(), right.time.to_bits(), "{mode} point {index}");
5607 assert_eq!(
5608 left.position.map(f64::to_bits),
5609 right.position.map(f64::to_bits),
5610 "{mode} point {index} position"
5611 );
5612 assert_eq!(
5613 left.velocity_magnitude.to_bits(),
5614 right.velocity_magnitude.to_bits(),
5615 "{mode} point {index} velocity"
5616 );
5617 assert_eq!(
5618 left.kinetic_energy.to_bits(),
5619 right.kinetic_energy.to_bits(),
5620 "{mode} point {index} energy"
5621 );
5622 }
5623 assert!(baseline.sampled_points.is_none());
5624 let samples = sampled
5625 .sampled_points
5626 .expect("sampling-enabled solve should return observations");
5627 assert_eq!(
5628 samples
5629 .iter()
5630 .map(|sample| sample.distance_m)
5631 .collect::<Vec<_>>(),
5632 vec![0.0, 0.5, 1.0, 1.5, 2.0],
5633 "{mode} normal sampling grid changed"
5634 );
5635 }
5636 }
5637}
5638
5639#[cfg(test)]
5640mod monte_carlo_result_tests {
5641 use super::*;
5642
5643 fn make_results(impact_positions: Vec<Vector3<f64>>) -> MonteCarloResults {
5644 let count = impact_positions.len();
5645 MonteCarloResults {
5646 ranges: vec![500.0; count],
5647 impact_velocities: vec![300.0; count],
5648 impact_positions,
5649 }
5650 }
5651
5652 #[test]
5653 fn target_plane_cep_excludes_shortfall_markers() {
5654 let mut positions: Vec<Vector3<f64>> = (1..=5)
5655 .map(|radius| Vector3::new(0.0, radius as f64, 0.0))
5656 .collect();
5657 positions.extend(
5658 (0..5).map(|_| Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)),
5659 );
5660 let results = make_results(positions);
5661
5662 assert_eq!(results.target_arrival_count(), 5);
5663 assert_eq!(results.target_shortfall_fraction(), 0.5);
5664 assert_eq!(results.target_plane_cep(), Some(3.0));
5665
5666 let one_shortfall = make_results(vec![
5667 Vector3::new(0.0, 1.0, 0.0),
5668 Vector3::new(0.0, 2.0, 0.0),
5669 Vector3::new(0.0, 3.0, 0.0),
5670 Vector3::new(0.0, 4.0, 0.0),
5671 Vector3::new(0.0, 5.0, 0.0),
5672 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5673 ]);
5674 assert_eq!(one_shortfall.target_plane_cep(), Some(3.0));
5675 }
5676
5677 #[test]
5678 fn all_shortfalls_have_no_cep_but_still_count_as_misses() {
5679 let all_shortfalls = make_results(vec![
5680 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5681 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5682 ]);
5683 assert_eq!(all_shortfalls.target_arrival_count(), 0);
5684 assert_eq!(all_shortfalls.target_shortfall_fraction(), 1.0);
5685 assert_eq!(all_shortfalls.target_plane_cep(), None);
5686 assert_eq!(all_shortfalls.hit_probability(0.3), 0.0);
5687
5688 let one_hit_one_shortfall = make_results(vec![
5689 Vector3::new(0.0, 0.1, 0.0),
5690 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5691 ]);
5692 assert_eq!(one_hit_one_shortfall.hit_probability(0.3), 0.5);
5693 }
5694
5695 #[test]
5697 fn rect_hit_probability_checks_independent_axis_halves() {
5698 let results = make_results(vec![
5699 Vector3::new(0.0, 0.1, 0.1),
5701 Vector3::new(0.0, 0.0, 0.2),
5703 Vector3::new(0.0, 0.0, 0.201),
5705 Vector3::new(0.0, 0.301, 0.0),
5707 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5709 ]);
5710 assert!((results.rect_hit_probability(0.4, 0.6) - 0.4).abs() < 1e-12);
5712 }
5713
5714 #[test]
5715 fn rect_hit_probability_matches_circular_hit_probability_for_a_centered_hit() {
5716 let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
5717 assert_eq!(results.rect_hit_probability(0.5, 0.5), 1.0);
5718 assert_eq!(results.hit_probability(0.3), 1.0);
5719 }
5720
5721 #[test]
5722 fn rect_hit_probability_is_zero_for_empty_or_nonpositive_dimensions() {
5723 let empty = make_results(vec![]);
5724 assert_eq!(empty.rect_hit_probability(1.0, 1.0), 0.0);
5725
5726 let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
5727 assert_eq!(results.rect_hit_probability(0.0, 1.0), 0.0);
5728 assert_eq!(results.rect_hit_probability(1.0, 0.0), 0.0);
5729 assert_eq!(results.rect_hit_probability(-1.0, 1.0), 0.0);
5730 }
5731}
5732
5733#[cfg(test)]
5734mod monte_carlo_seeded_tests {
5735 use super::*;
5736
5737 fn seeded_test_fixture() -> (BallisticInputs, WindConditions) {
5746 (
5747 BallisticInputs {
5748 muzzle_velocity: 800.0,
5749 ..BallisticInputs::default()
5750 },
5751 WindConditions::default(),
5752 )
5753 }
5754
5755 fn loose_params() -> MonteCarloParams {
5765 MonteCarloParams {
5766 num_simulations: 1, velocity_std_dev: 3.0,
5768 angle_std_dev: 3.5e-4,
5769 bc_std_dev: 0.01,
5770 wind_speed_std_dev: 1.0,
5771 target_distance: Some(300.0),
5772 base_wind_speed: 0.0,
5773 base_wind_direction: 0.0,
5774 azimuth_std_dev: 3.5e-4,
5775 }
5776 }
5777
5778 fn mixed_arrival_params() -> MonteCarloParams {
5786 MonteCarloParams {
5787 num_simulations: 1,
5788 target_distance: Some(1920.0),
5789 ..MonteCarloParams::default()
5790 }
5791 }
5792
5793 #[test]
5809 fn legacy_seeded_estimates_are_pinned_bit_for_bit() {
5810 let (inputs, wind) = seeded_test_fixture();
5811 let params = MonteCarloParams {
5812 num_simulations: 200,
5813 target_distance: Some(500.0),
5814 ..MonteCarloParams::default()
5815 };
5816
5817 let results = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5818 inputs,
5819 wind,
5820 params,
5821 0.01,
5822 0x1352_5EED,
5823 )
5824 .expect("seeded legacy run");
5825
5826 assert_eq!(results.ranges.len(), 200, "ranges length");
5831 assert_eq!(results.impact_velocities.len(), 200, "impact_velocities length");
5832 assert_eq!(results.impact_positions.len(), 200, "impact_positions length");
5833
5834 assert_eq!(
5836 results.hit_probability(DEFAULT_HIT_RADIUS_M).to_bits(),
5837 0.14_f64.to_bits(),
5838 "hit_probability = {:?}",
5839 results.hit_probability(DEFAULT_HIT_RADIUS_M)
5840 );
5841
5842 let expected_ranges: [f64; 3] =
5843 [1907.972891143359, 1936.408435469319, 1912.8150447617645];
5844 let expected_velocities: [f64; 3] =
5845 [238.6187151542299, 239.91923651600106, 243.14112455427164];
5846 let expected_positions: [(f64, f64, f64); 3] = [
5847 (0.0, -0.0643556039548101, 0.7344970252014579),
5848 (0.0, 0.5769422971539162, 0.27227201756386726),
5849 (0.0, -0.7440425792472842, 0.1541804446282822),
5850 ];
5851
5852 for (i, expected) in expected_ranges.iter().enumerate() {
5853 assert_eq!(
5854 results.ranges[i].to_bits(),
5855 expected.to_bits(),
5856 "ranges[{i}] = {:?}, pinned {expected:?}",
5857 results.ranges[i]
5858 );
5859 }
5860 for (i, expected) in expected_velocities.iter().enumerate() {
5861 assert_eq!(
5862 results.impact_velocities[i].to_bits(),
5863 expected.to_bits(),
5864 "impact_velocities[{i}] = {:?}, pinned {expected:?}",
5865 results.impact_velocities[i]
5866 );
5867 }
5868 for (i, (x, y, z)) in expected_positions.iter().enumerate() {
5869 let actual = results.impact_positions[i];
5870 assert_eq!(actual.x.to_bits(), x.to_bits(), "impact_positions[{i}].x = {:?}", actual.x);
5871 assert_eq!(actual.y.to_bits(), y.to_bits(), "impact_positions[{i}].y = {:?}", actual.y);
5872 assert_eq!(actual.z.to_bits(), z.to_bits(), "impact_positions[{i}].z = {:?}", actual.z);
5873 }
5874 }
5875
5876 #[test]
5877 fn seeded_runs_are_deterministic_and_match_the_using_rng_path() {
5878 let inputs = BallisticInputs {
5879 muzzle_velocity: 800.0,
5880 ..BallisticInputs::default()
5881 };
5882 let params = MonteCarloParams {
5883 num_simulations: 64,
5884 target_distance: Some(200.0),
5885 ..MonteCarloParams::default()
5886 };
5887
5888 let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5889 inputs.clone(),
5890 WindConditions::default(),
5891 params.clone(),
5892 0.01,
5893 42,
5894 )
5895 .expect("seeded run a");
5896 let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5897 inputs,
5898 WindConditions::default(),
5899 params,
5900 0.01,
5901 42,
5902 )
5903 .expect("seeded run b");
5904
5905 assert_eq!(a.ranges.len(), b.ranges.len());
5906 for (ra, rb) in a.ranges.iter().zip(b.ranges.iter()) {
5907 assert_eq!(ra.to_bits(), rb.to_bits());
5908 }
5909 for (pa, pb) in a.impact_positions.iter().zip(b.impact_positions.iter()) {
5910 assert_eq!(pa.x.to_bits(), pb.x.to_bits());
5911 assert_eq!(pa.y.to_bits(), pb.y.to_bits());
5912 assert_eq!(pa.z.to_bits(), pb.z.to_bits());
5913 }
5914 }
5915
5916 #[test]
5917 fn different_seeds_generally_produce_different_draws() {
5918 let inputs = BallisticInputs {
5919 muzzle_velocity: 800.0,
5920 ..BallisticInputs::default()
5921 };
5922 let params = MonteCarloParams {
5923 num_simulations: 32,
5924 velocity_std_dev: 5.0,
5925 target_distance: Some(200.0),
5926 ..MonteCarloParams::default()
5927 };
5928
5929 let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5930 inputs.clone(),
5931 WindConditions::default(),
5932 params.clone(),
5933 0.0,
5934 1,
5935 )
5936 .expect("seeded run a");
5937 let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5938 inputs,
5939 WindConditions::default(),
5940 params,
5941 0.0,
5942 2,
5943 )
5944 .expect("seeded run b");
5945
5946 assert_ne!(a.impact_velocities, b.impact_velocities);
5947 }
5948
5949 #[test]
5950 fn adaptive_stops_at_target_half_width_on_an_easy_case() {
5951 let (inputs, wind) = seeded_test_fixture();
5952 let conv = McConvergence {
5953 target_half_width: 0.05,
5954 ..Default::default()
5955 };
5956 let r = run_monte_carlo_adaptive_seeded(
5957 &inputs,
5958 &wind,
5959 &loose_params(),
5960 &conv,
5961 DEFAULT_HIT_RADIUS_M,
5962 0x1352_ADA9,
5963 )
5964 .unwrap();
5965
5966 assert_eq!(r.stop_reason, McStopReason::TargetHalfWidthMet);
5967 assert!(
5968 (r.ci_high - r.ci_low) / 2.0 <= 0.05 + 1e-12,
5969 "half-width {} exceeds the requested 0.05",
5970 (r.ci_high - r.ci_low) / 2.0
5971 );
5972 assert!(r.samples >= conv.min_samples, "stopped below min_samples");
5973 assert!(r.samples < conv.max_samples, "did not actually stop early");
5974 assert!(r.samples.is_multiple_of(conv.batch_size) || r.samples == conv.min_samples);
5975 assert!(r.ci_low <= r.hit_probability && r.hit_probability <= r.ci_high);
5976
5977 assert_eq!(r.hit_radius_m, DEFAULT_HIT_RADIUS_M);
5979 assert_eq!(r.target_distance_m, 300.0);
5980 assert_eq!(r.confidence_percent, 95);
5981 assert!(
5983 r.samples > 1,
5984 "params.num_simulations must be ignored by the adaptive driver"
5985 );
5986 assert!(
5989 r.mean_impact_velocity_mps > 0.0,
5990 "no impact velocity accumulated"
5991 );
5992 assert!(
5993 r.std_drop_at_target_m > 0.0 && r.std_wind_drift_at_target_m > 0.0,
5994 "dispersion collapsed: drop sd {} drift sd {}",
5995 r.std_drop_at_target_m,
5996 r.std_wind_drift_at_target_m
5997 );
5998 }
5999
6000 #[test]
6001 fn adaptive_caps_at_max_samples_on_an_impossible_target() {
6002 let (inputs, wind) = seeded_test_fixture();
6003 let conv = McConvergence {
6004 target_half_width: 1e-6,
6005 max_samples: 3_000,
6006 batch_size: 500,
6007 min_samples: 1_000,
6008 level: ConfidenceLevel::P95,
6009 };
6010 let r = run_monte_carlo_adaptive_seeded(
6011 &inputs,
6012 &wind,
6013 &loose_params(),
6014 &conv,
6015 DEFAULT_HIT_RADIUS_M,
6016 7,
6017 )
6018 .unwrap();
6019
6020 assert_eq!(r.stop_reason, McStopReason::MaxSamplesReached);
6021 assert_eq!(r.samples, 3_000);
6022 }
6023
6024 #[test]
6033 fn adaptive_stops_between_the_floor_and_the_ceiling() {
6034 let (inputs, wind) = seeded_test_fixture();
6035 let conv = McConvergence {
6036 level: ConfidenceLevel::P95,
6037 target_half_width: 0.03,
6038 min_samples: 0,
6039 max_samples: 10_000,
6040 batch_size: 100,
6041 };
6042 let r = run_monte_carlo_adaptive_seeded(
6043 &inputs,
6044 &wind,
6045 &loose_params(),
6046 &conv,
6047 DEFAULT_HIT_RADIUS_M,
6048 0x1352_5A1D,
6049 )
6050 .unwrap();
6051
6052 assert_eq!(r.stop_reason, McStopReason::TargetHalfWidthMet);
6053 assert!(r.samples > 0);
6054 assert!(
6055 r.samples.is_multiple_of(conv.batch_size),
6056 "samples {} is not a whole number of batches",
6057 r.samples
6058 );
6059 assert!(
6060 r.samples > conv.min_samples,
6061 "stopped on the floor, not on the data"
6062 );
6063 assert!(
6064 r.samples < conv.max_samples,
6065 "ran to the ceiling, so nothing adaptive was exercised"
6066 );
6067 assert!(
6068 r.samples > conv.batch_size,
6069 "stopped on the very first look ({} samples); the multi-batch path is untested",
6070 r.samples
6071 );
6072 assert!((r.ci_high - r.ci_low) / 2.0 <= 0.03 + 1e-12);
6073 assert!(r.ci_low <= r.hit_probability && r.hit_probability <= r.ci_high);
6074 assert_eq!(r.attempts, r.samples, "no trial should have been dropped");
6075 }
6076
6077 #[test]
6085 fn adaptive_runs_a_truncated_final_batch_up_to_max_samples() {
6086 let (inputs, wind) = seeded_test_fixture();
6087 let conv = McConvergence {
6088 level: ConfidenceLevel::P95,
6089 target_half_width: 1e-6,
6090 min_samples: 0,
6091 max_samples: 750,
6092 batch_size: 500,
6093 };
6094 let r = run_monte_carlo_adaptive_seeded(
6095 &inputs,
6096 &wind,
6097 &loose_params(),
6098 &conv,
6099 DEFAULT_HIT_RADIUS_M,
6100 0x1352_7B10,
6101 )
6102 .unwrap();
6103
6104 assert_eq!(r.stop_reason, McStopReason::MaxSamplesReached);
6105 assert_eq!(
6106 r.samples, 750,
6107 "the 250-trial final batch did not run, or was not truncated"
6108 );
6109 assert_eq!(r.attempts, 750);
6110 assert!(!r.samples.is_multiple_of(conv.batch_size));
6111 }
6112
6113 #[test]
6114 fn adaptive_is_deterministic_for_a_seed() {
6115 let (inputs, wind) = seeded_test_fixture();
6116 let conv = McConvergence::default();
6117 let a = run_monte_carlo_adaptive_seeded(
6118 &inputs,
6119 &wind,
6120 &loose_params(),
6121 &conv,
6122 DEFAULT_HIT_RADIUS_M,
6123 99,
6124 )
6125 .unwrap();
6126 let b = run_monte_carlo_adaptive_seeded(
6127 &inputs,
6128 &wind,
6129 &loose_params(),
6130 &conv,
6131 DEFAULT_HIT_RADIUS_M,
6132 99,
6133 )
6134 .unwrap();
6135
6136 assert_eq!(a.hit_probability.to_bits(), b.hit_probability.to_bits());
6137 assert_eq!(a.samples, b.samples);
6138 assert_eq!(a.ci_low.to_bits(), b.ci_low.to_bits());
6139 assert_eq!(a.ci_high.to_bits(), b.ci_high.to_bits());
6140 assert_eq!(
6143 a.mean_impact_velocity_mps.to_bits(),
6144 b.mean_impact_velocity_mps.to_bits()
6145 );
6146 assert_eq!(
6147 a.std_drop_at_target_m.to_bits(),
6148 b.std_drop_at_target_m.to_bits()
6149 );
6150 }
6151
6152 #[test]
6153 fn adaptive_report_carries_schema_method_and_all_four_assumptions() {
6154 let (inputs, wind) = seeded_test_fixture();
6155 let conv = McConvergence {
6159 min_samples: 0,
6160 max_samples: 50,
6161 batch_size: 50,
6162 target_half_width: 1.0,
6163 level: ConfidenceLevel::P90,
6164 };
6165 let r = run_monte_carlo_adaptive_seeded(
6166 &inputs,
6167 &wind,
6168 &mixed_arrival_params(),
6169 &conv,
6170 DEFAULT_HIT_RADIUS_M,
6171 0x1352_D0C5,
6172 )
6173 .unwrap();
6174
6175 assert_eq!(r.schema_version, MC_ADAPTIVE_SCHEMA_VERSION_V1);
6176 assert_eq!(r.schema_version, 1);
6177 assert_eq!(r.method, "anytime_beta_binomial_mixture_cs_v1");
6178 assert_eq!(r.confidence_percent, 90);
6179
6180 assert_eq!(r.attempts, 50, "one full batch was drawn");
6186 assert_eq!(r.samples, 50, "no trial was dropped by the solver");
6187 assert!(
6188 r.arrivals > 0 && r.arrivals < r.samples,
6189 "fixture must split the run: arrivals {} of samples {}",
6190 r.arrivals,
6191 r.samples
6192 );
6193 assert!(r.arrivals >= 2, "arrivals {} too few for a sample sd", r.arrivals);
6196 assert!(r.std_drop_at_target_m > 0.0 && r.std_impact_velocity_mps > 0.0);
6197 assert!(r.attempts >= r.samples && r.samples >= r.arrivals);
6199
6200 assert_eq!(r.assumptions.len(), 4, "exactly four assumptions expected");
6206 assert_eq!(
6207 r.assumptions[0],
6208 "Sampling uncertainty only: intervals cover Monte Carlo sampling error, not model error in the trajectory solver or its inputs."
6209 );
6210 assert_eq!(
6211 r.assumptions[1],
6212 "Anytime-valid stopping: the beta-binomial mixture confidence sequence keeps its coverage guarantee despite stopping the moment the target half-width is met."
6213 );
6214 assert_eq!(
6215 r.assumptions[2],
6216 "Input dispersions are the independent normal distributions declared in MonteCarloParams; correlations between inputs are not modeled."
6217 );
6218 assert_eq!(
6219 r.assumptions[3],
6220 "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."
6221 );
6222
6223 assert_eq!(
6225 serde_json::to_string(&McStopReason::TargetHalfWidthMet).unwrap(),
6226 "\"target_half_width_met\""
6227 );
6228 assert_eq!(
6229 serde_json::to_string(&McStopReason::MaxSamplesReached).unwrap(),
6230 "\"max_samples_reached\""
6231 );
6232 }
6233
6234 #[test]
6235 fn adaptive_rejects_nonsense_convergence() {
6236 let (inputs, wind) = seeded_test_fixture();
6237 let run = |conv: McConvergence| {
6238 run_monte_carlo_adaptive_seeded(
6239 &inputs,
6240 &wind,
6241 &loose_params(),
6242 &conv,
6243 DEFAULT_HIT_RADIUS_M,
6244 1,
6245 )
6246 .unwrap_err()
6247 };
6248
6249 for bad_width in [0.0, -0.01, f64::NAN] {
6250 let err = run(McConvergence {
6251 target_half_width: bad_width,
6252 ..Default::default()
6253 });
6254 assert!(
6255 err.contains("target_half_width"),
6256 "error must name the field, got: {err}"
6257 );
6258 }
6259
6260 let err = run(McConvergence {
6261 batch_size: 0,
6262 ..Default::default()
6263 });
6264 assert!(err.contains("batch_size"), "got: {err}");
6265
6266 let err = run(McConvergence {
6267 min_samples: 5_000,
6268 max_samples: 1_000,
6269 ..Default::default()
6270 });
6271 assert!(err.contains("max_samples"), "got: {err}");
6272 assert!(err.contains("min_samples"), "got: {err}");
6273
6274 let err = run(McConvergence {
6275 max_samples: 0,
6276 min_samples: 0,
6277 ..Default::default()
6278 });
6279 assert!(err.contains("max_samples"), "got: {err}");
6280
6281 let err = McConvergence {
6283 batch_size: 0,
6284 ..Default::default()
6285 }
6286 .validate()
6287 .unwrap_err();
6288 assert!(err.contains("batch_size"));
6289 }
6290
6291 #[test]
6292 fn wilson_companion_matches_hit_probability_and_wilson_interval() {
6293 let (inputs, wind) = seeded_test_fixture();
6294 let params = MonteCarloParams {
6295 num_simulations: 128,
6296 target_distance: Some(500.0),
6297 ..MonteCarloParams::default()
6298 };
6299 let results = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
6300 inputs,
6301 wind,
6302 params,
6303 0.01,
6304 0x1352_C0DE,
6305 )
6306 .expect("seeded legacy run");
6307
6308 for level in [
6309 ConfidenceLevel::P90,
6310 ConfidenceLevel::P95,
6311 ConfidenceLevel::P99,
6312 ] {
6313 let (p_hat, (lo, hi), n) =
6314 results.hit_probability_wilson(DEFAULT_HIT_RADIUS_M, level);
6315
6316 assert_eq!(
6319 p_hat.to_bits(),
6320 results.hit_probability(DEFAULT_HIT_RADIUS_M).to_bits()
6321 );
6322 assert_eq!(n, results.impact_positions.len() as u64);
6323
6324 let hits = results
6328 .impact_positions
6329 .iter()
6330 .filter(|p| MonteCarloResults::position_is_hit(p, DEFAULT_HIT_RADIUS_M))
6331 .count() as u64;
6332 let (want_lo, want_hi) = wilson_interval(hits, n, level);
6333 assert_eq!(lo.to_bits(), want_lo.to_bits());
6334 assert_eq!(hi.to_bits(), want_hi.to_bits());
6335 assert!(lo <= p_hat && p_hat <= hi, "interval excludes p_hat");
6336 }
6337
6338 let empty = MonteCarloResults {
6340 ranges: Vec::new(),
6341 impact_velocities: Vec::new(),
6342 impact_positions: Vec::new(),
6343 };
6344 assert_eq!(
6345 empty.hit_probability_wilson(DEFAULT_HIT_RADIUS_M, ConfidenceLevel::P95),
6346 (0.0, (0.0, 1.0), 0)
6347 );
6348 }
6349}
6350
6351#[cfg(test)]
6352mod monte_carlo_powder_curve_tests {
6353 use super::*;
6354 use rand::{rngs::StdRng, SeedableRng};
6355
6356 #[test]
6357 fn powder_curve_preserves_sampled_muzzle_velocity_dispersion() {
6358 let inputs = BallisticInputs {
6359 muzzle_velocity: 700.0,
6360 powder_temp_curve: Some(vec![(15.0, 800.0)]),
6361 powder_curve_temp_c: Some(15.0),
6362 ..BallisticInputs::default()
6363 };
6364 let params = MonteCarloParams {
6365 num_simulations: 16,
6366 velocity_std_dev: 20.0,
6367 angle_std_dev: 1e-12,
6368 bc_std_dev: 1e-12,
6369 wind_speed_std_dev: 1e-12,
6370 target_distance: Some(100.0),
6371 azimuth_std_dev: 1e-12,
6372 ..MonteCarloParams::default()
6373 };
6374
6375 let mut rng = StdRng::seed_from_u64(0x5EED_1176);
6376 let results = run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
6377 inputs,
6378 WindConditions::default(),
6379 params,
6380 0.0,
6381 &mut rng,
6382 )
6383 .expect("Monte Carlo solve");
6384 let min_velocity = results
6385 .impact_velocities
6386 .iter()
6387 .copied()
6388 .fold(f64::INFINITY, f64::min);
6389 let max_velocity = results
6390 .impact_velocities
6391 .iter()
6392 .copied()
6393 .fold(f64::NEG_INFINITY, f64::max);
6394
6395 assert!(
6396 max_velocity - min_velocity > 1.0,
6397 "20 m/s muzzle spread collapsed after curve resolution: impact-velocity span={} m/s",
6398 max_velocity - min_velocity
6399 );
6400 }
6401}
6402
6403#[cfg(test)]
6404mod monte_carlo_wind_sampling_tests {
6405 use super::*;
6406 use rand::{rngs::StdRng, SeedableRng};
6407
6408 #[test]
6409 fn wind_speed_sigma_does_not_change_seeded_direction_draws() {
6410 let base_wind = WindConditions {
6411 speed: 100.0,
6412 direction: 0.37,
6413 vertical_speed: 0.0,
6414 };
6415 let narrow_speed = MonteCarloWindSampler::new(&base_wind, 0.5, 0.2).unwrap();
6416 let wide_speed = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
6417 let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
6418 let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
6419 let mut speed_changed = false;
6420
6421 for _ in 0..32 {
6422 let narrow = narrow_speed.sample(&mut narrow_rng);
6423 let wide = wide_speed.sample(&mut wide_rng);
6424 assert!(narrow.speed > 0.0 && wide.speed > 0.0);
6425 assert_eq!(narrow.direction.to_bits(), wide.direction.to_bits());
6426 speed_changed |= narrow.speed.to_bits() != wide.speed.to_bits();
6427 }
6428 assert!(
6429 speed_changed,
6430 "different speed sigmas must still vary speed draws"
6431 );
6432 }
6433
6434 #[test]
6435 fn zero_direction_sigma_has_no_angular_jitter() {
6436 let base_wind = WindConditions {
6437 speed: 100.0,
6438 direction: 0.37,
6439 vertical_speed: 0.0,
6440 };
6441 let sampler = MonteCarloWindSampler::new(&base_wind, 4.0, 0.0).unwrap();
6442 let mut rng = StdRng::seed_from_u64(0x5EED_1223);
6443 let mut speed_changed = false;
6444
6445 for _ in 0..32 {
6446 let wind = sampler.sample(&mut rng);
6447 speed_changed |= wind.speed.to_bits() != base_wind.speed.to_bits();
6448 assert_eq!(wind.direction.to_bits(), base_wind.direction.to_bits());
6449 }
6450 assert!(speed_changed, "speed uncertainty should remain active");
6451 }
6452
6453 #[test]
6454 fn direction_sigma_controls_seeded_angular_spread_in_radians() {
6455 let base_wind = WindConditions {
6456 speed: 100.0,
6457 direction: 0.37,
6458 vertical_speed: 0.0,
6459 };
6460 let narrow = MonteCarloWindSampler::new(&base_wind, 4.0, 0.1).unwrap();
6461 let wide = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
6462 let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
6463 let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
6464 let mut nonzero_direction_draw = false;
6465
6466 for _ in 0..32 {
6467 let narrow_wind = narrow.sample(&mut narrow_rng);
6468 let wide_wind = wide.sample(&mut wide_rng);
6469 assert_eq!(narrow_wind.speed.to_bits(), wide_wind.speed.to_bits());
6470
6471 let narrow_delta = narrow_wind.direction - base_wind.direction;
6472 let wide_delta = wide_wind.direction - base_wind.direction;
6473 assert!((wide_delta - 2.0 * narrow_delta).abs() < 1e-12);
6474 nonzero_direction_draw |= narrow_delta.abs() > 1e-6;
6475 }
6476 assert!(
6477 nonzero_direction_draw,
6478 "positive radians sigma must vary direction"
6479 );
6480 }
6481
6482 #[test]
6483 fn direction_sigma_rejects_negative_or_nonfinite_values() {
6484 let base_wind = WindConditions::default();
6485 for sigma in [-0.1, f64::NAN, f64::INFINITY] {
6486 assert!(MonteCarloWindSampler::new(&base_wind, 1.0, sigma).is_err());
6487 }
6488 }
6489
6490 #[test]
6491 fn base_vertical_wind_rides_into_every_mc_sample() {
6492 use rand::SeedableRng;
6496 let base_wind = WindConditions { vertical_speed: 4.2, ..Default::default() };
6497 let sampler = MonteCarloWindSampler::new(&base_wind, 1.0, 0.2).unwrap();
6498 let mut rng = rand::rngs::StdRng::seed_from_u64(7);
6499 for _ in 0..32 {
6500 let w = sampler.sample(&mut rng);
6501 assert_eq!(w.vertical_speed, 4.2);
6502 }
6503 }
6504
6505 #[test]
6506 fn negative_speed_sample_reverses_wind_direction() {
6507 let direction = 0.25;
6508 let signed_speed = -2.5;
6509 let wind = wind_from_signed_speed_sample(signed_speed, direction, 0.0);
6510 let positive_wind = wind_from_signed_speed_sample(2.5, direction, 0.0);
6511
6512 assert_eq!(wind.speed, 2.5);
6513 assert!(
6514 (wind.direction - (direction + std::f64::consts::PI)).abs() < f64::EPSILON,
6515 "negative speed must reverse direction by pi: got {}",
6516 wind.direction
6517 );
6518 assert_eq!(positive_wind.speed, 2.5);
6519 assert_eq!(positive_wind.direction, direction);
6520
6521 let normalized_x = -wind.speed * wind.direction.cos();
6522 let normalized_z = -wind.speed * wind.direction.sin();
6523 let signed_x = -signed_speed * direction.cos();
6524 let signed_z = -signed_speed * direction.sin();
6525 assert!((normalized_x - signed_x).abs() < 1e-12);
6526 assert!((normalized_z - signed_z).abs() < 1e-12);
6527 }
6528}
6529
6530#[cfg(test)]
6531mod bc_fit_objective_tests {
6532 use super::*;
6533
6534 fn velocity_point(range_m: f64, velocity_mps: f64) -> TrajectoryPoint {
6535 TrajectoryPoint {
6536 time: 0.0,
6537 position: Vector3::new(range_m, 0.0, 0.0),
6538 velocity_magnitude: velocity_mps,
6539 kinetic_energy: 0.0,
6540 drag_coefficient: None,
6541 }
6542 }
6543
6544 #[test]
6545 fn candidate_that_misses_an_observation_has_no_score() {
6546 let trajectory = vec![velocity_point(0.0, 800.0), velocity_point(100.0, 700.0)];
6547 let observations = vec![(50.0, 750.0), (150.0, 600.0)];
6548
6549 assert!(
6550 fit_residual_sse(&trajectory, &observations, BcFitMode::Velocity, 0.0).is_none(),
6551 "a candidate that reaches only one of two observations must not compete on partial SSE"
6552 );
6553
6554 let complete_observations = vec![(50.0, 740.0), (100.0, 680.0)];
6555 assert_eq!(
6556 fit_residual_sse(
6557 &trajectory,
6558 &complete_observations,
6559 BcFitMode::Velocity,
6560 0.0,
6561 ),
6562 Some(500.0)
6563 );
6564 }
6565}
6566
6567#[cfg(test)]
6568mod cluster_bc_reference_space_tests {
6569 use super::*;
6570
6571 fn acceleration_at_1100_fps(inputs: BallisticInputs) -> Vector3<f64> {
6572 let solver = TrajectorySolver::new(
6573 inputs,
6574 WindConditions::default(),
6575 AtmosphericConditions::default(),
6576 );
6577 let position = Vector3::zeros();
6578 let velocity = Vector3::new(1100.0 / 3.28084, 0.0, 0.0);
6579 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6580 solver.calculate_acceleration(
6581 &position,
6582 &velocity,
6583 &Vector3::zeros(),
6584 (temp_c, pressure_hpa, density / 1.225),
6585 )
6586 }
6587
6588 #[test]
6589 fn solver_passes_g7_reference_model_to_cluster_classifier() {
6590 let inputs = BallisticInputs {
6591 bc_value: 0.190,
6592 bc_type: DragModel::G7,
6593 bullet_mass: 77.0 * crate::constants::GRAINS_TO_KG,
6594 bullet_diameter: 0.224 * 0.0254,
6595 use_cluster_bc: true,
6596 ..BallisticInputs::default()
6597 };
6598
6599 let solver = TrajectorySolver::new(
6600 inputs,
6601 WindConditions::default(),
6602 AtmosphericConditions::default(),
6603 );
6604 let corrected = solver.apply_cluster_bc_correction(0.190, 2800.0);
6605
6606 assert!(
6607 (corrected / 0.190 - 1.004).abs() < 1e-12,
6608 "solver selected the wrong G7 cluster multiplier: {}",
6609 corrected / 0.190
6610 );
6611 }
6612
6613 #[test]
6614 fn velocity_bc_segments_are_not_cluster_corrected_twice() {
6615 let segmented_clustered = BallisticInputs {
6616 bc_value: 0.5,
6617 bc_type: DragModel::G7,
6618 use_bc_segments: true,
6619 bc_segments_data: Some(vec![
6620 crate::BCSegmentData {
6621 velocity_min: 0.0,
6622 velocity_max: 1_600.0,
6623 bc_value: 0.4,
6624 },
6625 crate::BCSegmentData {
6626 velocity_min: 1_600.0,
6627 velocity_max: 5_000.0,
6628 bc_value: 0.45,
6629 },
6630 ]),
6631 use_cluster_bc: true,
6632 ..BallisticInputs::default()
6633 };
6634 let mut segmented_only = segmented_clustered.clone();
6635 segmented_only.use_cluster_bc = false;
6636 let mut constant_clustered = segmented_clustered.clone();
6637 constant_clustered.bc_value = 0.4;
6638 constant_clustered.bc_segments_data = None;
6639
6640 let stacked = acceleration_at_1100_fps(segmented_clustered);
6641 let segment_only = acceleration_at_1100_fps(segmented_only);
6642 let cluster_only = acceleration_at_1100_fps(constant_clustered);
6643
6644 assert!(
6645 (stacked.x - segment_only.x).abs() < 1e-12,
6646 "segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
6647 stacked.x,
6648 segment_only.x
6649 );
6650 assert!(
6651 (cluster_only.x - segment_only.x).abs() > 1e-6,
6652 "cluster correction must remain active for a constant BC"
6653 );
6654 }
6655
6656 #[test]
6657 fn mach_bc_segments_are_not_cluster_corrected_twice() {
6658 let mach_segmented_clustered = BallisticInputs {
6659 bc_value: 0.5,
6660 bc_type: DragModel::G7,
6661 use_bc_segments: false,
6662 bc_segments: Some(vec![(0.5, 0.3), (1.5, 0.5)]),
6663 use_cluster_bc: true,
6664 ..BallisticInputs::default()
6665 };
6666 let mut mach_segmented_only = mach_segmented_clustered.clone();
6667 mach_segmented_only.use_cluster_bc = false;
6668
6669 let stacked = acceleration_at_1100_fps(mach_segmented_clustered);
6670 let segment_only = acceleration_at_1100_fps(mach_segmented_only);
6671
6672 assert!(
6673 (stacked.x - segment_only.x).abs() < 1e-12,
6674 "Mach segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
6675 stacked.x,
6676 segment_only.x
6677 );
6678 }
6679}
6680
6681#[cfg(test)]
6682mod velocity_bc_flag_tests {
6683 use super::*;
6684
6685 fn acceleration_at_600_mps(inputs: BallisticInputs) -> Vector3<f64> {
6686 let solver = TrajectorySolver::new(
6687 inputs,
6688 WindConditions::default(),
6689 AtmosphericConditions::default(),
6690 );
6691 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6692 solver.calculate_acceleration(
6693 &Vector3::zeros(),
6694 &Vector3::new(600.0, 0.0, 0.0),
6695 &Vector3::zeros(),
6696 (temp_c, pressure_hpa, density / 1.225),
6697 )
6698 }
6699
6700 #[test]
6701 fn velocity_bc_data_requires_opt_in_in_trajectory_solver() {
6702 let scalar_inputs = BallisticInputs {
6703 bc_value: 0.5,
6704 bc_type: DragModel::G7,
6705 ..BallisticInputs::default()
6706 };
6707 let mut disabled_inputs = scalar_inputs.clone();
6708 disabled_inputs.bc_segments_data = Some(vec![crate::BCSegmentData {
6709 velocity_min: 0.0,
6710 velocity_max: 4_000.0,
6711 bc_value: 0.46,
6712 }]);
6713 disabled_inputs.use_bc_segments = false;
6714 let mut enabled_inputs = disabled_inputs.clone();
6715 enabled_inputs.use_bc_segments = true;
6716 let mut mach_only_inputs = scalar_inputs.clone();
6717 mach_only_inputs.bc_segments = Some(vec![(0.0, 0.4), (3.0, 0.4)]);
6718 let mut disabled_with_both = mach_only_inputs.clone();
6719 disabled_with_both.bc_segments_data = disabled_inputs.bc_segments_data.clone();
6720
6721 let scalar = acceleration_at_600_mps(scalar_inputs);
6722 let disabled = acceleration_at_600_mps(disabled_inputs);
6723 let enabled = acceleration_at_600_mps(enabled_inputs);
6724 let mach_only = acceleration_at_600_mps(mach_only_inputs);
6725 let disabled_with_both = acceleration_at_600_mps(disabled_with_both);
6726
6727 assert_eq!(
6728 disabled.x.to_bits(),
6729 scalar.x.to_bits(),
6730 "a populated velocity table must not change drag while use_bc_segments is false"
6731 );
6732 assert!(
6733 enabled.x < disabled.x - 1.0,
6734 "enabling the lower BC table must increase drag: disabled ax={} enabled ax={}",
6735 disabled.x,
6736 enabled.x
6737 );
6738 assert_eq!(
6739 disabled_with_both.x.to_bits(),
6740 mach_only.x.to_bits(),
6741 "disabling velocity data must fall through to an explicit Mach table"
6742 );
6743 }
6744}
6745
6746#[cfg(test)]
6747mod mach_bc_segment_tests {
6748 use super::*;
6749
6750 #[test]
6751 fn trajectory_solver_interpolates_explicit_mach_bc_segments() {
6752 let segmented_inputs = BallisticInputs {
6753 bc_value: 0.8,
6754 use_bc_segments: false,
6755 bc_segments: Some(vec![(1.0, 0.2), (2.0, 0.4)]),
6756 bc_segments_data: None,
6757 ..BallisticInputs::default()
6758 };
6759
6760 let mut expected_inputs = segmented_inputs.clone();
6761 expected_inputs.bc_value = 0.3;
6762 expected_inputs.bc_segments = None;
6763
6764 let atmosphere = AtmosphericConditions::default();
6765 let segmented_solver = TrajectorySolver::new(
6766 segmented_inputs,
6767 WindConditions::default(),
6768 atmosphere.clone(),
6769 );
6770 let expected_solver = TrajectorySolver::new(
6771 expected_inputs,
6772 WindConditions::default(),
6773 atmosphere,
6774 );
6775 let position = Vector3::zeros();
6776 let (density, _, temp_c, pressure_hpa) = segmented_solver.resolved_atmosphere();
6777 let (_, local_speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
6778 segmented_solver.atmosphere.altitude,
6779 segmented_solver.atmosphere.altitude,
6780 temp_c,
6781 pressure_hpa,
6782 density / 1.225,
6783 segmented_solver.atmosphere.humidity,
6784 );
6785 let velocity = Vector3::new(1.5 * local_speed_of_sound, 0.0, 0.0);
6786 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
6787
6788 let segmented_acceleration = segmented_solver.calculate_acceleration(
6789 &position,
6790 &velocity,
6791 &Vector3::zeros(),
6792 resolved_atmo,
6793 );
6794 let expected_acceleration = expected_solver.calculate_acceleration(
6795 &position,
6796 &velocity,
6797 &Vector3::zeros(),
6798 resolved_atmo,
6799 );
6800
6801 assert!(
6802 (segmented_acceleration.x - expected_acceleration.x).abs() < 1e-12,
6803 "Mach 1.5 must interpolate BC 0.3: segmented ax={} expected ax={}",
6804 segmented_acceleration.x,
6805 expected_acceleration.x
6806 );
6807 }
6808}
6809
6810#[cfg(test)]
6811mod custom_drag_table_validation_tests {
6812 use super::*;
6813
6814 #[test]
6815 fn solve_accepts_zero_bc_when_custom_table_present() {
6816 let inputs = BallisticInputs {
6817 bc_value: 0.0, bullet_mass: 0.0106,
6819 bullet_diameter: 0.00782,
6820 muzzle_velocity: 850.0,
6821 custom_drag_table: Some(crate::drag::DragTable::new(
6822 vec![0.5, 1.0, 2.0, 3.0],
6823 vec![0.23, 0.40, 0.30, 0.26],
6824 )),
6825 ..BallisticInputs::default()
6826 };
6827 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
6828 assert!(solver.solve().is_ok());
6830 }
6831
6832 #[test]
6833 fn solve_still_requires_bc_without_table() {
6834 let inputs = BallisticInputs {
6835 bc_value: 0.0,
6836 bullet_mass: 0.0106,
6837 bullet_diameter: 0.00782,
6838 muzzle_velocity: 850.0,
6839 ..BallisticInputs::default()
6840 };
6841 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
6842 assert!(solver.solve().is_err());
6843 }
6844}
6845
6846#[cfg(test)]
6848mod cd_scale_tests {
6849 use super::*;
6850
6851 fn deck() -> crate::drag::DragTable {
6852 crate::drag::DragTable::new(vec![0.5, 1.0, 2.0, 3.0], vec![0.23, 0.40, 0.30, 0.26])
6853 }
6854
6855 fn deck_inputs(cd_scale: f64) -> BallisticInputs {
6856 BallisticInputs {
6857 bullet_mass: 0.0106,
6858 bullet_diameter: 0.00782,
6859 muzzle_velocity: 850.0,
6860 custom_drag_table: Some(deck()),
6861 cd_scale,
6862 ..BallisticInputs::default()
6863 }
6864 }
6865
6866 #[test]
6867 fn default_cd_scale_is_one() {
6868 assert_eq!(BallisticInputs::default().cd_scale, 1.0);
6869 }
6870
6871 #[test]
6874 fn cd_scale_absent_is_byte_identical_to_explicit_one() {
6875 let omitted = BallisticInputs {
6876 bullet_mass: 0.0106,
6877 bullet_diameter: 0.00782,
6878 muzzle_velocity: 850.0,
6879 custom_drag_table: Some(deck()),
6880 ..BallisticInputs::default()
6881 };
6882 let explicit = BallisticInputs {
6883 cd_scale: 1.0,
6884 ..omitted.clone()
6885 };
6886
6887 let solver_omitted =
6888 TrajectorySolver::new(omitted, WindConditions::default(), AtmosphericConditions::default());
6889 let solver_explicit =
6890 TrajectorySolver::new(explicit, WindConditions::default(), AtmosphericConditions::default());
6891
6892 let cd_omitted = solver_omitted.calculate_drag_coefficient(700.0, 340.0);
6893 let cd_explicit = solver_explicit.calculate_drag_coefficient(700.0, 340.0);
6894 assert_eq!(
6895 cd_omitted.to_bits(),
6896 cd_explicit.to_bits(),
6897 "default cd_scale must be bit-identical to an explicit 1.0"
6898 );
6899
6900 let result = solver_omitted.solve();
6903 assert!(result.is_ok(), "existing custom-deck solves must pass unchanged");
6904 }
6905
6906 #[test]
6908 fn cd_scale_multiplies_the_interpolated_cd_exactly() {
6909 let velocity = 700.0;
6910 let speed_of_sound = 340.0;
6911 let mach = velocity / speed_of_sound;
6912 let expected_unscaled = deck().interpolate(mach);
6913
6914 for &scale in &[0.90, 1.0, 1.10, 1.5] {
6915 let solver = TrajectorySolver::new(
6916 deck_inputs(scale),
6917 WindConditions::default(),
6918 AtmosphericConditions::default(),
6919 );
6920 let cd = solver.calculate_drag_coefficient(velocity, speed_of_sound);
6921 assert!(
6922 (cd - expected_unscaled * scale).abs() < 1e-12,
6923 "scale={scale}: cd={cd} expected={}",
6924 expected_unscaled * scale
6925 );
6926 }
6927 }
6928
6929 #[test]
6933 fn cd_scale_direction_on_cli_api_solver() {
6934 let solve = |scale: f64| {
6935 TrajectorySolver::new(
6936 deck_inputs(scale),
6937 WindConditions::default(),
6938 AtmosphericConditions::default(),
6939 )
6940 .solve()
6941 .expect("custom-deck solve should succeed")
6942 };
6943
6944 let baseline = solve(1.0);
6945 let scaled_up = solve(1.10);
6946 let scaled_down = solve(0.90);
6947
6948 assert!(
6949 scaled_up.impact_velocity < baseline.impact_velocity,
6950 "cd_scale=1.10 must increase drag -> lower impact velocity: base={} up={}",
6951 baseline.impact_velocity,
6952 scaled_up.impact_velocity
6953 );
6954 assert!(
6955 scaled_down.impact_velocity > baseline.impact_velocity,
6956 "cd_scale=0.90 must decrease drag -> higher impact velocity: base={} down={}",
6957 baseline.impact_velocity,
6958 scaled_down.impact_velocity
6959 );
6960 }
6961
6962 #[test]
6964 fn validate_for_solve_rejects_invalid_cd_scale() {
6965 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
6966 let solver = TrajectorySolver::new(
6967 deck_inputs(bad),
6968 WindConditions::default(),
6969 AtmosphericConditions::default(),
6970 );
6971 assert!(
6972 solver.solve().is_err(),
6973 "cd_scale={bad} must be rejected by validate_for_solve"
6974 );
6975 }
6976 }
6977
6978 #[test]
6985 fn validate_for_solve_rejects_invalid_cd_scale_without_a_custom_drag_table() {
6986 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
6987 let inputs = BallisticInputs {
6988 bc_value: 0.5,
6989 bc_type: crate::DragModel::G1,
6990 bullet_mass: 0.0106,
6991 bullet_diameter: 0.00782,
6992 muzzle_velocity: 850.0,
6993 cd_scale: bad,
6994 ..BallisticInputs::default()
6995 };
6996 assert!(inputs.custom_drag_table.is_none(), "precondition: no custom deck");
6997 let solver = TrajectorySolver::new(
6998 inputs,
6999 WindConditions::default(),
7000 AtmosphericConditions::default(),
7001 );
7002 assert!(
7003 solver.solve().is_err(),
7004 "cd_scale={bad} must be rejected by validate_for_solve even without a custom \
7005 drag table"
7006 );
7007 }
7008 }
7009
7010 #[test]
7014 fn cd_scale_is_inert_without_a_custom_drag_table() {
7015 let make = |cd_scale: f64| BallisticInputs {
7016 bc_value: 0.5,
7017 bc_type: crate::DragModel::G1,
7018 bullet_mass: 0.0106,
7019 bullet_diameter: 0.00782,
7020 muzzle_velocity: 850.0,
7021 cd_scale,
7022 ..BallisticInputs::default()
7023 };
7024 let solver_neutral = TrajectorySolver::new(
7025 make(1.0),
7026 WindConditions::default(),
7027 AtmosphericConditions::default(),
7028 );
7029 let solver_far = TrajectorySolver::new(
7030 make(1.5),
7031 WindConditions::default(),
7032 AtmosphericConditions::default(),
7033 );
7034 let cd_neutral = solver_neutral.calculate_drag_coefficient(700.0, 340.0);
7035 let cd_far = solver_far.calculate_drag_coefficient(700.0, 340.0);
7036 assert_eq!(
7037 cd_neutral.to_bits(),
7038 cd_far.to_bits(),
7039 "cd_scale must not affect the G-model/BC drag path"
7040 );
7041 }
7042
7043 #[test]
7046 fn cd_scale_shifts_all_three_solver_paths_in_the_same_direction() {
7047 let cli_solve = |scale: f64| {
7049 TrajectorySolver::new(
7050 deck_inputs(scale),
7051 WindConditions::default(),
7052 AtmosphericConditions::default(),
7053 )
7054 .solve()
7055 .expect("cli_api custom-deck solve should succeed")
7056 };
7057 let cli_baseline = cli_solve(1.0);
7058 let cli_scaled = cli_solve(1.10);
7059 assert!(
7060 cli_scaled.impact_velocity < cli_baseline.impact_velocity,
7061 "cli_api: cd_scale=1.10 must lower impact velocity"
7062 );
7063
7064 let derivatives_accel_x = |scale: f64| {
7066 let inputs = deck_inputs(scale);
7067 crate::derivatives::compute_derivatives(
7068 nalgebra::Vector3::zeros(),
7069 nalgebra::Vector3::new(700.0, 0.0, 0.0),
7070 &inputs,
7071 nalgebra::Vector3::zeros(),
7072 (1.225, 340.0, 0.0, 0.0),
7073 inputs.bc_value,
7074 None,
7075 0.0,
7076 None,
7077 )[3]
7078 };
7079 let deriv_baseline = derivatives_accel_x(1.0);
7080 let deriv_scaled = derivatives_accel_x(1.10);
7081 assert!(
7082 deriv_scaled < deriv_baseline,
7083 "derivatives: cd_scale=1.10 must make x-acceleration more negative (more drag): \
7084 base={deriv_baseline} scaled={deriv_scaled}"
7085 );
7086
7087 let fast_final_speed = |scale: f64| {
7089 let inputs = deck_inputs(scale);
7090 let wind_sock = crate::wind::WindSock::new(vec![]);
7091 let params = crate::fast_trajectory::FastIntegrationParams {
7092 horiz: 500.0,
7093 vert: 0.0,
7094 initial_state: [0.0, 0.0, 0.0, 850.0, 0.0, 0.0],
7095 t_span: (0.0, 5.0),
7096 atmo_params: (0.0, 15.0, 1013.25, 1.0),
7097 atmo_sock: None,
7098 };
7099 let solution = crate::fast_trajectory::fast_integrate(&inputs, &wind_sock, params);
7100 assert!(solution.success, "fast_integrate must succeed for scale={scale}");
7101 let last = solution.t.len() - 1;
7102 let (vx, vy, vz) = (
7103 solution.y[3][last],
7104 solution.y[4][last],
7105 solution.y[5][last],
7106 );
7107 (vx * vx + vy * vy + vz * vz).sqrt()
7108 };
7109 let fast_baseline = fast_final_speed(1.0);
7110 let fast_scaled = fast_final_speed(1.10);
7111 assert!(
7112 fast_scaled < fast_baseline,
7113 "fast_trajectory: cd_scale=1.10 must lower final speed: base={fast_baseline} scaled={fast_scaled}"
7114 );
7115 }
7116}
7117
7118#[cfg(test)]
7119mod humid_local_mach_tests {
7120 use super::*;
7121
7122 fn solver_with_station_humidity(humidity_percent: f64) -> TrajectorySolver {
7123 let inputs = BallisticInputs {
7124 custom_drag_table: Some(crate::drag::DragTable::new(vec![0.5, 1.5], vec![0.1, 1.1])),
7125 ..BallisticInputs::default()
7126 };
7127 TrajectorySolver::new(
7128 inputs,
7129 WindConditions::default(),
7130 AtmosphericConditions {
7131 temperature: 30.0,
7132 pressure: 1013.25,
7133 humidity: humidity_percent,
7134 altitude: 0.0,
7135 },
7136 )
7137 }
7138
7139 fn acceleration(solver: &TrajectorySolver, base_ratio: f64) -> Vector3<f64> {
7140 solver.calculate_acceleration(
7141 &Vector3::zeros(),
7142 &Vector3::new(350.0, 0.0, 0.0),
7143 &Vector3::zeros(),
7144 (30.0, 1013.25, base_ratio),
7145 )
7146 }
7147
7148 #[test]
7149 fn local_mach_uses_station_humidity_when_density_is_held_constant() {
7150 let dry = acceleration(&solver_with_station_humidity(0.0), 1.0);
7151 let humid = acceleration(&solver_with_station_humidity(100.0), 1.0);
7152
7153 assert!(
7154 humid.x > dry.x,
7155 "humid sound speed should lower Mach and drag on the rising test curve: dry ax={} humid ax={}",
7156 dry.x,
7157 humid.x
7158 );
7159 }
7160
7161 #[test]
7162 fn active_atmosphere_zone_uses_zone_humidity_instead_of_station_humidity() {
7163 let zone_humidity = 80.0;
7164 let zone_ratio =
7165 crate::atmosphere::calculate_air_density_cimp(30.0, 1013.25, zone_humidity) / 1.225;
7166 let station_solver = solver_with_station_humidity(zone_humidity);
7167 let mut zoned_solver = solver_with_station_humidity(0.0);
7168 zoned_solver.set_atmo_segments(vec![(30.0, 1013.25, zone_humidity, 1_000.0)]);
7169
7170 let station = acceleration(&station_solver, zone_ratio);
7171 let zoned = acceleration(&zoned_solver, zone_ratio);
7172
7173 assert!(
7174 (zoned - station).norm() < 1e-12,
7175 "active zone T/P/RH should override the station atmosphere: station={station:?} zoned={zoned:?}"
7176 );
7177 }
7178}
7179
7180#[cfg(test)]
7181mod inclined_atmosphere_frame_tests {
7182 use super::*;
7183
7184 fn expected_shot_frame_vector(level: Vector3<f64>, angle: f64) -> Vector3<f64> {
7185 let (sin_angle, cos_angle) = angle.sin_cos();
7186 Vector3::new(
7187 level.x * cos_angle + level.y * sin_angle,
7188 -level.x * sin_angle + level.y * cos_angle,
7189 level.z,
7190 )
7191 }
7192
7193 #[test]
7194 fn inclined_positions_at_same_world_altitude_have_same_solver_acceleration() {
7195 let angle = std::f64::consts::FRAC_PI_6;
7196 let inputs = BallisticInputs {
7197 shooting_angle: angle,
7198 ..BallisticInputs::default()
7199 };
7200 let atmosphere = AtmosphericConditions {
7201 altitude: 100.0,
7202 ..AtmosphericConditions::default()
7203 };
7204 let solver = TrajectorySolver::new(inputs, WindConditions::default(), atmosphere);
7205 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7206 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
7207 let velocity = Vector3::new(600.0, 0.0, 0.0);
7208 let along_slant = Vector3::new(1_000.0, 0.0, 0.0);
7209 let across_slant = Vector3::new(0.0, 500.0 / angle.cos(), 0.0);
7210
7211 let a = solver.calculate_acceleration(
7212 &along_slant,
7213 &velocity,
7214 &Vector3::zeros(),
7215 resolved_atmo,
7216 );
7217 let b = solver.calculate_acceleration(
7218 &across_slant,
7219 &velocity,
7220 &Vector3::zeros(),
7221 resolved_atmo,
7222 );
7223
7224 assert!(
7225 (a - b).norm() < 1e-10,
7226 "solver acceleration differs at equal world altitude: {a:?} vs {b:?}"
7227 );
7228 }
7229
7230 #[test]
7231 fn inclined_headwind_is_rotated_into_solver_frame() {
7232 let angle = std::f64::consts::FRAC_PI_6;
7233 let inputs = BallisticInputs {
7234 shooting_angle: angle,
7235 ..BallisticInputs::default()
7236 };
7237 let solver = TrajectorySolver::new(
7238 inputs,
7239 WindConditions::default(),
7240 AtmosphericConditions::default(),
7241 );
7242 let level_headwind = Vector3::new(-100.0, 0.0, 0.0);
7243 let velocity = expected_shot_frame_vector(level_headwind, angle);
7244 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7245 let actual = solver.calculate_acceleration(
7246 &Vector3::zeros(),
7247 &velocity,
7248 &level_headwind,
7249 (temp_c, pressure_hpa, density / 1.225),
7250 );
7251
7252 assert!(
7253 (actual - solver.gravity_acceleration()).norm() < 1e-12,
7254 "co-moving horizontal wind must leave only shot-frame gravity: {actual:?}"
7255 );
7256 }
7257
7258 #[test]
7259 fn inclined_coriolis_is_rotated_into_solver_frame() {
7260 let angle = std::f64::consts::FRAC_PI_6;
7261 let latitude_deg = 45.0_f64;
7262 let shot_azimuth = 0.4_f64;
7263 let velocity = Vector3::new(600.0, 20.0, 5.0);
7264 let base_inputs = BallisticInputs {
7265 shooting_angle: angle,
7266 latitude: Some(latitude_deg),
7267 shot_azimuth,
7268 ..BallisticInputs::default()
7269 };
7270 let acceleration = |enable_coriolis| {
7271 let mut inputs = base_inputs.clone();
7272 inputs.enable_coriolis = enable_coriolis;
7273 let solver = TrajectorySolver::new(
7274 inputs,
7275 WindConditions::default(),
7276 AtmosphericConditions::default(),
7277 );
7278 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7279 solver.calculate_acceleration(
7280 &Vector3::zeros(),
7281 &velocity,
7282 &Vector3::zeros(),
7283 (temp_c, pressure_hpa, density / 1.225),
7284 )
7285 };
7286
7287 let omega_earth = 7.2921159e-5_f64;
7288 let latitude = latitude_deg.to_radians();
7289 let level_omega = Vector3::new(
7290 omega_earth * latitude.cos() * shot_azimuth.cos(),
7291 omega_earth * latitude.sin(),
7292 -omega_earth * latitude.cos() * shot_azimuth.sin(),
7293 );
7294 let expected = -2.0 * expected_shot_frame_vector(level_omega, angle).cross(&velocity);
7295 let actual = acceleration(true) - acceleration(false);
7296
7297 assert!(
7298 (actual - expected).norm() < 1e-12,
7299 "inclined Coriolis mismatch: actual={actual:?}, expected={expected:?}"
7300 );
7301 }
7302}
7303
7304#[cfg(test)]
7305mod terminal_range_interpolation_tests {
7306 use super::*;
7307
7308 #[test]
7309 fn terminal_finalizer_selects_the_earliest_crossed_boundary() {
7310 let inputs = BallisticInputs {
7311 ground_threshold: 0.0,
7312 ..BallisticInputs::default()
7313 };
7314 let mut solver = TrajectorySolver::new(
7315 inputs,
7316 WindConditions::default(),
7317 AtmosphericConditions::default(),
7318 );
7319 solver.set_max_range(120.0);
7320
7321 let previous_speed = 700.0;
7322 let mut points = vec![TrajectoryPoint {
7323 time: 99.0,
7324 position: Vector3::new(90.0, 1.0, -1.0),
7325 velocity_magnitude: previous_speed,
7326 kinetic_energy: 0.5 * solver.inputs.bullet_mass * previous_speed.powi(2),
7327 drag_coefficient: None,
7328 }];
7329 let mut max_height = 1.0;
7330 let termination = solver
7331 .append_terminal_endpoint(
7332 &mut points,
7333 Vector3::new(130.0, -3.0, 3.0),
7334 Vector3::new(600.0, 0.0, 0.0),
7335 101.0,
7336 &mut max_height,
7337 )
7338 .expect("the final step brackets supported boundaries");
7339
7340 assert_eq!(termination, TrajectoryTermination::GroundThreshold);
7341 assert_eq!(points.len(), 2);
7342 let terminal = points.last().expect("terminal point");
7343 assert_eq!(terminal.time, 99.5);
7344 assert_eq!(terminal.position, Vector3::new(100.0, 0.0, 0.0));
7345 assert_eq!(terminal.velocity_magnitude, 675.0);
7346 assert_eq!(
7347 terminal.kinetic_energy,
7348 0.5 * solver.inputs.bullet_mass * 675.0_f64.powi(2)
7349 );
7350
7351 solver.set_max_range(100.0);
7353 let mut tied_points = vec![points[0].clone()];
7354 assert_eq!(
7355 solver
7356 .append_terminal_endpoint(
7357 &mut tied_points,
7358 Vector3::new(130.0, -3.0, 3.0),
7359 Vector3::new(600.0, 0.0, 0.0),
7360 101.0,
7361 &mut max_height,
7362 )
7363 .expect("tied boundaries remain a valid terminal"),
7364 TrajectoryTermination::GroundThreshold
7365 );
7366 }
7367
7368 #[test]
7369 fn sub_ulp_terminal_crossing_replaces_instead_of_duplicating_range() {
7370 let ground_threshold = f64::from_bits(1.0_f64.to_bits() - 1);
7371 let inputs = BallisticInputs {
7372 ground_threshold,
7373 ..BallisticInputs::default()
7374 };
7375 let mut solver = TrajectorySolver::new(
7376 inputs,
7377 WindConditions::default(),
7378 AtmosphericConditions::default(),
7379 );
7380 solver.set_max_range(1_000.0);
7381
7382 let speed = 700.0;
7383 let mut points = vec![TrajectoryPoint {
7384 time: 0.0,
7385 position: Vector3::new(100.0, 1.0, 0.0),
7386 velocity_magnitude: speed,
7387 kinetic_energy: 0.5 * solver.inputs.bullet_mass * speed.powi(2),
7388 drag_coefficient: None,
7389 }];
7390 let mut max_height = 1.0;
7391 let termination = solver
7392 .append_terminal_endpoint(
7393 &mut points,
7394 Vector3::new(101.0, 0.0, 0.0),
7395 Vector3::new(699.0, 0.0, 0.0),
7396 1.0,
7397 &mut max_height,
7398 )
7399 .expect("sub-ULP ground crossing remains representable as one terminal state");
7400
7401 assert_eq!(termination, TrajectoryTermination::GroundThreshold);
7402 assert_eq!(points.len(), 1);
7403 assert_eq!(points[0].position.x, 100.0);
7404 assert_eq!(points[0].position.y.to_bits(), ground_threshold.to_bits());
7405 assert!(points[0].time > 0.0);
7406 }
7407
7408 #[test]
7409 fn every_solver_appends_an_exact_max_range_endpoint() {
7410 let target_range = 0.1;
7411 let modes = [
7412 ("Euler", false, false),
7413 ("RK4", true, false),
7414 ("RK45", true, true),
7415 ];
7416
7417 for (name, use_rk4, use_adaptive_rk45) in modes {
7418 let inputs = BallisticInputs {
7419 use_rk4,
7420 use_adaptive_rk45,
7421 ground_threshold: f64::NEG_INFINITY,
7422 enable_trajectory_sampling: true,
7423 sample_interval: target_range,
7424 ..BallisticInputs::default()
7425 };
7426 let mut solver = TrajectorySolver::new(
7427 inputs,
7428 WindConditions::default(),
7429 AtmosphericConditions::default(),
7430 );
7431 solver.set_max_range(target_range);
7432
7433 let result = solver.solve().expect("short-range solve should succeed");
7434 let terminal = result.points.last().expect("terminal point is missing");
7435 let muzzle = result.points.first().expect("muzzle point is missing");
7436
7437 assert_eq!(result.termination, TrajectoryTermination::MaxRange);
7438 assert_eq!(
7439 terminal.position.x.to_bits(),
7440 target_range.to_bits(),
7441 "{name} did not terminate exactly at max_range"
7442 );
7443 assert_eq!(result.max_range.to_bits(), target_range.to_bits());
7444 assert!(
7445 result.time_of_flight > 0.0 && result.time_of_flight < solver.time_step,
7446 "{name} terminal time was not interpolated within the crossing step: {}",
7447 result.time_of_flight
7448 );
7449 assert_eq!(result.time_of_flight.to_bits(), terminal.time.to_bits());
7450 assert_eq!(
7451 result.impact_velocity.to_bits(),
7452 terminal.velocity_magnitude.to_bits()
7453 );
7454 assert_eq!(
7455 result.impact_energy.to_bits(),
7456 terminal.kinetic_energy.to_bits()
7457 );
7458 let expected_energy = 0.5 * solver.inputs.bullet_mass * result.impact_velocity.powi(2);
7459 assert!((result.impact_energy - expected_energy).abs() < 1e-12);
7460 assert!(terminal.velocity_magnitude < muzzle.velocity_magnitude);
7461 assert!(terminal.kinetic_energy < muzzle.kinetic_energy);
7462
7463 let terminal_sample = result
7464 .sampled_points
7465 .as_ref()
7466 .and_then(|samples| samples.last())
7467 .expect("terminal trajectory sample is missing");
7468 assert_eq!(
7469 terminal_sample.distance_m.to_bits(),
7470 target_range.to_bits(),
7471 "{name} sampling did not include max_range"
7472 );
7473 assert_eq!(
7474 terminal_sample.time_s.to_bits(),
7475 result.time_of_flight.to_bits()
7476 );
7477 assert_eq!(
7478 terminal_sample.velocity_mps.to_bits(),
7479 result.impact_velocity.to_bits()
7480 );
7481 assert!((terminal_sample.energy_j - result.impact_energy).abs() < 1e-12);
7482 }
7483 }
7484}
7485
7486#[cfg(test)]
7487mod precession_inertia_wiring_tests {
7488 use super::*;
7489
7490 #[test]
7491 fn solver_uses_projectile_specific_moments_of_inertia() {
7492 let mass_kg = 55.0 * crate::constants::GRAINS_TO_KG;
7493 let caliber_m = 0.224 * 0.0254;
7494 let length_m = 0.75 * 0.0254;
7495 let inputs = BallisticInputs {
7496 bullet_mass: mass_kg,
7497 bullet_diameter: caliber_m,
7498 bullet_length: length_m,
7499 muzzle_velocity: 800.0,
7500 twist_rate: 7.0,
7501 enable_precession_nutation: true,
7502 use_rk4: false,
7503 use_adaptive_rk45: false,
7504 ..BallisticInputs::default()
7505 };
7506 let mut solver = TrajectorySolver::new(
7507 inputs,
7508 WindConditions::default(),
7509 AtmosphericConditions::default(),
7510 );
7511 solver.set_max_range(0.1);
7512
7513 let (air_density, speed_of_sound, _, _) = solver.resolved_atmosphere();
7514 let velocity_mps = solver.inputs.muzzle_velocity;
7515 let velocity_fps = velocity_mps * 3.28084;
7516 let twist_rate_ft = solver.inputs.twist_rate / 12.0;
7517 let spin_rate_rad_s = (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI;
7518 let initial_state = AngularState {
7519 pitch_angle: 0.001,
7520 yaw_angle: 0.001,
7521 pitch_rate: 0.0,
7522 yaw_rate: 0.0,
7523 precession_angle: 0.0,
7524 nutation_phase: 0.0,
7525 };
7526 let params = PrecessionNutationParams {
7527 mass_kg,
7528 caliber_m,
7529 length_m,
7530 spin_rate_rad_s,
7531 spin_inertia: crate::spin_decay::calculate_moment_of_inertia(
7532 mass_kg, caliber_m, length_m, "ogive",
7533 ),
7534 transverse_inertia: crate::pitch_damping::calculate_transverse_moment_of_inertia(
7535 mass_kg, caliber_m, length_m, "ogive",
7536 ),
7537 velocity_mps,
7538 air_density_kg_m3: air_density,
7539 mach: velocity_mps / speed_of_sound,
7540 pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
7541 nutation_damping_factor: 0.05,
7542 };
7543 let expected = calculate_combined_angular_motion(
7544 ¶ms,
7545 &initial_state,
7546 0.0,
7547 solver.time_step,
7548 0.001,
7549 );
7550 let actual = solver
7551 .solve()
7552 .expect("one-step solve should succeed")
7553 .angular_state
7554 .expect("precession/nutation was enabled");
7555
7556 assert!(
7557 (actual.precession_angle - expected.precession_angle).abs() < 1e-15,
7558 "precession phase used the wrong inertia: actual={}, expected={}",
7559 actual.precession_angle,
7560 expected.precession_angle
7561 );
7562 assert!(
7563 (actual.nutation_phase - expected.nutation_phase).abs() < 1e-15,
7564 "nutation phase used the wrong inertia: actual={}, expected={}",
7565 actual.nutation_phase,
7566 expected.nutation_phase
7567 );
7568 }
7569}
7570
7571#[cfg(test)]
7572mod form_factor_drag_tests {
7573 use super::*;
7574
7575 fn acceleration_with_form_factor_flag(enabled: bool) -> Vector3<f64> {
7576 let inputs = BallisticInputs {
7577 bc_value: 0.462,
7578 bc_type: DragModel::G1,
7579 bullet_model: Some("168gr SMK Match".to_string()),
7580 use_form_factor: enabled,
7581 ..BallisticInputs::default()
7582 };
7583 let solver = TrajectorySolver::new(
7584 inputs,
7585 WindConditions::default(),
7586 AtmosphericConditions::default(),
7587 );
7588 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7589 solver.calculate_acceleration(
7590 &Vector3::zeros(),
7591 &Vector3::new(600.0, 0.0, 0.0),
7592 &Vector3::zeros(),
7593 (temp_c, pressure_hpa, density / 1.225),
7594 )
7595 }
7596
7597 #[test]
7598 fn measured_bc_drag_does_not_apply_name_based_form_factor_again() {
7599 let baseline = acceleration_with_form_factor_flag(false);
7600 let flagged = acceleration_with_form_factor_flag(true);
7601
7602 assert!(
7603 (flagged - baseline).norm() < 1e-12,
7604 "published BC already encodes form factor: baseline={baseline:?} flagged={flagged:?}"
7605 );
7606 }
7607}
7608
7609#[cfg(test)]
7610mod rk45_adaptivity_tests {
7611 use super::*;
7612
7613 #[test]
7614 fn cli_rk45_error_norm_scales_components_independently() {
7615 let position = Vector3::new(1.0e9, 0.0, 0.0);
7616 let velocity = Vector3::new(800.0, 0.0, 0.0);
7617 let fifth_position = position;
7618 let fifth_velocity = velocity;
7619 let fourth_position = position;
7620 let fourth_velocity = Vector3::new(800.0, 1.0e-3, 0.0);
7621
7622 let error = cli_rk45_error_norm(
7623 &position,
7624 &velocity,
7625 &fifth_position,
7626 &fifth_velocity,
7627 &fourth_position,
7628 &fourth_velocity,
7629 );
7630 let expected = 1.0e-3 / 6.0_f64.sqrt();
7631
7632 assert!(
7633 (error - expected).abs() <= 1e-15,
7634 "large downrange position masked a velocity-component error: {error}"
7635 );
7636 }
7637
7638 fn discontinuous_wind_solver() -> TrajectorySolver {
7639 let inputs = BallisticInputs::default();
7640 let mut solver = TrajectorySolver::new(
7641 inputs,
7642 WindConditions::default(),
7643 AtmosphericConditions::default(),
7644 );
7645 solver.set_wind_segments(vec![
7646 crate::wind::WindSegment::new(0.0, 90.0, 4.0),
7647 crate::wind::WindSegment::new(1_000.0, 90.0, 10_000.0),
7648 ]);
7649 solver
7650 }
7651
7652 #[test]
7653 fn rk45_retries_discontinuous_trial_before_advancing() {
7654 let solver = discontinuous_wind_solver();
7655 let position = Vector3::new(0.0, solver.inputs.muzzle_height, 0.0);
7656 let velocity = Vector3::new(solver.inputs.muzzle_velocity, 0.0, 0.0);
7657 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7658 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
7659 let dt = 0.01;
7660
7661 let rejected_trial = solver.rk45_step(
7662 &position,
7663 &velocity,
7664 dt,
7665 &Vector3::zeros(),
7666 RK45_TOLERANCE,
7667 resolved_atmo,
7668 );
7669 assert!(
7670 rejected_trial.error > RK45_TOLERANCE,
7671 "discontinuous full step must exceed tolerance, got {}",
7672 rejected_trial.error
7673 );
7674
7675 let accepted = solver.adaptive_rk45_step(
7676 &position,
7677 &velocity,
7678 dt,
7679 &Vector3::zeros(),
7680 resolved_atmo,
7681 );
7682 assert!(accepted.used_dt < dt, "oversized trial was not retried");
7683 assert!(
7684 accepted.error <= RK45_TOLERANCE || accepted.used_dt <= RK45_MIN_DT,
7685 "accepted error {} exceeds tolerance at dt {}",
7686 accepted.error,
7687 accepted.used_dt
7688 );
7689
7690 let accepted_trial = solver.rk45_step(
7691 &position,
7692 &velocity,
7693 accepted.used_dt,
7694 &Vector3::zeros(),
7695 RK45_TOLERANCE,
7696 resolved_atmo,
7697 );
7698 assert_eq!(accepted.position, accepted_trial.position);
7699 assert_eq!(accepted.velocity, accepted_trial.velocity);
7700 assert!((RK45_MIN_DT..=RK45_MAX_DT).contains(&accepted.next_dt));
7701 }
7702}
7703
7704#[cfg(test)]
7705mod ground_termination_tests {
7706 use super::*;
7707 use crate::trajectory_observation::TrajectoryObservationFlag;
7708
7709 #[test]
7710 fn every_solver_reports_one_exact_early_ground_endpoint() {
7711 for (name, use_rk4, use_adaptive_rk45) in [
7712 ("Euler", false, false),
7713 ("RK4", true, false),
7714 ("RK45", true, true),
7715 ] {
7716 let inputs = BallisticInputs {
7717 muzzle_height: 1.0,
7718 muzzle_angle: -0.2,
7719 ground_threshold: 0.0,
7720 use_rk4,
7721 use_adaptive_rk45,
7722 ..BallisticInputs::default()
7723 };
7724 let mut solver = TrajectorySolver::new(
7725 inputs,
7726 WindConditions::default(),
7727 AtmosphericConditions::default(),
7728 );
7729 solver.set_max_range(1_000.0);
7730
7731 let result = solver.solve().expect("early-ground solve should succeed");
7732 let terminal = result.points.last().expect("terminal point is missing");
7733
7734 assert_eq!(result.termination, TrajectoryTermination::GroundThreshold);
7735 assert_eq!(terminal.position.y.to_bits(), 0.0_f64.to_bits());
7736 assert!(
7737 terminal.position.x < 1_000.0,
7738 "{name} incorrectly reached max range"
7739 );
7740 assert_eq!(result.max_range.to_bits(), terminal.position.x.to_bits());
7741 assert_eq!(
7742 result
7743 .points
7744 .iter()
7745 .filter(|point| point.position.y == 0.0)
7746 .count(),
7747 1,
7748 "{name} did not retain exactly one ground endpoint"
7749 );
7750
7751 let observations = result
7752 .sample_observations(1.0, 100)
7753 .expect("checked early-ground sampling should succeed");
7754 assert!(observations[..observations.len() - 1]
7755 .iter()
7756 .all(|observation| observation.distance_m < terminal.position.x));
7757 let terminal_observation = observations.last().expect("terminal observation");
7758 assert_eq!(
7759 terminal_observation.distance_m.to_bits(),
7760 terminal.position.x.to_bits()
7761 );
7762 assert!(terminal_observation
7763 .flags
7764 .contains(&TrajectoryObservationFlag::Terminal));
7765 assert!(terminal_observation
7766 .flags
7767 .contains(&TrajectoryObservationFlag::GroundThreshold));
7768 assert_eq!(
7769 observations
7770 .iter()
7771 .filter(|observation| observation
7772 .flags
7773 .contains(&TrajectoryObservationFlag::Terminal))
7774 .count(),
7775 1,
7776 "{name} repeated the terminal observation"
7777 );
7778 }
7779 }
7780
7781 #[test]
7786 fn rk4_and_rk45_descend_to_ground_threshold() {
7787 for adaptive in [false, true] {
7788 let inputs = BallisticInputs {
7789 muzzle_angle: 0.1, use_rk4: true,
7791 use_adaptive_rk45: adaptive,
7792 ..BallisticInputs::default()
7793 };
7794 assert_eq!(
7795 inputs.ground_threshold, -100.0,
7796 "default ground_threshold is -100 m"
7797 );
7798
7799 let mut solver = TrajectorySolver::new(
7800 inputs,
7801 WindConditions::default(),
7802 AtmosphericConditions::default(),
7803 );
7804 solver.set_max_range(1.0e7);
7806
7807 let result = solver.solve().expect("solve should succeed");
7808 let final_y = result
7809 .points
7810 .last()
7811 .expect("trajectory has points")
7812 .position
7813 .y;
7814 assert!(
7815 final_y < -1.0,
7816 "adaptive_rk45={adaptive}: final y = {final_y} m; a lofted shot should descend \
7817 past launch level toward the ground_threshold floor, not stop at y = 0"
7818 );
7819 }
7820 }
7821}
7822
7823#[cfg(test)]
7824mod magnus_stability_tests {
7825 use super::*;
7826
7827 #[test]
7828 fn yaw_of_repose_magnus_force_is_vertical_and_twist_signed() {
7829 let acceleration = |enable_magnus, is_twist_right| {
7830 let inputs = BallisticInputs {
7831 muzzle_velocity: 822.96,
7832 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
7833 bullet_diameter: 0.308 * 0.0254,
7834 bullet_length: 1.215 * 0.0254,
7835 twist_rate: 10.0,
7836 is_twist_right,
7837 enable_magnus,
7838 ..BallisticInputs::default()
7839 };
7840 let solver = TrajectorySolver::new(
7841 inputs,
7842 WindConditions::default(),
7843 AtmosphericConditions::default(),
7844 );
7845 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7846 solver.calculate_acceleration(
7847 &Vector3::zeros(),
7848 &Vector3::new(822.96, 0.0, 0.0),
7849 &Vector3::zeros(),
7850 (temp_c, pressure_hpa, density / 1.225),
7851 )
7852 };
7853
7854 let baseline = acceleration(false, true);
7855 let right_twist = acceleration(true, true) - baseline;
7856 let left_twist = acceleration(true, false) - baseline;
7857
7858 assert!(
7859 right_twist.y < 0.0,
7860 "right-hand Magnus must point down, got {right_twist:?}"
7861 );
7862 assert!(
7863 left_twist.y > 0.0,
7864 "left-hand Magnus must point up, got {left_twist:?}"
7865 );
7866 assert!((right_twist.y + left_twist.y).abs() < 1e-12);
7867 assert!(right_twist.x.abs() < 1e-12 && right_twist.z.abs() < 1e-12);
7868 assert!(left_twist.x.abs() < 1e-12 && left_twist.z.abs() < 1e-12);
7869 }
7870
7871 #[test]
7872 fn magnus_uses_velocity_corrected_muzzle_stability_gate() {
7873 let muzzle_velocity = 1_400.0 / 3.28084;
7874 let inputs = BallisticInputs {
7875 muzzle_velocity,
7876 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
7877 bullet_diameter: 0.308 * 0.0254,
7878 bullet_length: 1.215 * 0.0254,
7879 twist_rate: 15.0,
7880 enable_magnus: true,
7881 ..BallisticInputs::default()
7882 };
7883 let solver = TrajectorySolver::new(
7884 inputs.clone(),
7885 WindConditions::default(),
7886 AtmosphericConditions::default(),
7887 );
7888
7889 let bare_sg = crate::spin_drift::miller_stability(0.308, 168.0, 15.0, 1.215);
7890 let canonical_sg = solver.effective_spin_drift_sg();
7891 assert!(bare_sg > 1.0, "test requires bare Sg above the Magnus gate");
7892 assert!(
7893 canonical_sg < 1.0,
7894 "velocity-corrected Sg must be below the gate, got {canonical_sg}"
7895 );
7896
7897 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7898 let acceleration = solver.calculate_acceleration(
7899 &Vector3::zeros(),
7900 &Vector3::new(muzzle_velocity, 0.0, 0.0),
7901 &Vector3::zeros(),
7902 (temp_c, pressure_hpa, density / 1.225),
7903 );
7904 let mut baseline_inputs = inputs;
7905 baseline_inputs.enable_magnus = false;
7906 let baseline_solver = TrajectorySolver::new(
7907 baseline_inputs,
7908 WindConditions::default(),
7909 AtmosphericConditions::default(),
7910 );
7911 let baseline = baseline_solver.calculate_acceleration(
7912 &Vector3::zeros(),
7913 &Vector3::new(muzzle_velocity, 0.0, 0.0),
7914 &Vector3::zeros(),
7915 (temp_c, pressure_hpa, density / 1.225),
7916 );
7917
7918 assert_eq!(
7919 acceleration, baseline,
7920 "canonical Sg below 1 must suppress every Magnus acceleration component"
7921 );
7922 }
7923
7924 #[test]
7925 fn magnus_force_grows_as_fixed_spin_projectile_slows() {
7926 let inputs = BallisticInputs {
7927 muzzle_velocity: 800.0,
7928 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
7929 bullet_diameter: 0.308 * 0.0254,
7930 bullet_length: 1.215 * 0.0254,
7931 twist_rate: 12.0,
7932 enable_magnus: true,
7933 ..BallisticInputs::default()
7934 };
7935
7936 let magnus_acceleration = |speed_mps| {
7937 let evaluate = |enable_magnus| {
7938 let mut run_inputs = inputs.clone();
7939 run_inputs.enable_magnus = enable_magnus;
7940 let solver = TrajectorySolver::new(
7941 run_inputs,
7942 WindConditions::default(),
7943 AtmosphericConditions::default(),
7944 );
7945 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7946 solver
7947 .calculate_acceleration(
7948 &Vector3::zeros(),
7949 &Vector3::new(speed_mps, 0.0, 0.0),
7950 &Vector3::zeros(),
7951 (temp_c, pressure_hpa, density / 1.225),
7952 )
7953 .y
7954 };
7955 (evaluate(true) - evaluate(false)).abs()
7956 };
7957
7958 let fast = magnus_acceleration(200.0);
7959 let slow = magnus_acceleration(100.0);
7960 let ratio = slow / fast;
7961 let expected_ratio = 2.0_f64.powf(5.0 / 3.0);
7962
7963 assert!(fast > 0.0 && slow > 0.0, "fast={fast}, slow={slow}");
7964 assert!(
7965 (ratio - expected_ratio).abs() < 1e-3,
7966 "fixed-spin Magnus acceleration must grow downrange; slow/fast={ratio}, \
7967 expected={expected_ratio}"
7968 );
7969 }
7970}
7971
7972#[cfg(test)]
7973mod coriolis_direction_tests {
7974 use super::*;
7975 use std::f64::consts::FRAC_PI_2;
7976
7977 #[test]
7978 fn supersonic_crossing_flags_a_positive_range_sample() {
7979 use crate::trajectory_sampling::TrajectoryFlag;
7983
7984 for (solver_name, use_rk4, use_adaptive_rk45) in [
7985 ("Euler", false, false),
7986 ("RK4", true, false),
7987 ("RK45", true, true),
7988 ] {
7989 let inputs = BallisticInputs {
7990 muzzle_velocity: 850.0,
7991 bc_value: 0.2,
7992 bc_type: DragModel::G7,
7993 muzzle_angle: 0.03,
7994 enable_trajectory_sampling: true,
7995 sample_interval: 50.0,
7996 use_rk4,
7997 use_adaptive_rk45,
7998 ..BallisticInputs::default()
7999 };
8000 let mut solver = TrajectorySolver::new(
8001 inputs,
8002 WindConditions::default(),
8003 AtmosphericConditions::default(),
8004 );
8005 solver.set_max_range(2000.0);
8006 let samples = solver
8007 .solve()
8008 .expect("supersonic solve should succeed")
8009 .sampled_points
8010 .expect("sampling was enabled");
8011 let flagged_distances: Vec<_> = samples
8012 .iter()
8013 .filter(|sample| sample.flags.contains(&TrajectoryFlag::MachTransition))
8014 .map(|sample| sample.distance_m)
8015 .collect();
8016
8017 assert!(
8018 !flagged_distances.is_empty()
8019 && flagged_distances.iter().all(|distance| *distance > 0.0),
8020 "{solver_name} must flag genuine crossings only at positive range: {flagged_distances:?}"
8021 );
8022 }
8023 }
8024
8025 #[test]
8026 fn subsonic_launch_does_not_flag_a_muzzle_transition() {
8027 use crate::trajectory_sampling::TrajectoryFlag;
8028
8029 for (solver_name, use_rk4, use_adaptive_rk45) in [
8030 ("Euler", false, false),
8031 ("RK4", true, false),
8032 ("RK45", true, true),
8033 ] {
8034 let inputs = BallisticInputs {
8035 muzzle_velocity: 250.0,
8036 muzzle_angle: 0.02,
8037 enable_trajectory_sampling: true,
8038 sample_interval: 25.0,
8039 use_rk4,
8040 use_adaptive_rk45,
8041 ..BallisticInputs::default()
8042 };
8043 let mut solver = TrajectorySolver::new(
8044 inputs,
8045 WindConditions::default(),
8046 AtmosphericConditions::default(),
8047 );
8048 solver.set_max_range(300.0);
8049 let samples = solver
8050 .solve()
8051 .expect("subsonic solve should succeed")
8052 .sampled_points
8053 .expect("sampling was enabled");
8054
8055 assert!(
8056 samples
8057 .iter()
8058 .all(|sample| !sample.flags.contains(&TrajectoryFlag::MachTransition)),
8059 "{solver_name} marked a Mach transition for a launch already below Mach 1"
8060 );
8061 }
8062 }
8063
8064 #[test]
8065 fn mach_transition_tracker_requires_a_downward_crossing() {
8066 fn record(mach_values: &[f64]) -> Vec<f64> {
8067 let mut tracker = MachTransitionTracker::default();
8068 let mut distances = Vec::new();
8069 for (index, mach) in mach_values.iter().copied().enumerate() {
8070 tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
8071 }
8072 distances
8073 }
8074
8075 assert!(record(&[0.9, 0.8, 0.7]).is_empty());
8076 assert_eq!(record(&[1.1, 1.05, 0.99]), vec![20.0]);
8077 assert_eq!(record(&[1.2, 1.19, 1.0, 0.99]), vec![10.0, 30.0]);
8078 assert_eq!(record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]), vec![20.0, 30.0]);
8079 assert!(record(&[1.3, f64::NAN, 1.1]).is_empty());
8080 }
8081
8082 #[test]
8083 fn mach_transition_tracker_labels_0_9_without_touching_the_flat_vec() {
8084 fn record(mach_values: &[f64]) -> (Vec<f64>, MachTransitionTracker) {
8090 let mut tracker = MachTransitionTracker::default();
8091 let mut distances = Vec::new();
8092 for (index, mach) in mach_values.iter().copied().enumerate() {
8093 tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
8094 }
8095 (distances, tracker)
8096 }
8097
8098 let (distances, tracker) = record(&[0.9, 0.8, 0.7]);
8101 assert!(distances.is_empty()); assert_eq!(tracker.mach_1_2_distance_m, None);
8103 assert_eq!(tracker.mach_1_0_distance_m, None);
8104 assert_eq!(tracker.mach_0_9_distance_m, Some(10.0));
8105
8106 let (distances, tracker) = record(&[1.1, 1.05, 0.99]);
8108 assert_eq!(distances, vec![20.0]);
8109 assert_eq!(tracker.mach_1_2_distance_m, None);
8110 assert_eq!(tracker.mach_1_0_distance_m, Some(20.0));
8111 assert_eq!(tracker.mach_0_9_distance_m, None);
8112
8113 let (distances, tracker) = record(&[1.2, 1.19, 1.0, 0.99]);
8115 assert_eq!(distances, vec![10.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(10.0));
8117 assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
8118 assert_eq!(tracker.mach_0_9_distance_m, None);
8119
8120 let (distances, tracker) = record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]);
8123 assert_eq!(distances, vec![20.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(20.0));
8125 assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
8126 assert_eq!(tracker.mach_0_9_distance_m, Some(50.0));
8127 assert!(
8128 tracker.mach_1_2_distance_m < tracker.mach_1_0_distance_m
8129 && tracker.mach_1_0_distance_m < tracker.mach_0_9_distance_m,
8130 "labeled crossings must be strictly increasing downrange"
8131 );
8132
8133 let (distances, tracker) = record(&[1.3, f64::NAN, 1.1]);
8135 assert!(distances.is_empty());
8136 assert_eq!(tracker.mach_1_2_distance_m, None);
8137 assert_eq!(tracker.mach_1_0_distance_m, None);
8138 assert_eq!(tracker.mach_0_9_distance_m, None);
8139 }
8140
8141 #[test]
8142 fn humidity_percent_converts_and_clamps() {
8143 let mut i = BallisticInputs {
8145 humidity: 0.5,
8146 ..BallisticInputs::default()
8147 };
8148 assert!((i.humidity_percent() - 50.0).abs() < 1e-9, "0.5 -> 50%");
8149 i.humidity = 0.0;
8150 assert_eq!(i.humidity_percent(), 0.0);
8151 i.humidity = 1.0;
8152 assert_eq!(i.humidity_percent(), 100.0);
8153 i.humidity = 1.5; assert_eq!(i.humidity_percent(), 100.0);
8155 }
8156
8157 fn vertical_at(shot_azimuth: f64, range_m: f64) -> f64 {
8160 let inputs = BallisticInputs {
8161 muzzle_velocity: 800.0,
8162 bc_value: 0.5,
8163 bc_type: DragModel::G7,
8164 muzzle_angle: 0.02, enable_coriolis: true,
8166 latitude: Some(45.0),
8167 shot_azimuth,
8168 ground_threshold: f64::NEG_INFINITY, ..BallisticInputs::default()
8170 };
8171 let mut solver = TrajectorySolver::new(
8172 inputs,
8173 WindConditions::default(),
8174 AtmosphericConditions::default(),
8175 );
8176 solver.set_max_range(range_m + 50.0);
8177 let r = solver.solve().expect("solve");
8178 let pts = &r.points;
8179 for i in 1..pts.len() {
8180 if pts[i].position.x >= range_m {
8181 let p1 = &pts[i - 1];
8182 let p2 = &pts[i];
8183 let t = (range_m - p1.position.x) / (p2.position.x - p1.position.x);
8184 return p1.position.y + t * (p2.position.y - p1.position.y);
8185 }
8186 }
8187 panic!("range {range_m} not reached");
8188 }
8189
8190 #[test]
8195 fn eotvos_east_higher_than_west() {
8196 let range = 600.0;
8197 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!(
8201 east > west,
8202 "east ({east:.5}) must be higher than west ({west:.5}) at {range} m (Eötvös)"
8203 );
8204 assert!(
8205 east > north && north > west,
8206 "north ({north:.5}) must lie between east ({east:.5}) and west ({west:.5})"
8207 );
8208 assert!(
8209 (east - west) > 1e-3,
8210 "E-W vertical separation ({:.6} m) should be physically meaningful, not FP noise",
8211 east - west
8212 );
8213 }
8214
8215 #[test]
8223 fn labeled_mach_crossings_match_pinned_pre_change_flat_vec_across_solvers() {
8224 let cases = [
8226 ("Euler", false, false, 670.9878683238721_f64, 805.5274119916264_f64),
8227 ("RK4", true, false, 671.7257336844475_f64, 805.933409072171_f64),
8228 ("RK45", true, true, 672.4905711917901_f64, 806.5709746782849_f64),
8229 ];
8230
8231 for (solver_name, use_rk4, use_adaptive_rk45, expected_1_2, expected_1_0) in cases {
8232 let inputs = BallisticInputs {
8233 muzzle_velocity: 850.0,
8234 bc_value: 0.2,
8235 bc_type: DragModel::G7,
8236 muzzle_angle: 0.03,
8237 use_rk4,
8238 use_adaptive_rk45,
8239 ..BallisticInputs::default()
8240 };
8241 let mut solver = TrajectorySolver::new(
8242 inputs,
8243 WindConditions::default(),
8244 AtmosphericConditions::default(),
8245 );
8246 solver.set_max_range(2000.0);
8247 let result = solver.solve().expect("solve should succeed");
8248
8249 assert_eq!(
8250 result.mach_1_2_distance_m,
8251 Some(expected_1_2),
8252 "{solver_name}: mach_1_2_distance_m must match the pinned pre-change flat-Vec value"
8253 );
8254 assert_eq!(
8255 result.mach_1_0_distance_m,
8256 Some(expected_1_0),
8257 "{solver_name}: mach_1_0_distance_m must match the pinned pre-change flat-Vec value"
8258 );
8259
8260 let mach_1_2 = result.mach_1_2_distance_m.expect("crosses 1.2");
8261 let mach_1_0 = result.mach_1_0_distance_m.expect("crosses 1.0");
8262 let mach_0_9 = result
8263 .mach_0_9_distance_m
8264 .expect("this trajectory also goes past 0.9 within 2000 m");
8265 assert!(
8266 mach_1_2 < mach_1_0 && mach_1_0 < mach_0_9,
8267 "{solver_name}: labeled crossings must be strictly increasing downrange \
8268 (1.2={mach_1_2}, 1.0={mach_1_0}, 0.9={mach_0_9})"
8269 );
8270 }
8271 }
8272
8273 #[test]
8276 fn labeled_mach_crossings_are_none_for_a_fully_supersonic_trajectory() {
8277 for (solver_name, use_rk4, use_adaptive_rk45) in [
8278 ("Euler", false, false),
8279 ("RK4", true, false),
8280 ("RK45", true, true),
8281 ] {
8282 let inputs = BallisticInputs {
8283 muzzle_velocity: 850.0,
8284 bc_value: 0.2,
8285 bc_type: DragModel::G7,
8286 muzzle_angle: 0.03,
8287 use_rk4,
8288 use_adaptive_rk45,
8289 ..BallisticInputs::default()
8290 };
8291 let mut solver = TrajectorySolver::new(
8292 inputs,
8293 WindConditions::default(),
8294 AtmosphericConditions::default(),
8295 );
8296 solver.set_max_range(200.0);
8298 let result = solver.solve().expect("solve should succeed");
8299
8300 assert_eq!(
8301 result.mach_1_2_distance_m, None,
8302 "{solver_name}: must not report a 1.2 crossing that never happens"
8303 );
8304 assert_eq!(
8305 result.mach_1_0_distance_m, None,
8306 "{solver_name}: must not report a 1.0 crossing that never happens"
8307 );
8308 assert_eq!(
8309 result.mach_0_9_distance_m, None,
8310 "{solver_name}: must not report a 0.9 crossing that never happens"
8311 );
8312 }
8313 }
8314}
8315
8316#[cfg(test)]
8317mod cant_tests {
8318 use super::*;
8319
8320 fn base_inputs() -> BallisticInputs {
8321 BallisticInputs {
8322 muzzle_velocity: 800.0,
8323 bc_value: 0.5,
8324 bc_type: DragModel::G7,
8325 bullet_mass: 0.0109,
8326 bullet_diameter: 0.00782,
8327 bullet_length: 0.0309,
8328 sight_height: 0.05,
8329 twist_rate: 10.0,
8330 use_rk4: true,
8331 ..BallisticInputs::default()
8332 }
8333 }
8334
8335 fn solve_with(inputs: BallisticInputs, max_range: f64) -> TrajectoryResult {
8336 let mut s = TrajectorySolver::new(
8337 inputs,
8338 WindConditions::default(),
8339 AtmosphericConditions::default(),
8340 );
8341 s.set_max_range(max_range);
8342 s.solve().expect("solve")
8343 }
8344
8345 fn yz_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
8347 let pts = &result.points;
8348 for i in 1..pts.len() {
8349 if pts[i].position.x >= x {
8350 let (p1, p2) = (&pts[i - 1], &pts[i]);
8351 let dx = p2.position.x - p1.position.x;
8352 let t = if dx.abs() < 1e-12 { 0.0 } else { (x - p1.position.x) / dx };
8353 return (
8354 p1.position.y + t * (p2.position.y - p1.position.y),
8355 p1.position.z + t * (p2.position.z - p1.position.z),
8356 );
8357 }
8358 }
8359 panic!("trajectory never reached {x} m");
8360 }
8361
8362 #[test]
8363 fn cant_sign_clockwise_up_offset_goes_right_and_low() {
8364 let mut level = base_inputs();
8366 level.muzzle_angle = 0.003; let mut canted = level.clone();
8368 canted.cant_angle = 10f64.to_radians();
8369
8370 let (y0, z0) = yz_at(&solve_with(level, 400.0), 300.0);
8371 let (y1, z1) = yz_at(&solve_with(canted, 400.0), 300.0);
8372 assert!(z1 > z0 + 0.01, "clockwise cant must move POI right: z0={z0} z1={z1}");
8373 assert!(y1 < y0 - 0.001, "clockwise cant must move POI low: y0={y0} y1={y1}");
8374 }
8375
8376 #[test]
8377 fn pure_cant_shows_bore_offset_near_range() {
8378 let mut i = base_inputs();
8381 i.muzzle_angle = 0.0;
8382 i.cant_angle = 10f64.to_radians();
8383 let sh = i.sight_height;
8384 let r = solve_with(i, 60.0);
8385 let first = &r.points[1]; let expected = -sh * 10f64.to_radians().sin();
8387 assert!(
8388 (first.position.z - expected).abs() < 0.005,
8389 "near-muzzle lateral {} should be ~bore offset {expected}",
8390 first.position.z
8391 );
8392 }
8393
8394 #[test]
8395 fn zero_angle_is_independent_of_cant() {
8396 let a = base_inputs();
8397 let mut b = base_inputs();
8398 b.cant_angle = 15f64.to_radians();
8399 let za = calculate_zero_angle(a.clone(), 100.0, 0.0).expect("zero a");
8400 let zb = calculate_zero_angle(b.clone(), 100.0, 0.0).expect("zero b");
8401 assert_eq!(za.to_bits(), zb.to_bits(), "zeroing must ignore cant: {za} vs {zb}");
8402 let _ = (a.cant_angle, b.cant_angle);
8404 }
8405
8406 #[test]
8407 fn nonfinite_cant_is_rejected() {
8408 let mut i = base_inputs();
8409 i.cant_angle = f64::NAN;
8410 let s = TrajectorySolver::new(i, WindConditions::default(), AtmosphericConditions::default());
8411 assert!(s.solve().is_err());
8412 }
8413
8414 #[test]
8415 fn incline_and_cant_compose_without_breaking() {
8416 let mut flat = base_inputs();
8418 flat.muzzle_angle = 0.003;
8419 flat.shooting_angle = 15f64.to_radians();
8420 let mut canted = flat.clone();
8421 canted.cant_angle = 10f64.to_radians();
8422 let (_, z_flat) = yz_at(&solve_with(flat, 400.0), 300.0);
8423 let (_, z_cant) = yz_at(&solve_with(canted, 400.0), 300.0);
8424 assert!(z_cant > z_flat, "cant must still deflect right on an incline");
8425 }
8426}
8427
8428#[cfg(test)]
8429mod vertical_wind_tests {
8430 use super::*;
8431
8432 fn base_inputs() -> BallisticInputs {
8433 BallisticInputs {
8434 muzzle_velocity: 800.0,
8435 bc_value: 0.5,
8436 bc_type: DragModel::G7,
8437 bullet_mass: 0.0109,
8438 bullet_diameter: 0.00782,
8439 bullet_length: 0.0309,
8440 sight_height: 0.05,
8441 twist_rate: 10.0,
8442 use_rk4: true,
8443 ..BallisticInputs::default()
8444 }
8445 }
8446
8447 fn y_at(result: &TrajectoryResult, x: f64) -> f64 {
8449 let pts = &result.points;
8450 for i in 1..pts.len() {
8451 if pts[i].position.x >= x {
8452 let (p1, p2) = (&pts[i - 1], &pts[i]);
8453 let dx = p2.position.x - p1.position.x;
8454 let t = if dx.abs() < 1e-12 { 0.0 } else { (x - p1.position.x) / dx };
8455 return p1.position.y + t * (p2.position.y - p1.position.y);
8456 }
8457 }
8458 panic!("trajectory never reached {x} m");
8459 }
8460
8461 fn solve_with(inputs: BallisticInputs, wind: WindConditions, max_range: f64) -> TrajectoryResult {
8462 let mut s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
8463 s.set_max_range(max_range);
8464 s.solve().expect("solve")
8465 }
8466
8467 #[test]
8468 fn updraft_raises_poi_downrange() {
8469 let calm_inputs = base_inputs();
8472 let calm_wind = WindConditions::default();
8473 let updraft = WindConditions {
8474 vertical_speed: 5.0,
8475 ..Default::default()
8476 };
8477
8478 let calm = solve_with(calm_inputs.clone(), calm_wind, 500.0);
8479 let updraft_result = solve_with(calm_inputs, updraft, 500.0);
8480
8481 let y_calm = y_at(&calm, 400.0);
8482 let y_updraft = y_at(&updraft_result, 400.0);
8483 assert!(
8484 y_updraft > y_calm,
8485 "5 m/s updraft must raise POI at 400m: calm={y_calm}, updraft={y_updraft}"
8486 );
8487 }
8488
8489 #[test]
8490 fn zero_vertical_is_default_and_finite_required() {
8491 assert_eq!(WindConditions::default().vertical_speed, 0.0);
8492
8493 let inputs = base_inputs();
8494 let wind = WindConditions {
8495 vertical_speed: f64::NAN,
8496 ..Default::default()
8497 };
8498 let s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
8499 assert!(
8500 s.solve().is_err(),
8501 "NaN wind.vertical_speed must be rejected by validate_for_solve"
8502 );
8503 }
8504}
8505
8506#[cfg(test)]
8508mod bc_reference_standard_tests {
8509 use super::*;
8510
8511 fn base_inputs() -> BallisticInputs {
8512 BallisticInputs {
8513 muzzle_velocity: 800.0,
8514 bc_value: 0.5,
8515 bc_type: DragModel::G7,
8516 bullet_mass: 0.0109,
8517 bullet_diameter: 0.00782,
8518 bullet_length: 0.0309,
8519 sight_height: 0.05,
8520 twist_rate: 10.0,
8521 use_rk4: true,
8522 ..BallisticInputs::default()
8523 }
8524 }
8525
8526 fn y_and_speed_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
8528 let pts = &result.points;
8529 for i in 1..pts.len() {
8530 if pts[i].position.x >= x {
8531 let (p1, p2) = (&pts[i - 1], &pts[i]);
8532 let dx = p2.position.x - p1.position.x;
8533 let t = if dx.abs() < 1e-12 {
8534 0.0
8535 } else {
8536 (x - p1.position.x) / dx
8537 };
8538 return (
8539 p1.position.y + t * (p2.position.y - p1.position.y),
8540 p1.velocity_magnitude + t * (p2.velocity_magnitude - p1.velocity_magnitude),
8541 );
8542 }
8543 }
8544 panic!("trajectory never reached {x} m");
8545 }
8546
8547 #[test]
8550 fn asm_to_icao_ratio_matches_documented_value() {
8551 assert!(
8552 (crate::constants::ASM_TO_ICAO_BC - 0.98237).abs() < 1e-5,
8553 "ASM_TO_ICAO_BC = {} must equal 0.98237 to 5 decimal places",
8554 crate::constants::ASM_TO_ICAO_BC
8555 );
8556 assert_eq!(
8561 crate::constants::ASM_TO_ICAO_BC,
8562 crate::constants::ASM_DENSITY_LB_FT3 / crate::constants::ICAO_DENSITY_LB_FT3
8563 );
8564 }
8565
8566 #[test]
8569 fn default_bc_reference_standard_is_icao() {
8570 assert_eq!(
8571 BallisticInputs::default().bc_reference_standard,
8572 BcReferenceStandard::Icao
8573 );
8574 }
8575
8576 #[test]
8581 fn icao_reference_leaves_bc_value_bit_identical() {
8582 let raw_bc: f64 = 0.4372911; let inputs = BallisticInputs {
8584 bc_value: raw_bc,
8585 bc_reference_standard: BcReferenceStandard::Icao,
8586 ..base_inputs()
8587 };
8588 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8589 assert_eq!(solver.inputs.bc_value.to_bits(), raw_bc.to_bits());
8590 }
8591
8592 #[test]
8593 fn default_inputs_solve_is_unaffected_by_the_new_field_existing() {
8594 let a = TrajectorySolver::new(base_inputs(), WindConditions::default(), AtmosphericConditions::default())
8598 .solve()
8599 .expect("solve a");
8600 let b = TrajectorySolver::new(
8601 BallisticInputs { ..base_inputs() },
8602 WindConditions::default(),
8603 AtmosphericConditions::default(),
8604 )
8605 .solve()
8606 .expect("solve b");
8607 assert_eq!(a.impact_velocity.to_bits(), b.impact_velocity.to_bits());
8608 assert_eq!(a.max_range.to_bits(), b.max_range.to_bits());
8609 }
8610
8611 #[test]
8614 fn army_standard_metro_scales_bc_value_by_exactly_the_derived_ratio() {
8615 let raw_bc = 0.5;
8616 let inputs = BallisticInputs {
8617 bc_value: raw_bc,
8618 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8619 ..base_inputs()
8620 };
8621 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8622 assert_eq!(
8623 solver.inputs.bc_value,
8624 raw_bc * crate::constants::ASM_TO_ICAO_BC
8625 );
8626 }
8627
8628 #[test]
8629 fn army_standard_metro_scales_mach_keyed_bc_segments() {
8630 let inputs = BallisticInputs {
8631 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8632 bc_segments: Some(vec![(0.5, 0.40), (1.5, 0.30)]),
8633 ..base_inputs()
8634 };
8635 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8636 let segments = solver.inputs.bc_segments.as_ref().expect("segments");
8637 assert_eq!(segments[0], (0.5, 0.40 * crate::constants::ASM_TO_ICAO_BC));
8638 assert_eq!(segments[1], (1.5, 0.30 * crate::constants::ASM_TO_ICAO_BC));
8639 }
8640
8641 #[test]
8642 fn army_standard_metro_scales_velocity_keyed_bc_segments_data() {
8643 let inputs = BallisticInputs {
8644 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8645 bc_segments_data: Some(vec![
8646 crate::BCSegmentData {
8647 velocity_min: 0.0,
8648 velocity_max: 500.0,
8649 bc_value: 0.40,
8650 },
8651 crate::BCSegmentData {
8652 velocity_min: 500.0,
8653 velocity_max: 900.0,
8654 bc_value: 0.45,
8655 },
8656 ]),
8657 ..base_inputs()
8658 };
8659 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8660 let segments = solver.inputs.bc_segments_data.as_ref().expect("segments");
8661 assert_eq!(segments[0].bc_value, 0.40 * crate::constants::ASM_TO_ICAO_BC);
8662 assert_eq!(segments[1].bc_value, 0.45 * crate::constants::ASM_TO_ICAO_BC);
8663 assert_eq!(segments[0].velocity_min, 0.0);
8665 assert_eq!(segments[1].velocity_max, 900.0);
8666 }
8667
8668 #[test]
8673 fn army_standard_metro_moves_impact_in_the_more_drag_direction() {
8674 let solve_at = |standard: BcReferenceStandard| {
8675 let inputs = BallisticInputs {
8676 bc_value: 0.475,
8677 bc_reference_standard: standard,
8678 ..base_inputs()
8679 };
8680 let mut solver =
8681 TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8682 solver.set_max_range(500.0);
8683 solver.solve().expect("solve")
8684 };
8685
8686 let icao = solve_at(BcReferenceStandard::Icao);
8687 let asm = solve_at(BcReferenceStandard::ArmyStandardMetro);
8688
8689 let (y_icao, v_icao) = y_and_speed_at(&icao, 400.0);
8690 let (y_asm, v_asm) = y_and_speed_at(&asm, 400.0);
8691
8692 assert!(
8693 y_asm < y_icao,
8694 "ArmyStandardMetro must drop MORE (lower y) at 400m than Icao for the same raw \
8695 bc_value: icao_y={y_icao}, asm_y={y_asm}"
8696 );
8697 assert!(
8698 v_asm < v_icao,
8699 "ArmyStandardMetro must retain LESS velocity at 400m than Icao for the same raw \
8700 bc_value: icao_v={v_icao}, asm_v={v_asm}"
8701 );
8702 }
8703
8704 #[test]
8711 fn monte_carlo_inherits_the_normalized_bc_reference() {
8712 let base_inputs_asm = BallisticInputs {
8713 bc_value: 0.475,
8714 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8715 ..base_inputs()
8716 };
8717 let wind = WindConditions::default();
8718
8719 let mut direct_solver =
8725 TrajectorySolver::new(base_inputs_asm.clone(), wind.clone(), AtmosphericConditions::default());
8726 direct_solver.set_max_range(base_inputs_asm.target_distance.max(1000.0) * 2.0);
8727 let direct = direct_solver.solve().expect("direct solve");
8728
8729 let mc_params = MonteCarloParams {
8730 num_simulations: 1,
8731 velocity_std_dev: 0.0,
8732 angle_std_dev: 0.0,
8733 bc_std_dev: 0.0,
8734 wind_speed_std_dev: 0.0,
8735 target_distance: None,
8736 base_wind_speed: 0.0,
8737 base_wind_direction: 0.0,
8738 azimuth_std_dev: 0.0,
8739 };
8740 let mc = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
8741 base_inputs_asm,
8742 wind,
8743 mc_params,
8744 0.0,
8745 42,
8746 )
8747 .expect("monte carlo");
8748
8749 assert_eq!(mc.ranges.len(), 1);
8750 assert_eq!(
8751 mc.ranges[0].to_bits(),
8752 direct.max_range.to_bits(),
8753 "a zero-dispersion single MC sample must match a plain solve of the same \
8754 ASM-referenced inputs bit-for-bit"
8755 );
8756 assert_eq!(
8757 mc.impact_velocities[0].to_bits(),
8758 direct.impact_velocity.to_bits()
8759 );
8760 }
8761
8762 #[test]
8770 fn estimate_bc_fit_recovers_an_icao_referenced_bc() {
8771 let known_bc = 0.475;
8772 let velocity = 800.0;
8773 let mass = 0.0109;
8774 let diameter = 0.00782;
8775 let atmosphere = AtmosphericConditions::default();
8776
8777 let synth_inputs = BallisticInputs {
8778 muzzle_velocity: velocity,
8779 bc_value: known_bc,
8780 bc_type: DragModel::G7,
8781 bullet_mass: mass,
8782 bullet_diameter: diameter,
8783 bullet_length: 0.0309,
8784 sight_height: 0.05,
8785 twist_rate: 10.0,
8786 use_rk4: true,
8787 bc_reference_standard: BcReferenceStandard::Icao,
8788 ..BallisticInputs::default()
8789 };
8790 let mut solver = TrajectorySolver::new(synth_inputs, WindConditions::default(), atmosphere.clone());
8791 solver.set_max_range(500.0);
8792 let trajectory = solver.solve().expect("synthetic solve");
8793
8794 let points: Vec<(f64, f64)> = [100.0, 200.0, 300.0, 400.0]
8795 .iter()
8796 .map(|&d| {
8797 let (y, _) = {
8798 let pts = &trajectory.points;
8799 let mut found = None;
8800 for i in 1..pts.len() {
8801 if pts[i].position.x >= d {
8802 let (p1, p2) = (&pts[i - 1], &pts[i]);
8803 let dx = p2.position.x - p1.position.x;
8804 let t = if dx.abs() < 1e-12 {
8805 0.0
8806 } else {
8807 (d - p1.position.x) / dx
8808 };
8809 found = Some((
8810 p1.position.y + t * (p2.position.y - p1.position.y),
8811 0.0,
8812 ));
8813 break;
8814 }
8815 }
8816 found.expect("trajectory reached observation distance")
8817 };
8818 (d, -y) })
8820 .collect();
8821
8822 let estimate = estimate_bc_fit(
8823 velocity,
8824 mass,
8825 diameter,
8826 &points,
8827 DragModel::G7,
8828 BcFitMode::Drop,
8829 atmosphere,
8830 None,
8831 0.05,
8832 )
8833 .expect("fit should converge");
8834
8835 assert!(
8836 (estimate.bc - known_bc).abs() < 0.02,
8837 "fit should recover the known ICAO-referenced bc={known_bc}, got {}",
8838 estimate.bc
8839 );
8840 }
8841
8842 #[test]
8845 fn custom_drag_table_makes_bc_reference_standard_numerically_inert() {
8846 let table = crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0, 3.0], vec![0.3, 0.4, 0.3, 0.2])
8847 .expect("valid table");
8848
8849 let solve_with = |standard: BcReferenceStandard| {
8850 let inputs = BallisticInputs {
8851 bc_value: 0.5, bc_reference_standard: standard,
8853 custom_drag_table: Some(table.clone()),
8854 ..base_inputs()
8855 };
8856 let mut solver =
8857 TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8858 solver.set_max_range(500.0);
8859 solver.solve().expect("solve")
8860 };
8861
8862 let icao = solve_with(BcReferenceStandard::Icao);
8863 let asm = solve_with(BcReferenceStandard::ArmyStandardMetro);
8864
8865 assert_eq!(
8866 icao.impact_velocity.to_bits(),
8867 asm.impact_velocity.to_bits(),
8868 "a custom drag table must make bc_reference_standard fully inert"
8869 );
8870 assert_eq!(icao.max_range.to_bits(), asm.max_range.to_bits());
8871 }
8872
8873 #[test]
8874 fn custom_drag_table_inert_warning_fires_only_for_army_standard_metro_with_a_table() {
8875 let table = crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0], vec![0.3, 0.4, 0.3])
8876 .expect("valid table");
8877
8878 let no_table_icao = base_inputs();
8880 assert!(no_table_icao.bc_reference_standard_inert_warning().is_none());
8881 let no_table_asm = BallisticInputs {
8882 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8883 ..base_inputs()
8884 };
8885 assert!(no_table_asm.bc_reference_standard_inert_warning().is_none());
8886
8887 let table_icao = BallisticInputs {
8889 custom_drag_table: Some(table.clone()),
8890 ..base_inputs()
8891 };
8892 assert!(table_icao.bc_reference_standard_inert_warning().is_none());
8893
8894 let table_asm = BallisticInputs {
8896 custom_drag_table: Some(table),
8897 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8898 ..base_inputs()
8899 };
8900 let warning = table_asm
8901 .bc_reference_standard_inert_warning()
8902 .expect("must warn");
8903 assert!(warning.contains("--bc-reference"));
8904 assert!(warning.contains("--drag-table"));
8905 }
8906}
8907
8908#[cfg(test)]
8914mod effective_drag_coefficient_tests {
8915 use super::*;
8916
8917 fn inputs_175gr_g7() -> BallisticInputs {
8918 let mut inputs = BallisticInputs {
8919 bc_value: 0.243,
8920 bc_type: DragModel::G7,
8921 muzzle_velocity: 823.0,
8922 ..Default::default()
8923 };
8924 inputs.bullet_mass = 175.0 * crate::constants::GRAINS_TO_KG;
8928 inputs.bullet_diameter = 0.308 * 0.0254;
8929 inputs.weight_grains = 175.0;
8930 inputs.caliber_inches = 0.308;
8931 inputs
8932 }
8933
8934 fn solver(inputs: BallisticInputs) -> TrajectorySolver {
8935 TrajectorySolver::new(
8936 inputs,
8937 WindConditions::default(),
8938 AtmosphericConditions::default(),
8939 )
8940 }
8941
8942 #[test]
8946 fn reports_the_projectiles_own_cd_not_the_reference_tables() {
8947 let inputs = inputs_175gr_g7();
8948 let sd = inputs.sectional_density_lb_in2().expect("SD");
8949 let solver = solver(inputs);
8950
8951 let sos = 340.0;
8952 let velocity = 800.0;
8953 let mach = velocity / sos;
8954
8955 let reference = crate::drag::get_drag_coefficient(mach, &DragModel::G7);
8956 let reported = solver
8957 .effective_drag_coefficient(velocity, sos)
8958 .expect("mass and diameter are set");
8959
8960 let expected = reference * sd / 0.243;
8961 assert!(
8962 (reported - expected).abs() < 1e-12,
8963 "reported {reported} != Cd_ref * SD / BC {expected}"
8964 );
8965 assert!(
8968 (reported - reference).abs() > 1e-6,
8969 "form factor collapsed to 1; this fixture no longer distinguishes the two values"
8970 );
8971 }
8972
8973 #[test]
8976 fn a_custom_drag_table_passes_through_unscaled() {
8977 let mut inputs = inputs_175gr_g7();
8978 inputs.custom_drag_table = Some(crate::drag::DragTable::new(
8979 vec![0.5, 3.0],
8980 vec![0.15, 0.40],
8981 ));
8982 let solver = solver(inputs);
8983
8984 let sos = 340.0;
8985 let velocity = 0.9 * sos;
8986 let table_value = solver
8987 .inputs
8988 .custom_drag_table
8989 .as_ref()
8990 .expect("table")
8991 .interpolate(0.9);
8992
8993 let reported = solver
8994 .effective_drag_coefficient(velocity, sos)
8995 .expect("mass and diameter are set");
8996 assert!(
8997 (reported - table_value).abs() < 1e-12,
8998 "custom table Cd {table_value} was rescaled to {reported}"
8999 );
9000 }
9001
9002 #[test]
9005 fn a_velocity_segmented_bc_steps_the_reported_cd() {
9006 let mut inputs = inputs_175gr_g7();
9007 inputs.use_bc_segments = true;
9008 inputs.bc_segments_data = Some(vec![
9009 crate::BCSegmentData { velocity_min: 2400.0, velocity_max: 4000.0, bc_value: 0.243 },
9010 crate::BCSegmentData { velocity_min: 0.0, velocity_max: 2400.0, bc_value: 0.200 },
9011 ]);
9012 let solver = solver(inputs);
9013
9014 let sos = 340.0;
9015 let above = solver.effective_drag_coefficient(2500.0 / 3.28084, sos).expect("cd");
9017 let below = solver.effective_drag_coefficient(2300.0 / 3.28084, sos).expect("cd");
9018
9019 assert!(
9021 below > above,
9022 "expected the 0.200 band to report a higher Cd than the 0.243 band; got {below} vs {above}"
9023 );
9024 }
9025
9026 #[test]
9029 fn is_absent_when_sectional_density_is_unknown() {
9030 let mut inputs = inputs_175gr_g7();
9031 inputs.weight_grains = 0.0;
9032 inputs.bullet_mass = 0.0;
9033 let solver = solver(inputs);
9034 assert!(solver.effective_drag_coefficient(800.0, 340.0).is_none());
9035 }
9036
9037 #[test]
9042 fn the_json_emit_rule_is_flag_gated_and_absent_when_cd_is_unknown() {
9043 let mut point = TrajectoryPoint {
9044 time: 0.0,
9045 position: nalgebra::Vector3::new(0.0, 0.0, 0.0),
9046 velocity_magnitude: 800.0,
9047 kinetic_energy: 3000.0,
9048 drag_coefficient: Some(0.31),
9049 };
9050 assert_eq!(point.drag_coefficient_json_value(true), Some(0.31));
9051 assert_eq!(
9052 point.drag_coefficient_json_value(false),
9053 None,
9054 "without the flag the key must not exist, so default JSON stays byte-identical"
9055 );
9056 point.drag_coefficient = None;
9057 assert_eq!(
9058 point.drag_coefficient_json_value(true),
9059 None,
9060 "unknown sectional density must yield an ABSENT key, not null"
9061 );
9062 }
9063
9064 #[test]
9066 fn every_point_of_a_solved_trajectory_carries_the_value() {
9067 let mut solver = solver(inputs_175gr_g7());
9068 solver.set_max_range(300.0);
9069 let result = solver.solve().expect("solve");
9070 assert!(!result.points.is_empty());
9071 assert!(
9072 result.points.iter().all(|p| p.drag_coefficient.is_some()),
9073 "the post-integration pass missed at least one point"
9074 );
9075 }
9076}