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 cd_delta2: f64,
254 pub use_bc_segments: bool,
257 pub bc_segments: Option<Vec<(f64, f64)>>, pub bc_segments_data: Option<Vec<crate::BCSegmentData>>, pub use_enhanced_spin_drift: bool,
260 pub use_form_factor: bool,
263 pub enable_wind_shear: bool,
264 pub wind_shear_model: String,
265 pub enable_trajectory_sampling: bool,
266 pub sample_interval: f64, pub drops_reference: DropsReference,
276 pub enable_pitch_damping: bool,
277 pub enable_precession_nutation: bool,
278 pub enable_aerodynamic_jump: bool,
281 pub use_cluster_bc: bool, pub custom_drag_table: Option<crate::drag::DragTable>,
285 pub cd_scale: f64,
293
294 pub bc_type_str: Option<String>,
296}
297
298impl BallisticInputs {
299 pub fn humidity_percent(&self) -> f64 {
304 (self.humidity * 100.0).clamp(0.0, 100.0)
305 }
306
307 pub fn windage_zero_bias_rad(&self, zero_distance_m: f64) -> f64 {
322 if zero_distance_m > 0.0 {
323 (self.zero_poi_horizontal_m + self.sight_offset_lateral_m) / zero_distance_m
324 } else {
325 0.0
326 }
327 }
328
329 pub fn sectional_density_lb_in2(&self) -> Option<f64> {
335 let weight_gr = if self.weight_grains > 0.0 {
336 self.weight_grains
337 } else {
338 self.bullet_mass / crate::constants::GRAINS_TO_KG };
340 let diameter_in = if self.caliber_inches > 0.0 {
341 self.caliber_inches
342 } else {
343 self.bullet_diameter / 0.0254 };
345 if weight_gr > 0.0 && diameter_in > 0.0 {
346 Some(weight_gr / 7000.0 / (diameter_in * diameter_in))
347 } else {
348 None
349 }
350 }
351
352 pub fn custom_drag_denominator(&self, fallback_bc: f64) -> f64 {
364 match self.sectional_density_lb_in2() {
365 Some(sd) => sd,
366 None => {
367 static WARN_ONCE: std::sync::Once = std::sync::Once::new();
368 WARN_ONCE.call_once(|| {
369 eprintln!(
370 "Warning: custom drag table active but bullet mass/diameter are \
371 unavailable; falling back to bc_value for the retardation denominator"
372 );
373 });
374 fallback_bc
375 }
376 }
377 }
378
379 pub fn bc_reference_standard_inert_warning(&self) -> Option<&'static str> {
390 if self.custom_drag_table.is_some()
391 && matches!(self.bc_reference_standard, BcReferenceStandard::ArmyStandardMetro)
392 {
393 Some(BC_REFERENCE_STANDARD_INERT_WARNING)
394 } else {
395 None
396 }
397 }
398
399 pub fn normalize_for_solve(&mut self) {
418 if matches!(
429 self.bc_reference_standard,
430 BcReferenceStandard::ArmyStandardMetro
431 ) {
432 self.bc_value *= crate::constants::ASM_TO_ICAO_BC;
433 if let Some(segments) = self.bc_segments.as_mut() {
434 for (_mach, bc) in segments.iter_mut() {
435 *bc *= crate::constants::ASM_TO_ICAO_BC;
436 }
437 }
438 if let Some(segments) = self.bc_segments_data.as_mut() {
439 for segment in segments.iter_mut() {
440 segment.bc_value *= crate::constants::ASM_TO_ICAO_BC;
441 }
442 }
443 self.bc_reference_standard = BcReferenceStandard::Icao;
446 }
447
448 self.caliber_inches = self.bullet_diameter / 0.0254;
453 self.weight_grains = self.bullet_mass / crate::constants::GRAINS_TO_KG;
454
455 self.muzzle_velocity = resolve_powder_adjusted_velocity(
466 self.muzzle_velocity,
467 self.temperature,
468 self.use_powder_sensitivity,
469 self.powder_temp_sensitivity,
470 self.powder_temp,
471 self.powder_temp_curve.as_deref(),
472 self.powder_curve_temp_c,
473 );
474 }
475}
476
477impl Default for BallisticInputs {
478 fn default() -> Self {
479 let mass_kg = 0.01;
480 let diameter_m = 0.00762;
481 let bc = 0.5;
482 let muzzle_angle_rad = 0.0;
483 let bc_type = DragModel::G1;
484
485 Self {
486 bc_value: bc,
488 bc_type,
489 bc_reference_standard: BcReferenceStandard::Icao,
490 bullet_mass: mass_kg,
491 muzzle_velocity: 800.0,
492 bullet_diameter: diameter_m,
493 bullet_length: crate::stability::estimate_bullet_length_m(diameter_m, mass_kg),
497
498 muzzle_angle: muzzle_angle_rad,
500 target_distance: 100.0,
501 azimuth_angle: 0.0,
502 shot_azimuth: 0.0,
503 shooting_angle: 0.0,
504 cant_angle: 0.0,
505 sight_height: 0.05,
506 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,
511 ground_threshold: -100.0, altitude: 0.0,
515 temperature: 15.0,
516 pressure: 1013.25, humidity: 0.5, latitude: None,
519
520 wind_speed: 0.0,
522 wind_angle: 0.0,
523
524 twist_rate: 12.0, is_twist_right: true,
527 caliber_inches: diameter_m / 0.0254, weight_grains: mass_kg / crate::constants::GRAINS_TO_KG, manufacturer: None,
530 bullet_model: None,
531 bullet_id: None,
532 bullet_cluster: None,
533
534 use_rk4: true, use_adaptive_rk45: true, enable_advanced_effects: false,
540 enable_magnus: false,
541 enable_coriolis: false,
542 use_powder_sensitivity: false,
543 powder_temp_sensitivity: 0.0,
544 powder_temp: 15.0,
545 powder_temp_curve: None,
546 powder_curve_temp_c: None,
547 tipoff_yaw: 0.0,
548 tipoff_decay_distance: 50.0,
549 cd_delta2: 7.5,
550 use_bc_segments: false,
551 bc_segments: None,
552 bc_segments_data: None,
553 use_enhanced_spin_drift: false,
554 use_form_factor: false,
555 enable_wind_shear: false,
556 wind_shear_model: "none".to_string(),
557 enable_trajectory_sampling: false,
558 sample_interval: 10.0, drops_reference: DropsReference::Los, enable_pitch_damping: false,
561 enable_precession_nutation: false,
562 enable_aerodynamic_jump: false,
563 use_cluster_bc: false, custom_drag_table: None,
567 cd_scale: 1.0,
568
569 bc_type_str: None,
571 }
572 }
573}
574
575pub fn interpolate_powder_temp_curve(curve: &[(f64, f64)], temp_c: f64) -> f64 {
581 debug_assert!(!curve.is_empty());
582 if curve.is_empty() {
583 return 0.0;
584 }
585 let mut sorted;
588 let pts: &[(f64, f64)] = if curve.windows(2).all(|w| w[0].0 <= w[1].0) {
589 curve
590 } else {
591 sorted = curve.to_vec();
592 sorted.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
593 &sorted
594 };
595 let n = pts.len();
596 if temp_c <= pts[0].0 {
597 return pts[0].1; }
599 if temp_c >= pts[n - 1].0 {
600 return pts[n - 1].1; }
602 for i in 1..n {
603 let (t0, v0) = pts[i - 1];
604 let (t1, v1) = pts[i];
605 if temp_c <= t1 {
606 let span = t1 - t0;
607 if span.abs() < f64::EPSILON {
608 return v1; }
610 let f = (temp_c - t0) / span;
611 return v0 + f * (v1 - v0);
612 }
613 }
614 pts[n - 1].1
615}
616
617pub fn parse_powder_sweep(s: &str) -> Result<Vec<f64>, String> {
623 const MAX_SWEEP_ROWS: usize = 500;
624 let parts: Vec<&str> = s.split(':').collect();
625 if parts.len() != 3 {
626 return Err(format!(
627 "Invalid --sweep '{}': expected START:END:STEP (e.g. \"20:110:10\")",
628 s
629 ));
630 }
631 let parse = |p: &str, name: &str| -> Result<f64, String> {
632 p.trim()
633 .parse::<f64>()
634 .map_err(|_| format!("Invalid --sweep {}: '{}' is not a number", name, p.trim()))
635 };
636 let start = parse(parts[0], "START")?;
637 let end = parse(parts[1], "END")?;
638 let step = parse(parts[2], "STEP")?;
639 if !step.is_finite() || step <= 0.0 {
640 return Err(format!("Invalid --sweep STEP {}: must be positive", step));
641 }
642 if !start.is_finite() || !end.is_finite() || end < start {
643 return Err(format!(
644 "Invalid --sweep range {}:{}: END must be >= START",
645 start, end
646 ));
647 }
648 let n_f = ((end - start) / step + 1e-9).floor();
654 if !n_f.is_finite() || n_f + 1.0 > MAX_SWEEP_ROWS as f64 {
655 return Err(format!(
656 "--sweep would produce more than {} rows; use a larger STEP",
657 MAX_SWEEP_ROWS
658 ));
659 }
660 let n = n_f as usize + 1;
661 Ok((0..n).map(|i| start + step * i as f64).collect())
663}
664
665pub fn resolve_powder_adjusted_velocity(
674 nominal_velocity_mps: f64,
675 ambient_temperature_c: f64,
676 use_powder_sensitivity: bool,
677 powder_temp_sensitivity_mps_per_c: f64,
678 powder_reference_temp_c: f64,
679 powder_temp_curve: Option<&[(f64, f64)]>,
680 powder_curve_temp_c: Option<f64>,
681) -> f64 {
682 if let Some(curve) = powder_temp_curve {
683 if !curve.is_empty() {
684 let lookup_c = powder_curve_temp_c.unwrap_or(ambient_temperature_c);
685 return interpolate_powder_temp_curve(curve, lookup_c);
686 }
687 return nominal_velocity_mps;
690 }
691 if use_powder_sensitivity {
692 let temp_delta_c = ambient_temperature_c - powder_reference_temp_c;
693 return nominal_velocity_mps + powder_temp_sensitivity_mps_per_c * temp_delta_c;
694 }
695 nominal_velocity_mps
696}
697
698#[derive(Debug, Clone)]
700pub struct WindConditions {
701 pub speed: f64, pub direction: f64,
705 pub vertical_speed: f64,
713}
714
715impl Default for WindConditions {
716 fn default() -> Self {
717 Self {
718 speed: 0.0,
719 direction: 0.0,
720 vertical_speed: 0.0,
721 }
722 }
723}
724
725#[derive(Debug, Clone)]
727pub struct AtmosphericConditions {
728 pub temperature: f64, pub pressure: f64, pub humidity: f64,
734 pub altitude: f64, }
736
737impl Default for AtmosphericConditions {
738 fn default() -> Self {
739 Self {
740 temperature: 15.0,
741 pressure: 1013.25,
742 humidity: 50.0,
743 altitude: 0.0,
744 }
745 }
746}
747
748#[derive(Debug, Clone)]
750pub struct TrajectoryPoint {
751 pub time: f64,
752 pub position: Vector3<f64>,
753 pub velocity_magnitude: f64,
754 pub kinetic_energy: f64,
755 pub drag_coefficient: Option<f64>,
762}
763
764impl TrajectoryPoint {
765 pub fn drag_coefficient_json_value(&self, with_drag_coefficient: bool) -> Option<f64> {
773 if with_drag_coefficient {
774 self.drag_coefficient
775 } else {
776 None
777 }
778 }
779}
780
781#[derive(Debug, Clone)]
783pub struct TrajectoryResult {
784 pub max_range: f64,
785 pub max_height: f64,
786 pub time_of_flight: f64,
787 pub impact_velocity: f64,
788 pub impact_energy: f64,
789 pub projectile_mass_kg: f64,
791 pub line_of_sight_height_m: f64,
793 pub station_speed_of_sound_mps: f64,
795 pub termination: TrajectoryTermination,
797 pub points: Vec<TrajectoryPoint>,
798 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>,
807 pub mach_1_2_distance_m: Option<f64>,
812 pub mach_1_0_distance_m: Option<f64>,
816 pub mach_0_9_distance_m: Option<f64>,
824}
825
826const RK45_TOLERANCE: f64 = 1e-6;
827const RK45_SAFETY_FACTOR: f64 = 0.9;
828const RK45_MAX_DT: f64 = 0.01;
829const RK45_MIN_DT: f64 = 1e-6;
830const TRAJECTORY_TIME_LIMIT_S: f64 = 100.0;
831
832pub const MAX_TRAJECTORY_POINTS: usize = 250_000;
838
839fn cli_rk45_error_norm(
841 position: &Vector3<f64>,
842 velocity: &Vector3<f64>,
843 fifth_position: &Vector3<f64>,
844 fifth_velocity: &Vector3<f64>,
845 fourth_position: &Vector3<f64>,
846 fourth_velocity: &Vector3<f64>,
847) -> f64 {
848 let pack_state = |position: &Vector3<f64>, velocity: &Vector3<f64>| {
849 Vector6::new(
850 position.x, position.y, position.z, velocity.x, velocity.y, velocity.z,
851 )
852 };
853 let state = pack_state(position, velocity);
854 let fifth_order = pack_state(fifth_position, fifth_velocity);
855 let fourth_order = pack_state(fourth_position, fourth_velocity);
856
857 crate::trajectory_integration::rk45_error_norm(&state, &fifth_order, &fourth_order)
858}
859
860struct Rk45Trial {
861 position: Vector3<f64>,
862 velocity: Vector3<f64>,
863 suggested_dt: f64,
864 error: f64,
865}
866
867struct Rk45AcceptedStep {
868 position: Vector3<f64>,
869 velocity: Vector3<f64>,
870 used_dt: f64,
871 next_dt: f64,
872 error: f64,
873}
874
875#[derive(Default)]
889struct MachTransitionTracker {
890 previous_mach: Option<f64>,
891 crossed_transonic: bool,
892 crossed_subsonic: bool,
893 crossed_narrow: bool,
894 mach_1_2_distance_m: Option<f64>,
897 mach_1_0_distance_m: Option<f64>,
900 mach_0_9_distance_m: Option<f64>,
903}
904
905impl MachTransitionTracker {
906 fn record_downward_crossings(&mut self, mach: f64, downrange_m: f64, distances: &mut Vec<f64>) {
907 if !mach.is_finite() {
908 self.previous_mach = None;
909 return;
910 }
911
912 if let Some(previous_mach) = self.previous_mach {
913 if !self.crossed_transonic && previous_mach >= 1.2 && mach < 1.2 {
914 self.crossed_transonic = true;
915 distances.push(downrange_m);
916 self.mach_1_2_distance_m = Some(downrange_m);
917 }
918 if !self.crossed_subsonic && previous_mach >= 1.0 && mach < 1.0 {
919 self.crossed_subsonic = true;
920 distances.push(downrange_m);
921 self.mach_1_0_distance_m = Some(downrange_m);
922 }
923 if !self.crossed_narrow && previous_mach >= 0.9 && mach < 0.9 {
924 self.crossed_narrow = true;
925 self.mach_0_9_distance_m = Some(downrange_m);
927 }
928 }
929 self.previous_mach = Some(mach);
930 }
931}
932
933impl TrajectoryResult {
934 pub fn position_at_range(&self, target_range: f64) -> Option<Vector3<f64>> {
938 if self.points.is_empty() {
939 return None;
940 }
941
942 for i in 0..self.points.len() - 1 {
944 let p1 = &self.points[i];
945 let p2 = &self.points[i + 1];
946
947 if p1.position.x <= target_range && p2.position.x >= target_range {
949 let dx = p2.position.x - p1.position.x;
951 if dx.abs() < 1e-10 {
952 return Some(p1.position);
953 }
954 let t = (target_range - p1.position.x) / dx;
955
956 return Some(Vector3::new(
958 target_range,
959 p1.position.y + t * (p2.position.y - p1.position.y),
960 p1.position.z + t * (p2.position.z - p1.position.z),
961 ));
962 }
963 }
964
965 self.points.last().map(|p| p.position)
967 }
968}
969
970#[derive(Debug, Clone, Copy, PartialEq, Eq)]
972enum StationAtmosphereResolution {
973 LegacyDefaultSentinels,
976 Authoritative,
979}
980
981#[derive(Clone)]
982pub struct TrajectorySolver {
983 inputs: BallisticInputs,
984 wind: WindConditions,
985 atmosphere: AtmosphericConditions,
986 station_atmosphere_resolution: StationAtmosphereResolution,
987 max_range: f64,
988 time_step: f64,
989 max_trajectory_points: usize,
990 cluster_bc: Option<ClusterBCDegradation>,
991 precession_nutation_inertias: (f64, f64),
993 wind_sock: Option<crate::wind::WindSock>,
998 atmo_sock: Option<crate::atmosphere::AtmoSock>,
1005}
1006
1007#[derive(Debug, Clone, Copy, PartialEq, Eq)]
1019pub(crate) enum ZeroTargetFrame {
1020 SightLine,
1021 WorldVertical,
1022}
1023
1024#[derive(Debug, Clone, Copy, PartialEq)]
1037pub struct ZeroCrossings {
1038 pub near_m: Option<f64>,
1041 pub far_m: Option<f64>,
1044}
1045
1046impl TrajectorySolver {
1047 pub fn new(
1048 inputs: BallisticInputs,
1049 wind: WindConditions,
1050 atmosphere: AtmosphericConditions,
1051 ) -> Self {
1052 Self::new_with_station_atmosphere_resolution(
1053 inputs,
1054 wind,
1055 atmosphere,
1056 StationAtmosphereResolution::LegacyDefaultSentinels,
1057 )
1058 }
1059
1060 pub fn new_with_resolved_station_atmosphere(
1071 inputs: BallisticInputs,
1072 wind: WindConditions,
1073 atmosphere: AtmosphericConditions,
1074 ) -> Self {
1075 Self::new_with_station_atmosphere_resolution(
1076 inputs,
1077 wind,
1078 atmosphere,
1079 StationAtmosphereResolution::Authoritative,
1080 )
1081 }
1082
1083 fn new_with_station_atmosphere_resolution(
1084 mut inputs: BallisticInputs,
1085 wind: WindConditions,
1086 atmosphere: AtmosphericConditions,
1087 station_atmosphere_resolution: StationAtmosphereResolution,
1088 ) -> Self {
1089 inputs.normalize_for_solve();
1094
1095 let cluster_bc = if inputs.use_cluster_bc {
1097 Some(ClusterBCDegradation::new())
1098 } else {
1099 None
1100 };
1101 let precession_nutation_inertias = projectile_moments_of_inertia(
1102 inputs.bullet_mass,
1103 inputs.bullet_diameter,
1104 inputs.bullet_length,
1105 );
1106
1107 Self {
1108 inputs,
1109 wind,
1110 atmosphere,
1111 station_atmosphere_resolution,
1112 max_range: 1000.0,
1113 time_step: 0.001,
1114 max_trajectory_points: MAX_TRAJECTORY_POINTS,
1115 cluster_bc,
1116 precession_nutation_inertias,
1117 wind_sock: None,
1118 atmo_sock: None,
1119 }
1120 }
1121
1122 pub fn set_max_range(&mut self, range: f64) {
1123 self.max_range = range;
1124 }
1125
1126 pub fn set_time_step(&mut self, step: f64) {
1127 self.time_step = step;
1128 }
1129
1130 pub(crate) fn calculate_and_set_zero_angle(
1134 &mut self,
1135 target_distance_m: f64,
1136 target_height_m: f64,
1137 frame: ZeroTargetFrame,
1138 ) -> Result<f64, BallisticsError> {
1139 let angle = self.find_zero_angle(target_distance_m, target_height_m, frame)?;
1140 let angle = if target_distance_m > 0.0 {
1149 angle + self.inputs.zero_poi_vertical_m / target_distance_m
1150 } else {
1151 angle
1152 };
1153 self.inputs.muzzle_angle = angle;
1154 self.apply_windage_zero_bias(target_distance_m);
1155 Ok(angle)
1156 }
1157
1158 pub(crate) fn apply_windage_zero_bias(&mut self, target_distance_m: f64) {
1172 self.inputs.azimuth_angle += self.inputs.windage_zero_bias_rad(target_distance_m);
1173 }
1174
1175 fn find_zero_angle(
1176 &self,
1177 target_distance_m: f64,
1178 target_height_m: f64,
1179 frame: ZeroTargetFrame,
1180 ) -> Result<f64, BallisticsError> {
1181 let mut low_angle = 0.0;
1184 let mut high_angle = 0.2; let tolerance = 1e-7;
1186 let max_iterations = 60;
1187
1188 let low_height = self.zero_trial_height_at(low_angle, target_distance_m, frame)?;
1190 let high_height = self.zero_trial_height_at(high_angle, target_distance_m, frame)?;
1191
1192 match (low_height, high_height) {
1193 (Some(low_height), Some(high_height)) => {
1194 let low_error = low_height - target_height_m;
1195 let high_error = high_height - target_height_m;
1196
1197 if low_error > 0.0 && high_error > 0.0 {
1198 } else if low_error < 0.0 && high_error < 0.0 {
1201 let mut expanded = false;
1203 for multiplier in [2.0, 3.0, 4.0] {
1204 let new_high = (high_angle * multiplier).min(0.785);
1205 if let Ok(Some(height)) =
1206 self.zero_trial_height_at(new_high, target_distance_m, frame)
1207 {
1208 if height - target_height_m > 0.0 {
1209 high_angle = new_high;
1210 expanded = true;
1211 break;
1212 }
1213 }
1214 if new_high >= 0.785 {
1215 break;
1216 }
1217 }
1218 if !expanded {
1219 return Err("Cannot find zero angle: target beyond effective range even at maximum angle".into());
1220 }
1221 }
1222 }
1223 (None, Some(_)) => {
1224 }
1227 (Some(_), None) => {
1228 return Err(
1229 "Cannot find zero angle: high angle trajectory doesn't reach target distance"
1230 .into(),
1231 );
1232 }
1233 (None, None) => {
1234 return Err(
1235 "Cannot find zero angle: trajectory cannot reach target distance at any angle"
1236 .into(),
1237 );
1238 }
1239 }
1240
1241 for _ in 0..max_iterations {
1242 let mid_angle = (low_angle + high_angle) / 2.0;
1243 match self.zero_trial_height_at(mid_angle, target_distance_m, frame)? {
1244 Some(height) => {
1245 let error = height - target_height_m;
1246 if error.abs() < 0.0001 {
1249 return Ok(mid_angle);
1250 }
1251
1252 if (high_angle - low_angle).abs() < tolerance {
1255 if error.abs() < 0.01 {
1256 return Ok(mid_angle);
1257 }
1258 return Err("Zero angle did not converge: residual height error too large (target not reachable / not bracketed)".into());
1259 }
1260
1261 if error > 0.0 {
1262 high_angle = mid_angle;
1263 } else {
1264 low_angle = mid_angle;
1265 }
1266 }
1267 None => {
1268 low_angle = mid_angle;
1269 if (high_angle - low_angle).abs() < tolerance {
1270 return Err("Trajectory cannot reach target distance - angle converged without valid solution".into());
1271 }
1272 }
1273 }
1274 }
1275
1276 Err("Failed to find zero angle".into())
1277 }
1278
1279 fn zero_trial_height_at(
1282 &self,
1283 angle_rad: f64,
1284 target_distance_m: f64,
1285 frame: ZeroTargetFrame,
1286 ) -> Result<Option<f64>, BallisticsError> {
1287 let mut trial = self.clone();
1288 trial.inputs.muzzle_angle = angle_rad;
1289 trial.inputs.enable_aerodynamic_jump = false;
1292 trial.inputs.cant_angle = 0.0;
1295 if frame == ZeroTargetFrame::SightLine {
1302 trial.inputs.shooting_angle = 0.0;
1303 }
1304 trial.set_max_range(target_distance_m * 2.0);
1305 let result = trial.solve()?;
1306
1307 for (index, point) in result.points.iter().enumerate() {
1308 if point.position.x >= target_distance_m {
1309 let shot_y_m = if index == 0 {
1310 point.position.y
1311 } else {
1312 let previous = &result.points[index - 1];
1313 let span = point.position.x - previous.position.x;
1314 let fraction = (target_distance_m - previous.position.x) / span;
1315 previous.position.y + fraction * (point.position.y - previous.position.y)
1316 };
1317 return Ok(Some(crate::atmosphere::shot_frame_altitude(
1318 0.0,
1319 target_distance_m,
1320 shot_y_m,
1321 trial.inputs.shooting_angle,
1322 )));
1323 }
1324 }
1325 Ok(None)
1326 }
1327
1328 fn find_zero_range(
1355 &self,
1356 angle_rad: f64,
1357 target_height_m: f64,
1358 frame: ZeroTargetFrame,
1359 ) -> Result<ZeroCrossings, BallisticsError> {
1360 let mut trial = self.clone();
1361 trial.inputs.muzzle_angle = angle_rad;
1362 trial.inputs.enable_aerodynamic_jump = false;
1365 trial.inputs.cant_angle = 0.0;
1366 if frame == ZeroTargetFrame::SightLine {
1367 trial.inputs.shooting_angle = 0.0;
1368 }
1369 let result = trial.solve()?;
1370
1371 let mut near_crossing: Option<f64> = None;
1379 let mut far_crossing: Option<f64> = None;
1380 let mut previous: Option<(f64, f64)> = None; for point in &result.points {
1382 let height = crate::atmosphere::shot_frame_altitude(
1383 0.0,
1384 point.position.x,
1385 point.position.y,
1386 trial.inputs.shooting_angle,
1387 );
1388 let error = height - target_height_m;
1389 if let Some((prev_x, prev_error)) = previous {
1390 if prev_error == 0.0 {
1391 if near_crossing.is_none() {
1394 near_crossing = Some(prev_x);
1395 } else {
1396 far_crossing = Some(prev_x);
1397 }
1398 }
1399 if prev_error * error < 0.0 {
1400 let fraction = prev_error / (prev_error - error);
1401 let crossing = prev_x + fraction * (point.position.x - prev_x);
1402 if prev_error < 0.0 && error > 0.0 {
1403 if near_crossing.is_none() {
1405 near_crossing = Some(crossing);
1406 }
1407 } else {
1408 far_crossing = Some(crossing);
1410 }
1411 }
1412 }
1413 previous = Some((point.position.x, error));
1414 }
1415 if let Some((last_x, last_error)) = previous {
1418 if last_error == 0.0 {
1419 if near_crossing.is_none() {
1420 near_crossing = Some(last_x);
1421 } else {
1422 far_crossing = Some(last_x);
1423 }
1424 }
1425 }
1426
1427 if near_crossing.is_none() && far_crossing.is_none() {
1428 return Err(BallisticsError::from(
1429 "Cannot find zero range: this angle never crosses the target height within the \
1430 solved range (angle too shallow to reach it, or both crossings lie beyond \
1431 the solver's max range)."
1432 .to_string(),
1433 ));
1434 }
1435
1436 Ok(ZeroCrossings {
1437 near_m: near_crossing,
1438 far_m: far_crossing,
1439 })
1440 }
1441
1442 pub fn equivalent_horizontal_range(
1466 &self,
1467 target_range_m: f64,
1468 zero_distance_m: f64,
1469 ) -> Option<f64> {
1470 if !target_range_m.is_finite() || !zero_distance_m.is_finite() {
1471 return None;
1472 }
1473 if target_range_m <= zero_distance_m || target_range_m <= 0.0 {
1474 return None;
1475 }
1476
1477 fn path_y_at(points: &[TrajectoryPoint], distance_m: f64) -> Option<f64> {
1482 match bracket_param(points.len(), |i| points[i].position.x, distance_m) {
1483 Bracket::Below => Some(points[0].position.y),
1488 Bracket::Above | Bracket::Degenerate => None,
1489 Bracket::Inside { lo, t } => {
1490 let hi = lo + 1;
1491 Some(
1492 points[lo].position.y
1493 + t * (points[hi].position.y - points[lo].position.y),
1494 )
1495 }
1496 }
1497 }
1498
1499 let mut inclined = self.clone();
1503 inclined.inputs.enable_trajectory_sampling = false;
1504 let inclined_result = inclined.solve().ok()?;
1505 let los_height = inclined_result.line_of_sight_height_m;
1506 let inclined_drop = los_height - path_y_at(&inclined_result.points, target_range_m)?;
1507 let correction = inclined_drop / target_range_m;
1508 if correction <= 0.0 {
1509 return None;
1510 }
1511
1512 let mut flat = self.clone();
1514 flat.inputs.enable_trajectory_sampling = false;
1515 flat.inputs.shooting_angle = 0.0;
1516 let flat_result = flat.solve().ok()?;
1517 let flat_correction_at = |range_m: f64| -> Option<f64> {
1518 Some((los_height - path_y_at(&flat_result.points, range_m)?) / range_m)
1519 };
1520
1521 let flat_terminal_m = flat_result.points.last().map(|p| p.position.x)?;
1529 let mut low = zero_distance_m.max(1.0);
1530 let mut high = target_range_m.min(flat_terminal_m);
1531 if high <= low {
1532 return None;
1533 }
1534 if flat_correction_at(low)? - correction > 0.0 {
1535 return None; }
1537 if flat_correction_at(high)? - correction < 0.0 {
1538 return None; }
1540 for _ in 0..60 {
1541 let mid = 0.5 * (low + high);
1542 let error = flat_correction_at(mid)? - correction;
1543 if error.abs() == 0.0 {
1544 return Some(mid);
1545 }
1546 if error < 0.0 {
1547 low = mid;
1548 } else {
1549 high = mid;
1550 }
1551 if high - low < 0.01 {
1552 break;
1553 }
1554 }
1555 Some(0.5 * (low + high))
1556 }
1557
1558 fn validate_for_solve(&self) -> Result<(), BallisticsError> {
1564 let require_finite = |name: &str, value: f64| {
1565 if value.is_finite() {
1566 Ok(())
1567 } else {
1568 Err(BallisticsError::from(format!("{name} must be finite")))
1569 }
1570 };
1571 let require_positive = |name: &str, value: f64| {
1572 if value.is_finite() && value > 0.0 {
1573 Ok(())
1574 } else {
1575 Err(BallisticsError::from(format!(
1576 "{name} must be finite and greater than zero"
1577 )))
1578 }
1579 };
1580
1581 if self.inputs.custom_drag_table.is_none() {
1587 require_positive("bc_value", self.inputs.bc_value)?;
1588 }
1589 require_positive("bullet_mass", self.inputs.bullet_mass)?;
1590 require_positive("bullet_diameter", self.inputs.bullet_diameter)?;
1591 require_positive("muzzle_velocity", self.inputs.muzzle_velocity)?;
1592 require_positive("cd_scale", self.inputs.cd_scale)?;
1598
1599 require_finite("muzzle_angle", self.inputs.muzzle_angle)?;
1600 require_finite("azimuth_angle", self.inputs.azimuth_angle)?;
1601 require_finite("shooting_angle", self.inputs.shooting_angle)?;
1602 require_finite("cant_angle", self.inputs.cant_angle)?;
1603 require_finite("muzzle_height", self.inputs.muzzle_height)?;
1604
1605 for (name, value) in [
1609 ("zero_poi_vertical_m", self.inputs.zero_poi_vertical_m),
1610 ("zero_poi_horizontal_m", self.inputs.zero_poi_horizontal_m),
1611 ] {
1612 require_finite(name, value)?;
1613 if value.abs() >= 1.0 {
1614 return Err(BallisticsError::from(format!(
1615 "{name} must be smaller than 1.0 m in magnitude (it is a linear POI \
1616 offset at the zero range, in meters)"
1617 )));
1618 }
1619 }
1620
1621 require_finite(
1625 "sight_offset_lateral_m",
1626 self.inputs.sight_offset_lateral_m,
1627 )?;
1628 if self.inputs.sight_offset_lateral_m.abs() >= 0.5 {
1629 return Err(BallisticsError::from(
1630 "sight_offset_lateral_m must be smaller than 0.5 m in magnitude (it is \
1631 the lateral sight-to-bore mount offset, in meters)",
1632 ));
1633 }
1634
1635 if !(self.inputs.ground_threshold.is_finite()
1638 || self.inputs.ground_threshold == f64::NEG_INFINITY)
1639 {
1640 return Err(BallisticsError::from(
1641 "ground_threshold must be finite or negative infinity",
1642 ));
1643 }
1644
1645 match &self.wind_sock {
1646 Some(wind_sock) => wind_sock
1647 .validate_segments()
1648 .map_err(BallisticsError::from)?,
1649 None => {
1650 require_finite("wind.speed", self.wind.speed)?;
1651 require_finite("wind.direction", self.wind.direction)?;
1652 require_finite("wind.vertical_speed", self.wind.vertical_speed)?;
1653 }
1654 }
1655
1656 require_finite("atmosphere.temperature", self.atmosphere.temperature)?;
1657 require_finite("atmosphere.pressure", self.atmosphere.pressure)?;
1658 require_finite("atmosphere.humidity", self.atmosphere.humidity)?;
1659 require_finite("atmosphere.altitude", self.atmosphere.altitude)?;
1660
1661 require_positive("max_range", self.max_range)?;
1662 if !self.inputs.use_rk4 || !self.inputs.use_adaptive_rk45 {
1665 require_positive("time_step", self.time_step)?;
1666 }
1667
1668 if self.inputs.enable_trajectory_sampling {
1669 require_finite("sight_height", self.inputs.sight_height)?;
1670 require_positive("sample_interval", self.inputs.sample_interval)?;
1671 projected_sample_count(self.max_range, self.inputs.sample_interval)?;
1672 }
1673
1674 if self.inputs.drops_reference == DropsReference::Target {
1678 require_finite("target_height", self.inputs.target_height)?;
1679 if self.inputs.shooting_angle.cos() <= 1e-9 {
1680 return Err(BallisticsError::from(
1681 "drops reference 'target' is undefined for shooting angles at or beyond 90 degrees",
1682 ));
1683 }
1684 }
1685
1686 if self.inputs.enable_coriolis {
1687 require_finite("shot_azimuth", self.inputs.shot_azimuth)?;
1688 if let Some(latitude) = self.inputs.latitude {
1689 require_finite("latitude", latitude)?;
1690 }
1691 }
1692
1693 Ok(())
1694 }
1695
1696 fn validate_result_sanity(&self, result: &TrajectoryResult) -> Result<(), BallisticsError> {
1703 let require_finite = |name: &str, value: f64| {
1704 if value.is_finite() {
1705 Ok(())
1706 } else {
1707 Err(BallisticsError::from(format!(
1708 "trajectory result contains non-finite {name}"
1709 )))
1710 }
1711 };
1712 let require_non_negative = |name: &str, value: f64| {
1713 if value >= 0.0 {
1714 Ok(())
1715 } else {
1716 Err(BallisticsError::from(format!(
1717 "trajectory result contains non-physical negative {name} ({value})"
1718 )))
1719 }
1720 };
1721 let require_indexed_finite = |collection: &str, index: usize, field: &str, value: f64| {
1722 if value.is_finite() {
1723 Ok(())
1724 } else {
1725 Err(BallisticsError::from(format!(
1726 "trajectory result contains non-finite {collection}[{index}].{field}"
1727 )))
1728 }
1729 };
1730 let require_indexed_non_negative =
1731 |collection: &str, index: usize, field: &str, value: f64| {
1732 if value >= 0.0 {
1733 Ok(())
1734 } else {
1735 Err(BallisticsError::from(format!(
1736 "trajectory result contains non-physical negative {collection}[{index}].{field} ({value})"
1737 )))
1738 }
1739 };
1740
1741 require_finite("max_range", result.max_range)?;
1742 require_finite("max_height", result.max_height)?;
1743 require_finite("time_of_flight", result.time_of_flight)?;
1744 require_finite("impact_velocity", result.impact_velocity)?;
1745 require_finite("impact_energy", result.impact_energy)?;
1746 require_finite("projectile_mass_kg", result.projectile_mass_kg)?;
1747 require_finite(
1748 "line_of_sight_height_m",
1749 result.line_of_sight_height_m,
1750 )?;
1751 require_finite(
1752 "station_speed_of_sound_mps",
1753 result.station_speed_of_sound_mps,
1754 )?;
1755
1756 require_non_negative("max_range", result.max_range)?;
1760 require_non_negative("time_of_flight", result.time_of_flight)?;
1761 require_non_negative("impact_velocity", result.impact_velocity)?;
1762 require_non_negative("impact_energy", result.impact_energy)?;
1763 require_non_negative("projectile_mass_kg", result.projectile_mass_kg)?;
1764 require_non_negative(
1765 "station_speed_of_sound_mps",
1766 result.station_speed_of_sound_mps,
1767 )?;
1768
1769 for (index, point) in result.points.iter().enumerate() {
1770 require_indexed_finite("points", index, "time", point.time)?;
1771 require_indexed_finite("points", index, "position.x", point.position.x)?;
1772 require_indexed_finite("points", index, "position.y", point.position.y)?;
1773 require_indexed_finite("points", index, "position.z", point.position.z)?;
1774 require_indexed_finite(
1775 "points",
1776 index,
1777 "velocity_magnitude",
1778 point.velocity_magnitude,
1779 )?;
1780 require_indexed_finite("points", index, "kinetic_energy", point.kinetic_energy)?;
1781 require_indexed_non_negative("points", index, "time", point.time)?;
1782 require_indexed_non_negative(
1783 "points",
1784 index,
1785 "velocity_magnitude",
1786 point.velocity_magnitude,
1787 )?;
1788 require_indexed_non_negative("points", index, "kinetic_energy", point.kinetic_energy)?;
1789 }
1790
1791 if let Some(samples) = &result.sampled_points {
1792 for (index, sample) in samples.iter().enumerate() {
1793 require_indexed_finite("sampled_points", index, "distance_m", sample.distance_m)?;
1794 require_indexed_finite("sampled_points", index, "drop_m", sample.drop_m)?;
1795 require_indexed_finite(
1796 "sampled_points",
1797 index,
1798 "wind_drift_m",
1799 sample.wind_drift_m,
1800 )?;
1801 require_indexed_finite(
1802 "sampled_points",
1803 index,
1804 "velocity_mps",
1805 sample.velocity_mps,
1806 )?;
1807 require_indexed_finite("sampled_points", index, "energy_j", sample.energy_j)?;
1808 require_indexed_finite("sampled_points", index, "time_s", sample.time_s)?;
1809 }
1810 }
1811
1812 for (name, value) in [
1813 ("min_pitch_damping", result.min_pitch_damping),
1814 ("transonic_mach", result.transonic_mach),
1815 ("max_yaw_angle", result.max_yaw_angle),
1816 ("max_precession_angle", result.max_precession_angle),
1817 ] {
1818 if let Some(value) = value {
1819 require_finite(name, value)?;
1820 }
1821 }
1822
1823 if let Some(state) = result.angular_state {
1824 for (name, value) in [
1825 ("angular_state.pitch_angle", state.pitch_angle),
1826 ("angular_state.yaw_angle", state.yaw_angle),
1827 ("angular_state.pitch_rate", state.pitch_rate),
1828 ("angular_state.yaw_rate", state.yaw_rate),
1829 ("angular_state.precession_angle", state.precession_angle),
1830 ("angular_state.nutation_phase", state.nutation_phase),
1831 ] {
1832 require_finite(name, value)?;
1833 }
1834 }
1835
1836 if let Some(jump) = result.aerodynamic_jump {
1837 for (name, value) in [
1838 ("aerodynamic_jump.vertical_jump_moa", jump.vertical_jump_moa),
1839 (
1840 "aerodynamic_jump.horizontal_jump_moa",
1841 jump.horizontal_jump_moa,
1842 ),
1843 ("aerodynamic_jump.jump_angle_rad", jump.jump_angle_rad),
1844 (
1845 "aerodynamic_jump.magnus_component_moa",
1846 jump.magnus_component_moa,
1847 ),
1848 ("aerodynamic_jump.yaw_component_moa", jump.yaw_component_moa),
1849 (
1850 "aerodynamic_jump.stabilization_factor",
1851 jump.stabilization_factor,
1852 ),
1853 ] {
1854 require_finite(name, value)?;
1855 }
1856 }
1857
1858 Ok(())
1859 }
1860
1861 fn validate_integration_state(
1874 &self,
1875 position: &Vector3<f64>,
1876 velocity: &Vector3<f64>,
1877 time: f64,
1878 ) -> Result<(), BallisticsError> {
1879 if !(position.iter().all(|value| value.is_finite())
1880 && velocity.iter().all(|value| value.is_finite())
1881 && time.is_finite())
1882 {
1883 return Err(BallisticsError::from(
1884 "trajectory integration produced a non-finite state (often from physically \
1885 extreme inputs — e.g. an absurd bore/muzzle height placing the launch far \
1886 from sea level, or a degenerate atmosphere; check those inputs, or set \
1887 --altitude explicitly)",
1888 ));
1889 }
1890
1891 let speed = velocity.magnitude();
1892 let budget = self.speed_budget(time);
1893 if speed > budget {
1894 return Err(BallisticsError::from(format!(
1895 "trajectory integration diverged: speed {speed:.3e} m/s at t={time:.6}s exceeds \
1896 the physical budget of {budget:.3e} m/s"
1897 )));
1898 }
1899 Ok(())
1900 }
1901
1902 fn speed_budget(&self, time: f64) -> f64 {
1907 let scalar_wind = self.wind.speed.abs() + self.wind.vertical_speed.abs();
1908 let wind_bound = match &self.wind_sock {
1909 Some(sock) => scalar_wind.max(sock.max_speed_mps()),
1910 None => scalar_wind,
1911 };
1912 2.0 * (self.inputs.muzzle_velocity + wind_bound + 10.0)
1913 + crate::constants::G_ACCEL_MPS2 * time
1914 }
1915
1916 fn push_trajectory_point(
1918 &self,
1919 points: &mut Vec<TrajectoryPoint>,
1920 point: TrajectoryPoint,
1921 ) -> Result<(), BallisticsError> {
1922 if points.len() >= self.max_trajectory_points {
1923 return Err(BallisticsError::from(format!(
1924 "trajectory point limit of {} exceeded",
1925 self.max_trajectory_points
1926 )));
1927 }
1928 points.push(point);
1929 Ok(())
1930 }
1931
1932 pub fn set_wind_segments(&mut self, segments: Vec<crate::wind::WindSegment>) {
1939 self.wind_sock = if segments.is_empty() {
1940 None
1941 } else {
1942 Some(crate::wind::WindSock::new(segments))
1943 };
1944 }
1945
1946 pub fn set_atmo_segments(&mut self, segments: Vec<crate::atmosphere::AtmoSegment>) {
1955 self.atmo_sock = if segments.is_empty() {
1956 None
1957 } else {
1958 Some(crate::atmosphere::AtmoSock::new(segments))
1959 };
1960 }
1961
1962 fn launch_angles_from(
1971 &self,
1972 aj: Option<&crate::aerodynamic_jump::AerodynamicJumpComponents>,
1973 ) -> (f64, f64) {
1974 let (mut elev, mut azim) = (self.inputs.muzzle_angle, self.inputs.azimuth_angle);
1975 if self.inputs.cant_angle != 0.0 {
1982 let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
1983 let (e0, a0) = (elev, azim);
1984 elev = e0 * cos_c - a0 * sin_c;
1985 azim = a0 * cos_c + e0 * sin_c;
1986 }
1987 match aj {
1988 Some(c) => {
1989 const MOA_PER_RAD: f64 = 3437.7467707849;
1991 (
1992 elev + c.vertical_jump_moa / MOA_PER_RAD,
1993 azim + c.horizontal_jump_moa / MOA_PER_RAD,
1994 )
1995 }
1996 None => (elev, azim),
1997 }
1998 }
1999
2000 fn aerodynamic_jump_components(
2008 &self,
2009 ) -> Option<crate::aerodynamic_jump::AerodynamicJumpComponents> {
2010 if !self.inputs.enable_aerodynamic_jump {
2011 return None;
2012 }
2013 let diameter_m = self.inputs.bullet_diameter;
2017 if !(self.inputs.twist_rate.is_finite()
2018 && self.inputs.twist_rate != 0.0
2019 && diameter_m.is_finite()
2020 && diameter_m > 0.0
2021 && self.inputs.bullet_length.is_finite()
2022 && self.inputs.bullet_length > 0.0
2023 && self.inputs.muzzle_velocity.is_finite())
2024 {
2025 return None;
2026 }
2027
2028 let (_, _, temp_c, pressure_hpa) = self.resolved_atmosphere();
2030 let sg = crate::stability::compute_stability_coefficient(
2031 &self.inputs,
2032 (self.atmosphere.altitude, temp_c, pressure_hpa, 0.0),
2033 );
2034 if !(sg.is_finite() && sg > 0.0) {
2035 return None;
2036 }
2037 let length_calibers = self.inputs.bullet_length / diameter_m;
2038
2039 const MS_TO_MPH: f64 = 2.236_936_292_054_4;
2045 let crosswind_from_right_mps = if let Some(sock) = &self.wind_sock {
2046 -sock.vector_for_range_stateless(0.0)[2]
2047 } else {
2048 self.wind.speed * self.wind.direction.sin()
2049 };
2050 let crosswind_from_right_mph = crosswind_from_right_mps * MS_TO_MPH;
2051
2052 let vertical_jump_moa = crate::aerodynamic_jump::litz_crosswind_jump_moa(
2053 sg,
2054 length_calibers,
2055 crosswind_from_right_mph,
2056 self.inputs.is_twist_right,
2057 );
2058 if !vertical_jump_moa.is_finite() {
2059 return None;
2060 }
2061
2062 const MOA_PER_RAD: f64 = 3437.7467707849;
2063 Some(crate::aerodynamic_jump::AerodynamicJumpComponents {
2064 vertical_jump_moa,
2065 horizontal_jump_moa: 0.0,
2067 jump_angle_rad: vertical_jump_moa.abs() / MOA_PER_RAD,
2068 magnus_component_moa: 0.0,
2069 yaw_component_moa: 0.0,
2070 stabilization_factor: (sg / 1.5).clamp(0.0, 1.0),
2071 })
2072 }
2073
2074 fn resolved_atmosphere(&self) -> (f64, f64, f64, f64) {
2075 let (temp_c, pressure_hpa) = match self.station_atmosphere_resolution {
2076 StationAtmosphereResolution::LegacyDefaultSentinels => {
2077 crate::atmosphere::resolve_station_conditions(
2078 self.atmosphere.temperature,
2079 self.atmosphere.pressure,
2080 self.atmosphere.altitude,
2081 )
2082 }
2083 StationAtmosphereResolution::Authoritative => {
2084 (self.atmosphere.temperature, self.atmosphere.pressure)
2085 }
2086 };
2087 let (density, speed_of_sound) = crate::atmosphere::calculate_atmosphere(
2088 self.atmosphere.altitude,
2089 Some(temp_c),
2090 Some(pressure_hpa),
2091 self.atmosphere.humidity,
2092 );
2093 (density, speed_of_sound, temp_c, pressure_hpa)
2094 }
2095
2096 fn precession_nutation_params(
2097 &self,
2098 velocity_mps: f64,
2099 air_density_kg_m3: f64,
2100 speed_of_sound_mps: f64,
2101 ) -> PrecessionNutationParams {
2102 let (spin_inertia, transverse_inertia) = self.precession_nutation_inertias;
2103 let spin_rate_rad_s = if self.inputs.twist_rate > 0.0 {
2104 let velocity_fps = velocity_mps * 3.28084;
2105 let twist_rate_ft = self.inputs.twist_rate / 12.0;
2106 (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI
2107 } else {
2108 0.0
2109 };
2110
2111 PrecessionNutationParams {
2112 mass_kg: self.inputs.bullet_mass,
2113 caliber_m: self.inputs.bullet_diameter,
2114 length_m: self.inputs.bullet_length,
2115 spin_rate_rad_s,
2116 spin_inertia,
2117 transverse_inertia,
2118 velocity_mps,
2119 air_density_kg_m3,
2120 mach: velocity_mps / speed_of_sound_mps,
2121 pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
2122 nutation_damping_factor: 0.05,
2123 }
2124 }
2125
2126 fn append_terminal_endpoint(
2133 &self,
2134 points: &mut Vec<TrajectoryPoint>,
2135 post_position: Vector3<f64>,
2136 post_velocity: Vector3<f64>,
2137 post_time: f64,
2138 max_height: &mut f64,
2139 ) -> Result<TrajectoryTermination, BallisticsError> {
2140 let previous = points
2141 .last()
2142 .cloned()
2143 .ok_or_else(|| BallisticsError::from("No trajectory points generated"))?;
2144
2145 let mut crossings = Vec::with_capacity(3);
2146 if previous.position.x < self.max_range && post_position.x >= self.max_range {
2147 let span = post_position.x - previous.position.x;
2148 if span.is_finite() && span > 0.0 {
2149 crossings.push((
2150 (self.max_range - previous.position.x) / span,
2151 TrajectoryTermination::MaxRange,
2152 ));
2153 }
2154 }
2155 if self.inputs.ground_threshold.is_finite()
2156 && previous.position.y > self.inputs.ground_threshold
2157 && post_position.y <= self.inputs.ground_threshold
2158 {
2159 let span = post_position.y - previous.position.y;
2160 if span.is_finite() && span < 0.0 {
2161 crossings.push((
2162 (self.inputs.ground_threshold - previous.position.y) / span,
2163 TrajectoryTermination::GroundThreshold,
2164 ));
2165 }
2166 }
2167 if previous.time < TRAJECTORY_TIME_LIMIT_S && post_time >= TRAJECTORY_TIME_LIMIT_S {
2168 let span = post_time - previous.time;
2169 if span.is_finite() && span > 0.0 {
2170 crossings.push((
2171 (TRAJECTORY_TIME_LIMIT_S - previous.time) / span,
2172 TrajectoryTermination::TimeLimit,
2173 ));
2174 }
2175 }
2176
2177 let (fraction, termination) = crossings
2178 .into_iter()
2179 .filter(|(fraction, _)| fraction.is_finite() && (0.0..=1.0).contains(fraction))
2180 .min_by(|left, right| {
2181 let priority = |termination: TrajectoryTermination| match termination {
2182 TrajectoryTermination::GroundThreshold => 0,
2183 TrajectoryTermination::MaxRange => 1,
2184 TrajectoryTermination::TimeLimit => 2,
2185 TrajectoryTermination::VelocityFloor => 3,
2186 };
2187 left.0
2188 .total_cmp(&right.0)
2189 .then_with(|| priority(left.1).cmp(&priority(right.1)))
2190 })
2191 .ok_or_else(|| {
2192 BallisticsError::from(
2193 "trajectory integration stopped without crossing a supported boundary",
2194 )
2195 })?;
2196
2197 let mut position = previous.position + (post_position - previous.position) * fraction;
2198 match termination {
2199 TrajectoryTermination::MaxRange => position.x = self.max_range,
2200 TrajectoryTermination::GroundThreshold => {
2201 position.y = self.inputs.ground_threshold;
2202 }
2203 TrajectoryTermination::TimeLimit | TrajectoryTermination::VelocityFloor => {}
2204 }
2205 let velocity_magnitude = previous.velocity_magnitude
2206 + (post_velocity.magnitude() - previous.velocity_magnitude) * fraction;
2207 let mut time = previous.time + (post_time - previous.time) * fraction;
2208 if termination == TrajectoryTermination::TimeLimit {
2209 time = TRAJECTORY_TIME_LIMIT_S;
2210 }
2211 let kinetic_energy =
2212 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2213
2214 if position.y > *max_height {
2215 *max_height = position.y;
2216 }
2217 let terminal_point = TrajectoryPoint {
2218 time,
2219 position,
2220 velocity_magnitude,
2221 kinetic_energy,
2222 drag_coefficient: None,
2223 };
2224 if terminal_point.position.x < previous.position.x {
2225 return Err(BallisticsError::from(
2226 "trajectory terminal state reversed downrange before the crossed boundary",
2227 ));
2228 }
2229 if terminal_point.position.x == previous.position.x {
2230 let last = points.last_mut().ok_or_else(|| {
2235 BallisticsError::from("trajectory points disappeared during terminal finalization")
2236 })?;
2237 *last = terminal_point;
2238 } else {
2239 self.push_trajectory_point(points, terminal_point)?;
2240 }
2241 Ok(termination)
2242 }
2243
2244 fn gravity_acceleration(&self) -> Vector3<f64> {
2245 let theta = self.inputs.shooting_angle;
2246 Vector3::new(
2247 -crate::constants::G_ACCEL_MPS2 * theta.sin(),
2248 -crate::constants::G_ACCEL_MPS2 * theta.cos(),
2249 0.0,
2250 )
2251 }
2252
2253 fn get_wind_at_altitude(&self, altitude_m: f64) -> Vector3<f64> {
2254 let model = match self.inputs.wind_shear_model.as_str() {
2269 "logarithmic" => WindShearModel::Logarithmic,
2270 "power_law" | "powerlaw" | "exponential" => WindShearModel::PowerLaw,
2271 "ekman_spiral" | "ekman" => WindShearModel::EkmanSpiral,
2272 "custom_layers" | "custom" => WindShearModel::CustomLayers,
2273 _ => WindShearModel::PowerLaw,
2274 };
2275 let speed_ratio = crate::wind_shear::boundary_layer_speed_ratio(altitude_m, model);
2276
2277 crate::wind::wind_vector(self.wind.speed, self.wind.direction, 0.0) * speed_ratio
2284 + Vector3::new(0.0, self.wind.vertical_speed, 0.0)
2285 }
2286
2287 pub fn solve(&self) -> Result<TrajectoryResult, BallisticsError> {
2288 self.validate_for_solve()?;
2289 let mut result = if self.inputs.use_rk4 {
2290 if self.inputs.use_adaptive_rk45 {
2291 self.solve_rk45()?
2292 } else {
2293 self.solve_rk4()?
2294 }
2295 } else {
2296 self.solve_euler()?
2297 };
2298 self.apply_spin_drift(&mut result);
2299 self.validate_result_sanity(&result)?;
2300 Ok(result)
2301 }
2302
2303 fn apply_spin_drift(&self, result: &mut TrajectoryResult) {
2309 if !self.inputs.use_enhanced_spin_drift {
2310 return;
2311 }
2312 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 {
2316 return;
2317 }
2318
2319 let sg = self.effective_spin_drift_sg();
2326
2327 for p in result.points.iter_mut() {
2328 if p.time <= 0.0 {
2329 continue;
2330 }
2331 p.position.z +=
2333 crate::spin_drift::litz_drift_meters(sg, p.time, self.inputs.is_twist_right);
2334 }
2335
2336 if let Some(samples) = result.sampled_points.as_mut() {
2340 for s in samples.iter_mut() {
2341 if s.time_s <= 0.0 {
2342 continue;
2343 }
2344 s.wind_drift_m +=
2345 crate::spin_drift::litz_drift_meters(sg, s.time_s, self.inputs.is_twist_right);
2346 }
2347 }
2348 }
2349
2350 fn effective_spin_drift_sg(&self) -> f64 {
2355 let (_, _, temp_c, press_hpa) = self.resolved_atmosphere();
2356 crate::spin_drift::effective_sg_from_inputs(&self.inputs, temp_c, press_hpa)
2357 }
2358
2359 fn initial_position(&self) -> Vector3<f64> {
2377 if self.inputs.cant_angle == 0.0 && self.inputs.sight_offset_lateral_m == 0.0 {
2378 return Vector3::new(0.0, self.inputs.muzzle_height, 0.0);
2379 }
2380 let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
2381 let sh = self.inputs.sight_height;
2382 let off = self.inputs.sight_offset_lateral_m;
2383 Vector3::new(
2384 0.0,
2385 self.inputs.muzzle_height + sh * (1.0 - cos_c) + off * sin_c,
2386 -sh * sin_c - off * cos_c,
2387 )
2388 }
2389
2390 fn build_sampled_points(
2397 &self,
2398 points: &[TrajectoryPoint],
2399 max_height: f64,
2400 transonic_distances: Vec<f64>,
2401 mach_transitions: &MachTransitionTracker,
2402 ) -> Result<Option<Vec<TrajectorySample>>, BallisticsError> {
2403 if !self.inputs.enable_trajectory_sampling {
2404 return Ok(None);
2405 }
2406
2407 let last_point = points.last().ok_or("No trajectory points generated")?;
2408 let trajectory_data = TrajectoryData {
2409 times: points.iter().map(|p| p.time).collect(),
2410 positions: points.iter().map(|p| p.position).collect(),
2411 velocities: points
2412 .iter()
2413 .map(|p| {
2414 Vector3::new(0.0, 0.0, p.velocity_magnitude)
2416 })
2417 .collect(),
2418 transonic_distances, mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2420 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2421 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2422 };
2423
2424 let sight_position_m = self.inputs.muzzle_height + self.inputs.sight_height;
2429 let target_reference = self.inputs.drops_reference == DropsReference::Target;
2430 let target_vertical_height_m = if target_reference && self.inputs.target_height != 0.0 {
2437 self.inputs.target_height
2438 } else {
2439 sight_position_m
2440 };
2441 let outputs = TrajectoryOutputs {
2442 target_distance_horiz_m: last_point.position.x, target_vertical_height_m,
2444 time_of_flight_s: last_point.time,
2445 max_ord_dist_horiz_m: max_height,
2446 sight_height_m: sight_position_m,
2447 };
2448
2449 let mut samples = sample_trajectory(
2451 &trajectory_data,
2452 &outputs,
2453 self.inputs.sample_interval,
2454 self.inputs.bullet_mass,
2455 )?;
2456
2457 if target_reference {
2464 let cos_theta = self.inputs.shooting_angle.cos();
2465 for sample in &mut samples {
2466 sample.drop_m /= cos_theta;
2467 }
2468 }
2469 Ok(Some(samples))
2470 }
2471
2472 fn solve_euler(&self) -> Result<TrajectoryResult, BallisticsError> {
2473 let mut time = 0.0;
2475 let mut position = self.initial_position();
2479 let aj_components = self.aerodynamic_jump_components();
2485 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2486 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2487 let mut velocity = Vector3::new(
2488 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2492
2493 let mut points = Vec::new();
2494 let mut max_height = position.y;
2495 let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
2501 let mut mach_transitions = MachTransitionTracker::default();
2502
2503 let mut angular_state = if self.inputs.enable_precession_nutation {
2505 Some(AngularState {
2506 pitch_angle: 0.001, yaw_angle: 0.001,
2508 pitch_rate: 0.0,
2509 yaw_rate: 0.0,
2510 precession_angle: 0.0,
2511 nutation_phase: 0.0,
2512 })
2513 } else {
2514 None
2515 };
2516 let mut max_yaw_angle = 0.0;
2517 let mut max_precession_angle = 0.0;
2518
2519 let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2521 self.resolved_atmosphere();
2522 let base_ratio = air_density / 1.225;
2527
2528 let wind_vector =
2534 crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
2535
2536 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2539 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2540 );
2541
2542 while position.x < self.max_range
2544 && position.y > self.inputs.ground_threshold
2545 && time < TRAJECTORY_TIME_LIMIT_S
2546 {
2547 let velocity_magnitude = velocity.magnitude();
2549 let kinetic_energy =
2550 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2551
2552 self.push_trajectory_point(
2553 &mut points,
2554 TrajectoryPoint {
2555 time,
2556 position,
2557 velocity_magnitude,
2558 kinetic_energy,
2559 drag_coefficient: None,
2560 },
2561 )?;
2562
2563 {
2566 let mach_here = if speed_of_sound > 0.0 {
2567 velocity_magnitude / speed_of_sound
2568 } else {
2569 0.0
2570 };
2571 mach_transitions.record_downward_crossings(
2572 mach_here,
2573 position.x,
2574 &mut transonic_distances,
2575 );
2576 }
2577
2578 if position.y > max_height {
2580 max_height = position.y;
2581 }
2582
2583 if self.inputs.enable_pitch_damping {
2585 let mach = velocity_magnitude / speed_of_sound;
2586
2587 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
2589 transonic_mach = Some(mach);
2590 }
2591
2592 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
2594
2595 if pitch_damping < min_pitch_damping {
2597 min_pitch_damping = pitch_damping;
2598 }
2599 }
2600
2601 if self.inputs.enable_precession_nutation {
2603 if let Some(ref mut state) = angular_state {
2604 let velocity_magnitude = velocity.magnitude();
2605 let params = self.precession_nutation_params(
2606 velocity_magnitude,
2607 air_density,
2608 speed_of_sound,
2609 );
2610
2611 *state = calculate_combined_angular_motion(
2613 ¶ms,
2614 state,
2615 time,
2616 self.time_step,
2617 0.001, );
2619
2620 if state.yaw_angle.abs() > max_yaw_angle {
2622 max_yaw_angle = state.yaw_angle.abs();
2623 }
2624 if state.precession_angle.abs() > max_precession_angle {
2625 max_precession_angle = state.precession_angle.abs();
2626 }
2627 }
2628 }
2629
2630 let acceleration = self.calculate_acceleration(
2637 &position,
2638 &velocity,
2639 &wind_vector,
2640 (resolved_temp_c, resolved_press_hpa, base_ratio),
2641 );
2642
2643 velocity += acceleration * self.time_step;
2645 position += velocity * self.time_step;
2646 time += self.time_step;
2647 self.validate_integration_state(&position, &velocity, time)?;
2648 }
2649
2650 let termination =
2651 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
2652
2653 self.annotate_drag_coefficients(&mut points, speed_of_sound);
2658
2659 let last_point = points.last().ok_or("No trajectory points generated")?;
2660
2661 let sampled_points = self.build_sampled_points(
2663 &points,
2664 max_height,
2665 transonic_distances,
2666 &mach_transitions,
2667 )?;
2668
2669 Ok(TrajectoryResult {
2670 max_range: last_point.position.x, max_height,
2672 time_of_flight: last_point.time,
2673 impact_velocity: last_point.velocity_magnitude,
2674 impact_energy: last_point.kinetic_energy,
2675 projectile_mass_kg: self.inputs.bullet_mass,
2676 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
2677 station_speed_of_sound_mps: speed_of_sound,
2678 termination,
2679 points,
2680 sampled_points,
2681 min_pitch_damping: if self.inputs.enable_pitch_damping {
2682 Some(min_pitch_damping)
2683 } else {
2684 None
2685 },
2686 transonic_mach,
2687 angular_state,
2688 max_yaw_angle: if self.inputs.enable_precession_nutation {
2689 Some(max_yaw_angle)
2690 } else {
2691 None
2692 },
2693 max_precession_angle: if self.inputs.enable_precession_nutation {
2694 Some(max_precession_angle)
2695 } else {
2696 None
2697 },
2698 aerodynamic_jump: aj_components,
2699 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2700 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2701 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2702 })
2703 }
2704
2705 fn solve_rk4(&self) -> Result<TrajectoryResult, BallisticsError> {
2706 let mut time = 0.0;
2708 let mut position = self.initial_position();
2713
2714 let aj_components = self.aerodynamic_jump_components();
2720 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2721 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2722 let mut velocity = Vector3::new(
2723 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2727
2728 let mut points = Vec::new();
2729 let mut max_height = position.y;
2730 let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
2736 let mut mach_transitions = MachTransitionTracker::default();
2737
2738 let mut angular_state = if self.inputs.enable_precession_nutation {
2740 Some(AngularState {
2741 pitch_angle: 0.001, yaw_angle: 0.001,
2743 pitch_rate: 0.0,
2744 yaw_rate: 0.0,
2745 precession_angle: 0.0,
2746 nutation_phase: 0.0,
2747 })
2748 } else {
2749 None
2750 };
2751 let mut max_yaw_angle = 0.0;
2752 let mut max_precession_angle = 0.0;
2753
2754 let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2756 self.resolved_atmosphere();
2757 let base_ratio = air_density / 1.225;
2762
2763 let wind_vector =
2769 crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
2770
2771 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2774 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2775 );
2776
2777 while position.x < self.max_range
2779 && position.y > self.inputs.ground_threshold
2780 && time < TRAJECTORY_TIME_LIMIT_S
2781 {
2782 let velocity_magnitude = velocity.magnitude();
2784 let kinetic_energy =
2785 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2786
2787 self.push_trajectory_point(
2788 &mut points,
2789 TrajectoryPoint {
2790 time,
2791 position,
2792 velocity_magnitude,
2793 kinetic_energy,
2794 drag_coefficient: None,
2795 },
2796 )?;
2797
2798 {
2801 let mach_here = if speed_of_sound > 0.0 {
2802 velocity_magnitude / speed_of_sound
2803 } else {
2804 0.0
2805 };
2806 mach_transitions.record_downward_crossings(
2807 mach_here,
2808 position.x,
2809 &mut transonic_distances,
2810 );
2811 }
2812
2813 if position.y > max_height {
2814 max_height = position.y;
2815 }
2816
2817 if self.inputs.enable_pitch_damping {
2819 let mach = velocity_magnitude / speed_of_sound;
2820
2821 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
2823 transonic_mach = Some(mach);
2824 }
2825
2826 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
2828
2829 if pitch_damping < min_pitch_damping {
2831 min_pitch_damping = pitch_damping;
2832 }
2833 }
2834
2835 if self.inputs.enable_precession_nutation {
2837 if let Some(ref mut state) = angular_state {
2838 let velocity_magnitude = velocity.magnitude();
2839 let params = self.precession_nutation_params(
2840 velocity_magnitude,
2841 air_density,
2842 speed_of_sound,
2843 );
2844
2845 *state = calculate_combined_angular_motion(
2847 ¶ms,
2848 state,
2849 time,
2850 self.time_step,
2851 0.001, );
2853
2854 if state.yaw_angle.abs() > max_yaw_angle {
2856 max_yaw_angle = state.yaw_angle.abs();
2857 }
2858 if state.precession_angle.abs() > max_precession_angle {
2859 max_precession_angle = state.precession_angle.abs();
2860 }
2861 }
2862 }
2863
2864 let dt = self.time_step;
2866
2867 let acc1 = self.calculate_acceleration(
2869 &position,
2870 &velocity,
2871 &wind_vector,
2872 (resolved_temp_c, resolved_press_hpa, base_ratio),
2873 );
2874
2875 let pos2 = position + velocity * (dt * 0.5);
2877 let vel2 = velocity + acc1 * (dt * 0.5);
2878 let acc2 = self.calculate_acceleration(
2879 &pos2,
2880 &vel2,
2881 &wind_vector,
2882 (resolved_temp_c, resolved_press_hpa, base_ratio),
2883 );
2884
2885 let pos3 = position + vel2 * (dt * 0.5);
2887 let vel3 = velocity + acc2 * (dt * 0.5);
2888 let acc3 = self.calculate_acceleration(
2889 &pos3,
2890 &vel3,
2891 &wind_vector,
2892 (resolved_temp_c, resolved_press_hpa, base_ratio),
2893 );
2894
2895 let pos4 = position + vel3 * dt;
2897 let vel4 = velocity + acc3 * dt;
2898 let acc4 = self.calculate_acceleration(
2899 &pos4,
2900 &vel4,
2901 &wind_vector,
2902 (resolved_temp_c, resolved_press_hpa, base_ratio),
2903 );
2904
2905 position += (velocity + vel2 * 2.0 + vel3 * 2.0 + vel4) * (dt / 6.0);
2907 velocity += (acc1 + acc2 * 2.0 + acc3 * 2.0 + acc4) * (dt / 6.0);
2908 time += dt;
2909 self.validate_integration_state(&position, &velocity, time)?;
2910 }
2911
2912 let termination =
2913 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
2914
2915 self.annotate_drag_coefficients(&mut points, speed_of_sound);
2920
2921 let last_point = points.last().ok_or("No trajectory points generated")?;
2922
2923 let sampled_points = self.build_sampled_points(
2925 &points,
2926 max_height,
2927 transonic_distances,
2928 &mach_transitions,
2929 )?;
2930
2931 Ok(TrajectoryResult {
2932 max_range: last_point.position.x, max_height,
2934 time_of_flight: last_point.time,
2935 impact_velocity: last_point.velocity_magnitude,
2936 impact_energy: last_point.kinetic_energy,
2937 projectile_mass_kg: self.inputs.bullet_mass,
2938 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
2939 station_speed_of_sound_mps: speed_of_sound,
2940 termination,
2941 points,
2942 sampled_points,
2943 min_pitch_damping: if self.inputs.enable_pitch_damping {
2944 Some(min_pitch_damping)
2945 } else {
2946 None
2947 },
2948 transonic_mach,
2949 angular_state,
2950 max_yaw_angle: if self.inputs.enable_precession_nutation {
2951 Some(max_yaw_angle)
2952 } else {
2953 None
2954 },
2955 max_precession_angle: if self.inputs.enable_precession_nutation {
2956 Some(max_precession_angle)
2957 } else {
2958 None
2959 },
2960 aerodynamic_jump: aj_components,
2961 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2962 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2963 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2964 })
2965 }
2966
2967 fn solve_rk45(&self) -> Result<TrajectoryResult, BallisticsError> {
2968 let mut time = 0.0;
2970 let mut position = self.initial_position();
2974
2975 let aj_components = self.aerodynamic_jump_components();
2981 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2982 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2983 let mut velocity = Vector3::new(
2984 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2988
2989 let mut points = Vec::new();
2990 let mut max_height = position.y;
2991 let mut dt = 0.001; let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2996 self.resolved_atmosphere();
2997 let base_ratio = air_density / 1.225;
3002 let wind_vector =
3007 crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
3008
3009 let mut transonic_distances: Vec<f64> = Vec::new();
3011 let mut mach_transitions = MachTransitionTracker::default();
3012
3013 let mut min_pitch_damping = f64::INFINITY;
3018 let mut transonic_mach: Option<f64> = None;
3019 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
3020 self.inputs.bullet_model.as_deref().unwrap_or("default"),
3021 );
3022 let mut angular_state = if self.inputs.enable_precession_nutation {
3023 Some(AngularState {
3024 pitch_angle: 0.001,
3025 yaw_angle: 0.001,
3026 pitch_rate: 0.0,
3027 yaw_rate: 0.0,
3028 precession_angle: 0.0,
3029 nutation_phase: 0.0,
3030 })
3031 } else {
3032 None
3033 };
3034 let mut max_yaw_angle = 0.0;
3035 let mut max_precession_angle = 0.0;
3036
3037 while position.x < self.max_range
3038 && position.y > self.inputs.ground_threshold
3039 && time < TRAJECTORY_TIME_LIMIT_S
3040 {
3041 let velocity_magnitude = velocity.magnitude();
3043 let kinetic_energy = 0.5 * self.inputs.bullet_mass * velocity_magnitude.powi(2);
3044
3045 self.push_trajectory_point(
3046 &mut points,
3047 TrajectoryPoint {
3048 time,
3049 position,
3050 velocity_magnitude,
3051 kinetic_energy,
3052 drag_coefficient: None,
3053 },
3054 )?;
3055
3056 {
3059 let mach_here = if speed_of_sound > 0.0 {
3060 velocity_magnitude / speed_of_sound
3061 } else {
3062 0.0
3063 };
3064 mach_transitions.record_downward_crossings(
3065 mach_here,
3066 position.x,
3067 &mut transonic_distances,
3068 );
3069 }
3070
3071 if position.y > max_height {
3072 max_height = position.y;
3073 }
3074
3075 if self.inputs.enable_pitch_damping {
3078 let mach = velocity_magnitude / speed_of_sound;
3079 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
3080 transonic_mach = Some(mach);
3081 }
3082 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
3083 if pitch_damping < min_pitch_damping {
3084 min_pitch_damping = pitch_damping;
3085 }
3086 }
3087
3088 let accepted_step = self.adaptive_rk45_step(
3091 &position,
3092 &velocity,
3093 dt,
3094 &wind_vector,
3095 (resolved_temp_c, resolved_press_hpa, base_ratio),
3096 );
3097 debug_assert!(
3098 accepted_step.error <= RK45_TOLERANCE || accepted_step.used_dt <= RK45_MIN_DT
3099 );
3100
3101 if self.inputs.enable_precession_nutation {
3105 if let Some(ref mut state) = angular_state {
3106 let params = self.precession_nutation_params(
3107 velocity_magnitude,
3108 air_density,
3109 speed_of_sound,
3110 );
3111
3112 *state = calculate_combined_angular_motion(
3113 ¶ms,
3114 state,
3115 time,
3116 accepted_step.used_dt,
3117 0.001,
3118 );
3119
3120 if state.yaw_angle.abs() > max_yaw_angle {
3121 max_yaw_angle = state.yaw_angle.abs();
3122 }
3123 if state.precession_angle.abs() > max_precession_angle {
3124 max_precession_angle = state.precession_angle.abs();
3125 }
3126 }
3127 }
3128
3129 position = accepted_step.position;
3130 velocity = accepted_step.velocity;
3131 time += accepted_step.used_dt;
3132 self.validate_integration_state(&position, &velocity, time)?;
3133
3134 dt = accepted_step.next_dt;
3136 }
3137
3138 if points.is_empty() {
3140 return Err(BallisticsError::from("No trajectory points calculated"));
3141 }
3142
3143 let termination =
3145 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
3146
3147 self.annotate_drag_coefficients(&mut points, speed_of_sound);
3149
3150 let last_point = points.last().unwrap();
3151
3152 let sampled_points = self.build_sampled_points(
3154 &points,
3155 max_height,
3156 transonic_distances,
3157 &mach_transitions,
3158 )?;
3159
3160 Ok(TrajectoryResult {
3161 max_range: last_point.position.x, max_height,
3163 time_of_flight: last_point.time,
3164 impact_velocity: last_point.velocity_magnitude,
3165 impact_energy: last_point.kinetic_energy,
3166 projectile_mass_kg: self.inputs.bullet_mass,
3167 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
3168 station_speed_of_sound_mps: speed_of_sound,
3169 termination,
3170 points,
3171 sampled_points,
3172 min_pitch_damping: if self.inputs.enable_pitch_damping {
3173 Some(min_pitch_damping)
3174 } else {
3175 None
3176 },
3177 transonic_mach,
3178 angular_state,
3179 max_yaw_angle: if self.inputs.enable_precession_nutation {
3180 Some(max_yaw_angle)
3181 } else {
3182 None
3183 },
3184 max_precession_angle: if self.inputs.enable_precession_nutation {
3185 Some(max_precession_angle)
3186 } else {
3187 None
3188 },
3189 aerodynamic_jump: aj_components,
3190 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
3191 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
3192 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
3193 })
3194 }
3195
3196 fn adaptive_rk45_step(
3197 &self,
3198 position: &Vector3<f64>,
3199 velocity: &Vector3<f64>,
3200 initial_dt: f64,
3201 wind_vector: &Vector3<f64>,
3202 resolved_atmo: (f64, f64, f64),
3203 ) -> Rk45AcceptedStep {
3204 let mut trial_dt = initial_dt;
3205
3206 loop {
3207 let trial = self.rk45_step(
3208 position,
3209 velocity,
3210 trial_dt,
3211 wind_vector,
3212 RK45_TOLERANCE,
3213 resolved_atmo,
3214 );
3215 let next_dt = if trial.suggested_dt.is_finite() {
3220 (RK45_SAFETY_FACTOR * trial.suggested_dt).clamp(RK45_MIN_DT, RK45_MAX_DT)
3221 } else {
3222 RK45_MIN_DT
3223 };
3224
3225 if trial.error <= RK45_TOLERANCE || trial_dt <= RK45_MIN_DT {
3226 return Rk45AcceptedStep {
3227 position: trial.position,
3228 velocity: trial.velocity,
3229 used_dt: trial_dt,
3230 next_dt,
3231 error: trial.error,
3232 };
3233 }
3234
3235 trial_dt = next_dt;
3236 }
3237 }
3238
3239 fn rk45_step(
3240 &self,
3241 position: &Vector3<f64>,
3242 velocity: &Vector3<f64>,
3243 dt: f64,
3244 wind_vector: &Vector3<f64>,
3245 tolerance: f64,
3246 resolved_atmo: (f64, f64, f64), ) -> Rk45Trial {
3248 const A21: f64 = 1.0 / 5.0;
3250 const A31: f64 = 3.0 / 40.0;
3251 const A32: f64 = 9.0 / 40.0;
3252 const A41: f64 = 44.0 / 45.0;
3253 const A42: f64 = -56.0 / 15.0;
3254 const A43: f64 = 32.0 / 9.0;
3255 const A51: f64 = 19372.0 / 6561.0;
3256 const A52: f64 = -25360.0 / 2187.0;
3257 const A53: f64 = 64448.0 / 6561.0;
3258 const A54: f64 = -212.0 / 729.0;
3259 const A61: f64 = 9017.0 / 3168.0;
3260 const A62: f64 = -355.0 / 33.0;
3261 const A63: f64 = 46732.0 / 5247.0;
3262 const A64: f64 = 49.0 / 176.0;
3263 const A65: f64 = -5103.0 / 18656.0;
3264 const A71: f64 = 35.0 / 384.0;
3265 const A73: f64 = 500.0 / 1113.0;
3266 const A74: f64 = 125.0 / 192.0;
3267 const A75: f64 = -2187.0 / 6784.0;
3268 const A76: f64 = 11.0 / 84.0;
3269
3270 const B1: f64 = 35.0 / 384.0;
3272 const B3: f64 = 500.0 / 1113.0;
3273 const B4: f64 = 125.0 / 192.0;
3274 const B5: f64 = -2187.0 / 6784.0;
3275 const B6: f64 = 11.0 / 84.0;
3276
3277 const B1_ERR: f64 = 5179.0 / 57600.0;
3279 const B3_ERR: f64 = 7571.0 / 16695.0;
3280 const B4_ERR: f64 = 393.0 / 640.0;
3281 const B5_ERR: f64 = -92097.0 / 339200.0;
3282 const B6_ERR: f64 = 187.0 / 2100.0;
3283 const B7_ERR: f64 = 1.0 / 40.0;
3284
3285 let k1_v = self.calculate_acceleration(position, velocity, wind_vector, resolved_atmo);
3287 let k1_p = *velocity;
3288
3289 let p2 = position + dt * A21 * k1_p;
3290 let v2 = velocity + dt * A21 * k1_v;
3291 let k2_v = self.calculate_acceleration(&p2, &v2, wind_vector, resolved_atmo);
3292 let k2_p = v2;
3293
3294 let p3 = position + dt * (A31 * k1_p + A32 * k2_p);
3295 let v3 = velocity + dt * (A31 * k1_v + A32 * k2_v);
3296 let k3_v = self.calculate_acceleration(&p3, &v3, wind_vector, resolved_atmo);
3297 let k3_p = v3;
3298
3299 let p4 = position + dt * (A41 * k1_p + A42 * k2_p + A43 * k3_p);
3300 let v4 = velocity + dt * (A41 * k1_v + A42 * k2_v + A43 * k3_v);
3301 let k4_v = self.calculate_acceleration(&p4, &v4, wind_vector, resolved_atmo);
3302 let k4_p = v4;
3303
3304 let p5 = position + dt * (A51 * k1_p + A52 * k2_p + A53 * k3_p + A54 * k4_p);
3305 let v5 = velocity + dt * (A51 * k1_v + A52 * k2_v + A53 * k3_v + A54 * k4_v);
3306 let k5_v = self.calculate_acceleration(&p5, &v5, wind_vector, resolved_atmo);
3307 let k5_p = v5;
3308
3309 let p6 = position + dt * (A61 * k1_p + A62 * k2_p + A63 * k3_p + A64 * k4_p + A65 * k5_p);
3310 let v6 = velocity + dt * (A61 * k1_v + A62 * k2_v + A63 * k3_v + A64 * k4_v + A65 * k5_v);
3311 let k6_v = self.calculate_acceleration(&p6, &v6, wind_vector, resolved_atmo);
3312 let k6_p = v6;
3313
3314 let p7 = position + dt * (A71 * k1_p + A73 * k3_p + A74 * k4_p + A75 * k5_p + A76 * k6_p);
3315 let v7 = velocity + dt * (A71 * k1_v + A73 * k3_v + A74 * k4_v + A75 * k5_v + A76 * k6_v);
3316 let k7_v = self.calculate_acceleration(&p7, &v7, wind_vector, resolved_atmo);
3317 let k7_p = v7;
3318
3319 let new_pos = position + dt * (B1 * k1_p + B3 * k3_p + B4 * k4_p + B5 * k5_p + B6 * k6_p);
3321 let new_vel = velocity + dt * (B1 * k1_v + B3 * k3_v + B4 * k4_v + B5 * k5_v + B6 * k6_v);
3322
3323 let pos_err = position
3325 + dt * (B1_ERR * k1_p
3326 + B3_ERR * k3_p
3327 + B4_ERR * k4_p
3328 + B5_ERR * k5_p
3329 + B6_ERR * k6_p
3330 + B7_ERR * k7_p);
3331 let vel_err = velocity
3332 + dt * (B1_ERR * k1_v
3333 + B3_ERR * k3_v
3334 + B4_ERR * k4_v
3335 + B5_ERR * k5_v
3336 + B6_ERR * k6_v
3337 + B7_ERR * k7_v);
3338
3339 let error = cli_rk45_error_norm(position, velocity, &new_pos, &new_vel, &pos_err, &vel_err);
3341
3342 let dt_new = if error < tolerance {
3344 dt * (tolerance / error).powf(0.2).min(2.0)
3345 } else {
3346 dt * (tolerance / error).powf(0.25).max(0.1)
3347 };
3348
3349 Rk45Trial {
3350 position: new_pos,
3351 velocity: new_vel,
3352 suggested_dt: dt_new,
3353 error,
3354 }
3355 }
3356
3357 fn apply_cluster_bc_correction(&self, base_bc: f64, velocity_fps: f64) -> f64 {
3358 if let Some(ref cluster_bc) = self.cluster_bc {
3359 cluster_bc.apply_correction_for_drag_model(
3360 base_bc,
3361 self.inputs.caliber_inches,
3362 self.inputs.weight_grains,
3363 velocity_fps,
3364 self.inputs.bc_type,
3365 )
3366 } else {
3367 base_bc
3368 }
3369 }
3370
3371 fn calculate_acceleration(
3372 &self,
3373 position: &Vector3<f64>,
3374 velocity: &Vector3<f64>,
3375 wind_vector: &Vector3<f64>,
3376 resolved_atmo: (f64, f64, f64), ) -> Vector3<f64> {
3378 let actual_wind = if let Some(ref sock) = self.wind_sock {
3384 sock.vector_for_range_stateless(position.x)
3385 } else if self.inputs.enable_wind_shear {
3386 self.get_wind_at_altitude(position.y)
3387 } else {
3388 *wind_vector
3389 };
3390 let actual_wind =
3391 crate::derivatives::level_vector_to_shot_frame(actual_wind, self.inputs.shooting_angle);
3392
3393 let relative_velocity = velocity - actual_wind;
3394 let velocity_magnitude = relative_velocity.magnitude();
3395
3396 if velocity_magnitude < 0.001 {
3397 return self.gravity_acceleration();
3398 }
3399
3400 let (base_temp_c, base_press_hpa, station_ratio) = resolved_atmo;
3411
3412 let (drag_base_temp_c, drag_base_press_hpa, drag_base_ratio, drag_humidity_percent) =
3420 if let Some(ref sock) = self.atmo_sock {
3421 let (zone_temp_c, zone_press_hpa, zone_humidity) = sock.atmo_for_range(position.x);
3422 let zone_base_ratio = crate::atmosphere::calculate_air_density_cimp(
3423 zone_temp_c,
3424 zone_press_hpa,
3425 zone_humidity,
3426 ) / 1.225;
3427 (zone_temp_c, zone_press_hpa, zone_base_ratio, zone_humidity)
3428 } else {
3429 (
3430 base_temp_c,
3431 base_press_hpa,
3432 station_ratio,
3433 self.atmosphere.humidity,
3434 )
3435 };
3436 let local_alt = crate::atmosphere::shot_frame_altitude(
3437 self.atmosphere.altitude,
3438 position.x,
3439 position.y,
3440 self.inputs.shooting_angle,
3441 );
3442 let (air_density, speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
3443 local_alt,
3444 self.atmosphere.altitude,
3445 drag_base_temp_c,
3446 drag_base_press_hpa,
3447 drag_base_ratio,
3448 drag_humidity_percent,
3449 );
3450
3451 let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
3455
3456 let velocity_fps = velocity_magnitude * 3.28084;
3458
3459 let cd_to_retard = crate::constants::CD_TO_RETARD;
3464 let standard_factor = cd * cd_to_retard;
3465 let density_scale = air_density / 1.225; let a_drag_ft_s2 =
3469 (velocity_fps * velocity_fps) * standard_factor * density_scale / retard_denom;
3470 let a_drag_m_s2 = a_drag_ft_s2 * 0.3048; let drag_acceleration = -a_drag_m_s2 * (relative_velocity / velocity_magnitude);
3474
3475 let mut accel = drag_acceleration + self.gravity_acceleration();
3478
3479 if self.inputs.enable_coriolis {
3482 if let Some(lat_deg) = self.inputs.latitude {
3483 let omega_earth = 7.2921159e-5_f64; let lat = lat_deg.to_radians();
3485 let az = self.inputs.shot_azimuth; let omega = Vector3::new(
3492 omega_earth * lat.cos() * az.cos(), omega_earth * lat.sin(), -omega_earth * lat.cos() * az.sin(), );
3496 let omega = crate::derivatives::level_vector_to_shot_frame(
3497 omega,
3498 self.inputs.shooting_angle,
3499 );
3500 accel += -2.0 * omega.cross(velocity);
3505 }
3506 }
3507
3508 if self.inputs.enable_magnus
3515 && !self.inputs.use_enhanced_spin_drift
3516 && self.inputs.bullet_diameter > 0.0
3517 && self.inputs.twist_rate > 0.0
3518 {
3519 let diameter_m = self.inputs.bullet_diameter;
3520 let (spin_rad_s, spin_param) = crate::spin_drift::calculate_magnus_spin_state(
3521 self.inputs.muzzle_velocity,
3522 velocity_magnitude,
3523 self.inputs.twist_rate,
3524 diameter_m,
3525 );
3526 let mach = velocity_magnitude / speed_of_sound;
3528
3529 let d_in = self.inputs.bullet_diameter / 0.0254;
3531 let m_gr = self.inputs.bullet_mass / crate::constants::GRAINS_TO_KG;
3532 let l_in = if self.inputs.bullet_length > 0.0 {
3533 self.inputs.bullet_length / 0.0254
3534 } else {
3535 let est_m = crate::stability::estimate_bullet_length_m(
3537 self.inputs.bullet_diameter,
3538 self.inputs.bullet_mass,
3539 );
3540 if est_m > 0.0 {
3541 est_m / 0.0254
3542 } else {
3543 4.5 * d_in
3544 }
3545 };
3546 let sg = crate::spin_drift::calculate_dynamic_stability(
3550 m_gr,
3551 velocity_magnitude,
3552 spin_rad_s,
3553 d_in,
3554 l_in,
3555 air_density,
3556 );
3557
3558 let (yaw_rad, _) = crate::spin_drift::calculate_yaw_of_repose(
3560 sg,
3561 velocity_magnitude,
3562 spin_rad_s,
3563 0.0, 0.0, air_density,
3566 d_in,
3567 l_in,
3568 m_gr,
3569 mach,
3570 "match",
3571 false,
3572 );
3573
3574 let c_np = crate::derivatives::calculate_magnus_moment_coefficient(mach);
3576 let area = std::f64::consts::PI * (diameter_m / 2.0).powi(2);
3577 let magnus_force = 0.5
3578 * air_density
3579 * velocity_magnitude.powi(2)
3580 * area
3581 * c_np
3582 * spin_param
3583 * yaw_rad.sin();
3584
3585 if magnus_force.abs() > 1e-12 {
3589 if let Some(dir) = crate::derivatives::yaw_of_repose_magnus_direction(
3590 relative_velocity,
3591 self.gravity_acceleration(),
3592 self.inputs.is_twist_right,
3593 ) {
3594 accel += (magnus_force / self.inputs.bullet_mass) * dir;
3595 }
3596 }
3597 }
3598
3599 accel
3600 }
3601
3602 fn drag_terms(&self, velocity_magnitude: f64, speed_of_sound: f64) -> (f64, f64) {
3614 let cd = self.calculate_drag_coefficient(velocity_magnitude, speed_of_sound);
3615
3616 let velocity_fps = velocity_magnitude * 3.28084;
3617
3618 let (base_bc, bc_from_segments) = if let Some(segments) = self
3623 .inputs
3624 .bc_segments_data
3625 .as_ref()
3626 .filter(|segments| self.inputs.use_bc_segments && !segments.is_empty())
3627 {
3628 (
3630 crate::bc_estimation::velocity_segment_bc(
3631 velocity_fps,
3632 segments,
3633 self.inputs.bc_value,
3634 ),
3635 true,
3636 )
3637 } else if let Some(segments) = self
3638 .inputs
3639 .bc_segments
3640 .as_ref()
3641 .filter(|segments| !segments.is_empty())
3642 {
3643 (
3644 crate::derivatives::interpolated_bc(
3645 velocity_magnitude / speed_of_sound,
3646 segments,
3647 Some(&self.inputs),
3648 ),
3649 true,
3650 )
3651 } else {
3652 (self.inputs.bc_value, false)
3653 };
3654
3655 let effective_bc = if bc_from_segments {
3660 base_bc
3661 } else {
3662 self.apply_cluster_bc_correction(base_bc, velocity_fps)
3663 };
3664 let effective_bc = effective_bc.max(1e-6);
3667
3668 let retard_denom = if self.inputs.custom_drag_table.is_some() {
3673 self.inputs.custom_drag_denominator(effective_bc)
3674 } else {
3675 effective_bc
3676 };
3677
3678 (cd, retard_denom)
3679 }
3680
3681 pub fn effective_drag_coefficient(
3702 &self,
3703 velocity_magnitude: f64,
3704 speed_of_sound: f64,
3705 ) -> Option<f64> {
3706 if !velocity_magnitude.is_finite() || speed_of_sound <= 1e-9 {
3707 return None;
3708 }
3709 let sectional_density = self.inputs.sectional_density_lb_in2()?;
3710 let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
3711 if retard_denom <= 0.0 {
3712 return None;
3713 }
3714 let effective = cd * sectional_density / retard_denom;
3715 effective.is_finite().then_some(effective)
3716 }
3717
3718 fn annotate_drag_coefficients(&self, points: &mut [TrajectoryPoint], speed_of_sound: f64) {
3728 for point in points.iter_mut() {
3729 point.drag_coefficient =
3730 self.effective_drag_coefficient(point.velocity_magnitude, speed_of_sound);
3731 }
3732 }
3733
3734 fn calculate_drag_coefficient(&self, velocity: f64, speed_of_sound: f64) -> f64 {
3735 let mach = velocity / speed_of_sound;
3736
3737 if let Some(ref table) = self.inputs.custom_drag_table {
3741 return table.interpolate(mach) * self.inputs.cd_scale;
3746 }
3747
3748 crate::drag::get_drag_coefficient(mach, &self.inputs.bc_type)
3751 }
3752}
3753
3754#[derive(Debug, Clone)]
3756pub struct MonteCarloParams {
3757 pub num_simulations: usize,
3758 pub velocity_std_dev: f64,
3759 pub angle_std_dev: f64,
3760 pub bc_std_dev: f64,
3761 pub wind_speed_std_dev: f64,
3762 pub target_distance: Option<f64>,
3763 pub base_wind_speed: f64,
3764 pub base_wind_direction: f64,
3765 pub azimuth_std_dev: f64, }
3767
3768impl Default for MonteCarloParams {
3769 fn default() -> Self {
3770 Self {
3771 num_simulations: 1000,
3772 velocity_std_dev: 1.0,
3773 angle_std_dev: 0.001,
3774 bc_std_dev: 0.01,
3775 wind_speed_std_dev: 1.0,
3776 target_distance: None,
3777 base_wind_speed: 0.0,
3778 base_wind_direction: 0.0,
3779 azimuth_std_dev: 0.001, }
3781 }
3782}
3783
3784#[derive(Debug, Clone)]
3786pub struct MonteCarloResults {
3787 pub ranges: Vec<f64>,
3788 pub impact_velocities: Vec<f64>,
3789 pub impact_positions: Vec<Vector3<f64>>,
3795}
3796
3797pub const DEFAULT_HIT_RADIUS_M: f64 = 0.3;
3800
3801pub const TARGET_NOT_REACHED_SENTINEL_M: f64 = -1.0e9;
3807
3808impl MonteCarloResults {
3809 pub fn position_reached_target(position: &Vector3<f64>) -> bool {
3811 position.iter().all(|component| component.is_finite())
3812 && position.y != TARGET_NOT_REACHED_SENTINEL_M
3813 }
3814
3815 pub fn target_arrival_count(&self) -> usize {
3817 self.impact_positions
3818 .iter()
3819 .filter(|position| Self::position_reached_target(position))
3820 .count()
3821 }
3822
3823 pub fn target_shortfall_fraction(&self) -> f64 {
3826 if self.impact_positions.is_empty() {
3827 return 0.0;
3828 }
3829 (self.impact_positions.len() - self.target_arrival_count()) as f64
3830 / self.impact_positions.len() as f64
3831 }
3832
3833 pub fn target_plane_cep(&self) -> Option<f64> {
3839 let mut radial_misses: Vec<f64> = self
3840 .impact_positions
3841 .iter()
3842 .filter(|position| Self::position_reached_target(position))
3843 .map(Vector3::norm)
3844 .filter(|miss| miss.is_finite())
3845 .collect();
3846 radial_misses.sort_by(f64::total_cmp);
3847 if radial_misses.is_empty() {
3848 None
3849 } else {
3850 Some(radial_misses[radial_misses.len() / 2])
3851 }
3852 }
3853
3854 pub fn hit_probability(&self, hit_radius_m: f64) -> f64 {
3863 if self.impact_positions.is_empty() {
3864 return 0.0;
3865 }
3866 let hits = self
3867 .impact_positions
3868 .iter()
3869 .filter(|position| Self::position_is_hit(position, hit_radius_m))
3870 .count();
3871 hits as f64 / self.impact_positions.len() as f64
3872 }
3873
3874 pub fn position_is_hit(position: &Vector3<f64>, hit_radius_m: f64) -> bool {
3886 Self::position_reached_target(position) && position.norm() < hit_radius_m
3887 }
3888
3889 pub fn hit_probability_wilson(
3911 &self,
3912 hit_radius_m: f64,
3913 level: ConfidenceLevel,
3914 ) -> (f64, (f64, f64), u64) {
3915 let trials = self.impact_positions.len() as u64;
3916 let hits = self
3917 .impact_positions
3918 .iter()
3919 .filter(|position| Self::position_is_hit(position, hit_radius_m))
3920 .count() as u64;
3921 (
3922 self.hit_probability(hit_radius_m),
3923 wilson_interval(hits, trials, level),
3924 trials,
3925 )
3926 }
3927
3928 pub fn rect_hit_probability(&self, width_m: f64, height_m: f64) -> f64 {
3940 let dimensions_invalid = width_m.is_nan()
3941 || width_m <= 0.0
3942 || height_m.is_nan()
3943 || height_m <= 0.0;
3944 if self.impact_positions.is_empty() || dimensions_invalid {
3945 return 0.0;
3946 }
3947 let half_width = width_m / 2.0;
3948 let half_height = height_m / 2.0;
3949 let hits = self
3950 .impact_positions
3951 .iter()
3952 .filter(|position| {
3953 Self::position_reached_target(position)
3954 && position.z.abs() <= half_width
3955 && position.y.abs() <= half_height
3956 })
3957 .count();
3958 hits as f64 / self.impact_positions.len() as f64
3959 }
3960}
3961
3962fn wind_from_signed_speed_sample(
3963 signed_speed: f64,
3964 sampled_direction: f64,
3965 vertical_speed: f64,
3966) -> WindConditions {
3967 if signed_speed < 0.0 {
3972 WindConditions {
3973 speed: -signed_speed,
3974 direction: sampled_direction + std::f64::consts::PI,
3975 vertical_speed,
3976 }
3977 } else {
3978 WindConditions {
3979 speed: signed_speed,
3980 direction: sampled_direction,
3981 vertical_speed,
3982 }
3983 }
3984}
3985
3986struct MonteCarloWindSampler {
3987 speed: rand_distr::Normal<f64>,
3988 direction: rand_distr::Normal<f64>,
3989 vertical_speed: f64,
3991}
3992
3993impl MonteCarloWindSampler {
3994 fn new(
3995 base_wind: &WindConditions,
3996 wind_speed_std_dev: f64,
3997 wind_direction_std_dev: f64,
3998 ) -> Result<Self, BallisticsError> {
3999 use rand_distr::Normal;
4000
4001 if !wind_direction_std_dev.is_finite() || wind_direction_std_dev < 0.0 {
4002 return Err("Wind direction standard deviation must be finite and non-negative".into());
4003 }
4004
4005 let speed = Normal::new(base_wind.speed, wind_speed_std_dev)
4006 .map_err(|e| format!("Invalid wind speed distribution: {e}"))?;
4007 let direction = Normal::new(base_wind.direction, wind_direction_std_dev)
4008 .map_err(|e| format!("Invalid wind direction distribution: {e}"))?;
4009 Ok(Self { speed, direction, vertical_speed: base_wind.vertical_speed })
4010 }
4011
4012 fn sample<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> WindConditions {
4013 use rand_distr::Distribution;
4014
4015 wind_from_signed_speed_sample(
4016 self.speed.sample(rng),
4017 self.direction.sample(rng),
4018 self.vertical_speed,
4019 )
4020 }
4021}
4022
4023#[derive(Debug, Clone, Copy)]
4026struct TrialOutcome {
4027 range: f64,
4030 impact_velocity: f64,
4032 impact_position: Vector3<f64>,
4035}
4036
4037struct MonteCarloTrialSampler {
4057 base_inputs: BallisticInputs,
4058 atmosphere: AtmosphericConditions,
4059 solver_max_range: f64,
4060 target_distance: f64,
4063 baseline_at_target: Vector3<f64>,
4066 velocity_delta_dist: rand_distr::Normal<f64>,
4067 angle_dist: rand_distr::Normal<f64>,
4068 bc_dist: rand_distr::Normal<f64>,
4069 wind_sampler: MonteCarloWindSampler,
4070 azimuth_dist: rand_distr::Normal<f64>,
4071}
4072
4073impl MonteCarloTrialSampler {
4074 fn new(
4078 base_inputs: BallisticInputs,
4079 base_wind: &WindConditions,
4080 params: &MonteCarloParams,
4081 wind_direction_std_dev: f64,
4082 ) -> Result<Self, BallisticsError> {
4083 use rand_distr::Normal;
4084
4085 let atmosphere = AtmosphericConditions {
4086 temperature: base_inputs.temperature,
4087 pressure: base_inputs.pressure,
4088 humidity: base_inputs.humidity_percent(),
4089 altitude: base_inputs.altitude,
4090 };
4091 let target_hint = params
4092 .target_distance
4093 .unwrap_or(base_inputs.target_distance);
4094 let solver_max_range = target_hint.max(1000.0) * 2.0;
4095
4096 let mut baseline_solver =
4098 TrajectorySolver::new(base_inputs.clone(), base_wind.clone(), atmosphere.clone());
4099 baseline_solver.set_max_range(solver_max_range);
4100 let baseline_result = baseline_solver.solve()?;
4101
4102 let target_distance = params.target_distance.unwrap_or(baseline_result.max_range);
4104
4105 let baseline_at_target = baseline_result
4107 .position_at_range(target_distance)
4108 .ok_or("Could not interpolate baseline at target distance")?;
4109
4110 let velocity_delta_dist = Normal::new(0.0, params.velocity_std_dev)
4115 .map_err(|e| format!("Invalid velocity distribution: {}", e))?;
4116 let angle_dist = Normal::new(base_inputs.muzzle_angle, params.angle_std_dev)
4117 .map_err(|e| format!("Invalid angle distribution: {}", e))?;
4118 let bc_dist = Normal::new(base_inputs.bc_value, params.bc_std_dev)
4119 .map_err(|e| format!("Invalid BC distribution: {}", e))?;
4120 let wind_sampler = MonteCarloWindSampler::new(
4123 base_wind,
4124 params.wind_speed_std_dev,
4125 wind_direction_std_dev,
4126 )?;
4127 let azimuth_dist = Normal::new(base_inputs.azimuth_angle, params.azimuth_std_dev)
4128 .map_err(|e| format!("Invalid azimuth distribution: {}", e))?;
4129
4130 Ok(Self {
4131 base_inputs,
4132 atmosphere,
4133 solver_max_range,
4134 target_distance,
4135 baseline_at_target,
4136 velocity_delta_dist,
4137 angle_dist,
4138 bc_dist,
4139 wind_sampler,
4140 azimuth_dist,
4141 })
4142 }
4143
4144 fn sample_one_trial<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> Option<TrialOutcome> {
4154 use rand_distr::Distribution;
4155
4156 let mut inputs = self.base_inputs.clone();
4158 let muzzle_velocity_delta = self.velocity_delta_dist.sample(&mut *rng);
4159 inputs.muzzle_angle = self.angle_dist.sample(&mut *rng);
4160 inputs.bc_value = self.bc_dist.sample(&mut *rng).max(0.01);
4161 inputs.azimuth_angle = self.azimuth_dist.sample(&mut *rng); let wind = self.wind_sampler.sample(&mut *rng);
4165
4166 let mut solver = TrajectorySolver::new(inputs, wind, self.atmosphere.clone());
4170 solver.inputs.muzzle_velocity =
4171 (solver.inputs.muzzle_velocity + muzzle_velocity_delta).max(0.0);
4172 solver.set_max_range(self.solver_max_range);
4173 let result = solver.solve().ok()?;
4175
4176 let impact_position = if result.max_range < self.target_distance {
4182 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)
4185 } else {
4186 let pos_at_target = result.position_at_range(self.target_distance)?;
4188 Vector3::new(
4193 0.0,
4194 pos_at_target.y - self.baseline_at_target.y,
4195 pos_at_target.z - self.baseline_at_target.z,
4196 )
4197 };
4198
4199 Some(TrialOutcome {
4200 range: result.max_range,
4201 impact_velocity: result.impact_velocity,
4202 impact_position,
4203 })
4204 }
4205}
4206
4207pub fn run_monte_carlo(
4209 base_inputs: BallisticInputs,
4210 params: MonteCarloParams,
4211) -> Result<MonteCarloResults, BallisticsError> {
4212 run_monte_carlo_with_direction_std_dev(base_inputs, params, 0.0)
4213}
4214
4215pub fn run_monte_carlo_with_direction_std_dev(
4220 base_inputs: BallisticInputs,
4221 params: MonteCarloParams,
4222 wind_direction_std_dev: f64,
4223) -> Result<MonteCarloResults, BallisticsError> {
4224 let base_wind = WindConditions {
4225 speed: params.base_wind_speed,
4226 direction: params.base_wind_direction,
4227 vertical_speed: 0.0,
4228 };
4229 run_monte_carlo_with_wind_and_direction_std_dev(
4230 base_inputs,
4231 base_wind,
4232 params,
4233 wind_direction_std_dev,
4234 )
4235}
4236
4237pub fn run_monte_carlo_with_wind(
4239 base_inputs: BallisticInputs,
4240 base_wind: WindConditions,
4241 params: MonteCarloParams,
4242) -> Result<MonteCarloResults, BallisticsError> {
4243 run_monte_carlo_with_wind_and_direction_std_dev(base_inputs, base_wind, params, 0.0)
4244}
4245
4246pub fn run_monte_carlo_with_wind_and_direction_std_dev(
4251 base_inputs: BallisticInputs,
4252 base_wind: WindConditions,
4253 params: MonteCarloParams,
4254 wind_direction_std_dev: f64,
4255) -> Result<MonteCarloResults, BallisticsError> {
4256 let mut rng = rand::rng();
4257 run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
4258 base_inputs,
4259 base_wind,
4260 params,
4261 wind_direction_std_dev,
4262 &mut rng,
4263 )
4264}
4265
4266pub fn run_monte_carlo_with_wind_and_direction_std_dev_seeded(
4273 base_inputs: BallisticInputs,
4274 base_wind: WindConditions,
4275 params: MonteCarloParams,
4276 wind_direction_std_dev: f64,
4277 seed: u64,
4278) -> Result<MonteCarloResults, BallisticsError> {
4279 use rand::{rngs::StdRng, SeedableRng};
4280 let mut rng = StdRng::seed_from_u64(seed);
4281 run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
4282 base_inputs,
4283 base_wind,
4284 params,
4285 wind_direction_std_dev,
4286 &mut rng,
4287 )
4288}
4289
4290fn run_monte_carlo_with_wind_and_direction_std_dev_using_rng<R: rand::Rng + ?Sized>(
4291 base_inputs: BallisticInputs,
4292 base_wind: WindConditions,
4293 params: MonteCarloParams,
4294 wind_direction_std_dev: f64,
4295 rng: &mut R,
4296) -> Result<MonteCarloResults, BallisticsError> {
4297 let mut ranges = Vec::new();
4298 let mut impact_velocities = Vec::new();
4299 let mut impact_positions = Vec::new();
4300
4301 let sampler = MonteCarloTrialSampler::new(
4302 base_inputs,
4303 &base_wind,
4304 ¶ms,
4305 wind_direction_std_dev,
4306 )?;
4307
4308 for _ in 0..params.num_simulations {
4309 if let Some(outcome) = sampler.sample_one_trial(rng) {
4313 ranges.push(outcome.range);
4314 impact_velocities.push(outcome.impact_velocity);
4315 impact_positions.push(outcome.impact_position);
4316 }
4317 }
4318
4319 if ranges.is_empty() {
4320 return Err("No successful simulations".into());
4321 }
4322
4323 Ok(MonteCarloResults {
4324 ranges,
4325 impact_velocities,
4326 impact_positions,
4327 })
4328}
4329
4330pub const MC_ADAPTIVE_SCHEMA_VERSION_V1: u32 = 1;
4332
4333pub const MC_ADAPTIVE_METHOD_V1: &str = "anytime_beta_binomial_mixture_cs_v1";
4340
4341pub const MC_ADAPTIVE_ASSUMPTIONS_V1: [&str; 4] = [
4355 "Sampling uncertainty only: intervals cover Monte Carlo sampling error, not model error in the trajectory solver or its inputs.",
4356 "Anytime-valid stopping: the beta-binomial mixture confidence sequence keeps its coverage guarantee despite stopping the moment the target half-width is met.",
4357 "Input dispersions are the independent normal distributions declared in MonteCarloParams; correlations between inputs are not modeled.",
4358 "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.",
4359];
4360
4361#[derive(Debug, Clone)]
4368pub struct McConvergence {
4369 pub level: ConfidenceLevel,
4371 pub target_half_width: f64,
4374 pub min_samples: u64,
4378 pub max_samples: u64,
4381 pub batch_size: u64,
4384}
4385
4386impl Default for McConvergence {
4387 fn default() -> Self {
4388 Self {
4389 level: ConfidenceLevel::P95,
4390 target_half_width: 0.02,
4391 min_samples: 1_000,
4392 max_samples: 100_000,
4393 batch_size: 500,
4394 }
4395 }
4396}
4397
4398impl McConvergence {
4399 pub fn validate(&self) -> Result<(), String> {
4408 if !self.target_half_width.is_finite() || self.target_half_width <= 0.0 {
4409 return Err(format!(
4410 "McConvergence.target_half_width must be a finite value greater than zero (got {})",
4411 self.target_half_width
4412 ));
4413 }
4414 if self.batch_size == 0 {
4415 return Err("McConvergence.batch_size must be greater than zero".to_string());
4416 }
4417 if self.max_samples == 0 {
4418 return Err("McConvergence.max_samples must be greater than zero".to_string());
4419 }
4420 if self.max_samples < self.min_samples {
4421 return Err(format!(
4422 "McConvergence.max_samples ({}) must be at least McConvergence.min_samples ({})",
4423 self.max_samples, self.min_samples
4424 ));
4425 }
4426 Ok(())
4427 }
4428}
4429
4430#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize)]
4432#[serde(rename_all = "snake_case")]
4433pub enum McStopReason {
4434 TargetHalfWidthMet,
4437 MaxSamplesReached,
4440}
4441
4442#[derive(Debug, Clone, serde::Serialize)]
4449pub struct AdaptiveMcReportV1 {
4450 pub schema_version: u32,
4452 pub method: String,
4454 pub assumptions: Vec<String>,
4456 pub confidence_percent: u32,
4458 pub hit_probability: f64,
4461 pub ci_low: f64,
4463 pub ci_high: f64,
4465 pub samples: u64,
4467 pub attempts: u64,
4478 pub arrivals: u64,
4493 pub stop_reason: McStopReason,
4495 pub hit_radius_m: f64,
4497 pub target_distance_m: f64,
4499 pub mean_impact_velocity_mps: f64,
4504 pub std_impact_velocity_mps: f64,
4506 pub mean_drop_at_target_m: f64,
4514 pub std_drop_at_target_m: f64,
4517 pub mean_wind_drift_at_target_m: f64,
4520 pub std_wind_drift_at_target_m: f64,
4523}
4524
4525pub fn run_monte_carlo_adaptive_seeded(
4600 base_inputs: &BallisticInputs,
4601 base_wind: &WindConditions,
4602 params: &MonteCarloParams,
4603 convergence: &McConvergence,
4604 hit_radius_m: f64,
4605 seed: u64,
4606) -> Result<AdaptiveMcReportV1, String> {
4607 use rand::{rngs::StdRng, SeedableRng};
4608
4609 convergence.validate()?;
4610
4611 let sampler = MonteCarloTrialSampler::new(base_inputs.clone(), base_wind, params, 0.0)
4614 .map_err(|e| e.to_string())?;
4615
4616 let mut rng = StdRng::seed_from_u64(seed);
4617 let mut hits_cs = BernoulliConfidenceSequence::new(convergence.level);
4618 let mut impact_velocity = Welford::new();
4619 let mut drop_at_target = Welford::new();
4620 let mut drift_at_target = Welford::new();
4621
4622 let mut attempts: u64 = 0;
4623 let mut stop_reason = McStopReason::MaxSamplesReached;
4624
4625 while attempts < convergence.max_samples {
4626 let batch = convergence.batch_size.min(convergence.max_samples - attempts);
4629 let mut batch_hits: u64 = 0;
4630 let mut batch_trials: u64 = 0;
4631
4632 for _ in 0..batch {
4633 attempts += 1;
4634 let Some(outcome) = sampler.sample_one_trial(&mut rng) else {
4635 continue; };
4637 batch_trials += 1;
4638 if MonteCarloResults::position_is_hit(&outcome.impact_position, hit_radius_m) {
4639 batch_hits += 1;
4640 }
4641 if MonteCarloResults::position_reached_target(&outcome.impact_position) {
4645 impact_velocity.push(outcome.impact_velocity);
4646 drop_at_target.push(outcome.impact_position.y);
4647 drift_at_target.push(outcome.impact_position.z);
4648 }
4649 }
4650
4651 hits_cs.update_batch(batch_hits, batch_trials);
4652
4653 if hits_cs.trials() >= convergence.min_samples
4654 && hits_cs.half_width() <= convergence.target_half_width
4655 {
4656 stop_reason = McStopReason::TargetHalfWidthMet;
4657 break;
4658 }
4659 }
4660
4661 let samples = hits_cs.trials();
4662 if samples == 0 {
4663 return Err("No successful simulations".to_string());
4664 }
4665 let (ci_low, ci_high) = hits_cs.bounds();
4666
4667 Ok(AdaptiveMcReportV1 {
4668 schema_version: MC_ADAPTIVE_SCHEMA_VERSION_V1,
4669 method: MC_ADAPTIVE_METHOD_V1.to_string(),
4670 assumptions: MC_ADAPTIVE_ASSUMPTIONS_V1
4671 .iter()
4672 .map(|s| s.to_string())
4673 .collect(),
4674 confidence_percent: convergence.level.as_percent(),
4675 hit_probability: hits_cs.successes() as f64 / samples as f64,
4676 ci_low,
4677 ci_high,
4678 samples,
4679 attempts,
4680 arrivals: drop_at_target.count(),
4684 stop_reason,
4685 hit_radius_m,
4686 target_distance_m: sampler.target_distance,
4687 mean_impact_velocity_mps: impact_velocity.mean(),
4688 std_impact_velocity_mps: impact_velocity.sample_std(),
4689 mean_drop_at_target_m: drop_at_target.mean(),
4690 std_drop_at_target_m: drop_at_target.sample_std(),
4691 mean_wind_drift_at_target_m: drift_at_target.mean(),
4692 std_wind_drift_at_target_m: drift_at_target.sample_std(),
4693 })
4694}
4695
4696pub fn calculate_zero_angle(
4698 inputs: BallisticInputs,
4699 target_distance: f64,
4700 target_height: f64,
4701) -> Result<f64, BallisticsError> {
4702 calculate_zero_angle_with_conditions(
4703 inputs,
4704 target_distance,
4705 target_height,
4706 WindConditions::default(),
4707 AtmosphericConditions::default(),
4708 )
4709}
4710
4711pub fn calculate_zero_angle_with_conditions(
4712 inputs: BallisticInputs,
4713 target_distance: f64,
4714 target_height: f64,
4715 wind: WindConditions,
4716 atmosphere: AtmosphericConditions,
4717) -> Result<f64, BallisticsError> {
4718 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
4719 solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
4720}
4721
4722pub fn calculate_zero_angle_with_resolved_conditions(
4728 inputs: BallisticInputs,
4729 target_distance: f64,
4730 target_height: f64,
4731 wind: WindConditions,
4732 atmosphere: AtmosphericConditions,
4733) -> Result<f64, BallisticsError> {
4734 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
4735 solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
4736}
4737
4738pub const ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M: f64 = 2000.0;
4746
4747pub fn calculate_zero_range_from_angle_with_conditions(
4762 inputs: BallisticInputs,
4763 zero_angle_rad: f64,
4764 target_height: f64,
4765 wind: WindConditions,
4766 atmosphere: AtmosphericConditions,
4767) -> Result<ZeroCrossings, BallisticsError> {
4768 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
4769 solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
4770 solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
4771}
4772
4773pub fn calculate_zero_range_from_angle_with_resolved_conditions(
4779 inputs: BallisticInputs,
4780 zero_angle_rad: f64,
4781 target_height: f64,
4782 wind: WindConditions,
4783 atmosphere: AtmosphericConditions,
4784) -> Result<ZeroCrossings, BallisticsError> {
4785 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
4786 solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
4787 solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
4788}
4789
4790#[derive(Debug, Clone, Copy, PartialEq, Eq)]
4792pub enum BcFitMode {
4793 Drop,
4795 Velocity,
4798}
4799
4800#[derive(Debug, Clone, Copy)]
4802pub struct BcEstimate {
4803 pub bc: f64,
4805 pub rms_error: f64,
4807 pub drag_model: DragModel,
4809 pub mode: BcFitMode,
4811 pub at_bound: bool,
4815}
4816
4817fn fit_value_at(
4825 points: &[TrajectoryPoint],
4826 target_dist: f64,
4827 mode: BcFitMode,
4828 drop_offset: f64,
4829) -> Option<f64> {
4830 let val = |p: &TrajectoryPoint| match mode {
4831 BcFitMode::Drop => drop_offset - p.position.y,
4832 BcFitMode::Velocity => p.velocity_magnitude,
4833 };
4834 for i in 0..points.len() {
4835 if points[i].position.x >= target_dist {
4836 if i == 0 {
4837 return Some(val(&points[0]));
4838 }
4839 let p1 = &points[i - 1];
4840 let p2 = &points[i];
4841 let dx = p2.position.x - p1.position.x;
4842 if dx.abs() < 1e-9 {
4843 return Some(val(p2));
4844 }
4845 let t = (target_dist - p1.position.x) / dx;
4846 return Some(val(p1) + t * (val(p2) - val(p1)));
4847 }
4848 }
4849 None
4850}
4851
4852fn fit_residual_sse(
4853 trajectory: &[TrajectoryPoint],
4854 observations: &[(f64, f64)],
4855 mode: BcFitMode,
4856 drop_offset: f64,
4857) -> Option<f64> {
4858 if observations.is_empty() {
4859 return None;
4860 }
4861 let mut total = 0.0;
4862 for (target_dist, target_val) in observations {
4863 let value = fit_value_at(trajectory, *target_dist, mode, drop_offset)?;
4866 let error = value - target_val;
4867 total += error * error;
4868 }
4869 Some(total)
4870}
4871
4872#[allow(clippy::too_many_arguments)] pub fn estimate_bc_fit(
4890 velocity: f64,
4891 mass: f64,
4892 diameter: f64,
4893 points: &[(f64, f64)],
4894 drag_model: DragModel,
4895 mode: BcFitMode,
4896 atmosphere: AtmosphericConditions,
4897 zero_range: Option<f64>,
4898 sight_height: f64,
4899) -> Result<BcEstimate, BallisticsError> {
4900 if points.is_empty() {
4901 return Err(BallisticsError::from(
4902 "No data points provided for BC estimation.".to_string(),
4903 ));
4904 }
4905 let max_dist = points.iter().map(|(d, _)| *d).fold(0.0_f64, f64::max);
4906 let drop_offset = if zero_range.is_some() { sight_height } else { 0.0 };
4909
4910 let sse = |bc_value: f64| -> Option<f64> {
4912 let mut inputs = BallisticInputs {
4913 muzzle_velocity: velocity,
4914 bc_value,
4915 bc_type: drag_model,
4916 bullet_mass: mass,
4917 bullet_diameter: diameter,
4918 sight_height,
4919 ..Default::default()
4920 };
4921 if let Some(zr) = zero_range {
4924 let za = calculate_zero_angle_with_conditions(
4930 inputs.clone(),
4931 zr,
4932 sight_height,
4933 WindConditions::default(),
4934 atmosphere.clone(),
4935 )
4936 .ok()?;
4937 inputs.muzzle_angle = za;
4938 }
4939 let mut solver =
4940 TrajectorySolver::new(inputs, WindConditions::default(), atmosphere.clone());
4941 solver.set_max_range(max_dist * 1.5);
4942 let result = solver.solve().ok()?;
4943 fit_residual_sse(&result.points, points, mode, drop_offset)
4944 };
4945
4946 let (bc_min, bc_max) = match drag_model {
4950 DragModel::G7 => (0.05, 0.70),
4951 _ => (0.10, 1.20),
4952 };
4953
4954 let mut best_bc = f64::NAN;
4956 let mut best_sse = f64::MAX;
4957 let mut bc = bc_min;
4958 while bc <= bc_max + 1e-9 {
4959 if let Some(s) = sse(bc) {
4960 if s < best_sse {
4961 best_sse = s;
4962 best_bc = bc;
4963 }
4964 }
4965 bc += 0.01;
4966 }
4967 if !best_bc.is_finite() {
4968 return Err(BallisticsError::from(
4969 "Unable to estimate BC from provided data. Check that the values and units are correct."
4970 .to_string(),
4971 ));
4972 }
4973
4974 let lo = (best_bc - 0.01).max(bc_min);
4976 let hi = (best_bc + 0.01).min(bc_max);
4977 let mut bc = lo;
4978 while bc <= hi + 1e-9 {
4979 if let Some(s) = sse(bc) {
4980 if s < best_sse {
4981 best_sse = s;
4982 best_bc = bc;
4983 }
4984 }
4985 bc += 0.001;
4986 }
4987
4988 let at_bound = best_bc <= bc_min + 0.011 || best_bc >= bc_max - 0.011;
4991 let rms_error = (best_sse / points.len() as f64).sqrt();
4994 Ok(BcEstimate {
4995 bc: best_bc,
4996 rms_error,
4997 drag_model,
4998 mode,
4999 at_bound,
5000 })
5001}
5002
5003pub fn estimate_bc_from_trajectory(
5006 velocity: f64,
5007 mass: f64,
5008 diameter: f64,
5009 points: &[(f64, f64)], ) -> Result<f64, BallisticsError> {
5011 estimate_bc_fit(
5012 velocity,
5013 mass,
5014 diameter,
5015 points,
5016 DragModel::G1,
5017 BcFitMode::Drop,
5018 AtmosphericConditions::default(),
5019 None,
5020 0.05,
5021 )
5022 .map(|e| e.bc)
5023}
5024
5025use rand;
5027use rand_distr;
5028
5029#[cfg(test)]
5030mod mba737_powder_resolution_tests {
5031 use super::*;
5032
5033 #[test]
5034 fn linear_model_cold_powder_subtracts() {
5035 let v = resolve_powder_adjusted_velocity(823.0, 11.1, true, 0.5486, 21.1, None, None);
5037 assert!((v - (823.0 + 0.5486 * (11.1 - 21.1))).abs() < 1e-12);
5038 assert!(v < 823.0);
5039 }
5040
5041 #[test]
5042 fn linear_model_hot_powder_adds() {
5043 let v = resolve_powder_adjusted_velocity(823.0, 31.1, true, 0.5486, 21.1, None, None);
5044 assert!((v - (823.0 + 0.5486 * 10.0)).abs() < 1e-12);
5045 }
5046
5047 #[test]
5048 fn disabled_flag_is_passthrough() {
5049 let v = resolve_powder_adjusted_velocity(823.0, 40.0, false, 0.5486, 21.1, None, None);
5050 assert_eq!(v, 823.0);
5051 }
5052
5053 #[test]
5054 fn curve_overrides_linear_and_interpolates_at_powder_temp() {
5055 let curve = [(4.4, 798.6), (21.1, 823.0), (37.8, 841.2)];
5056 let v = resolve_powder_adjusted_velocity(823.0, 30.0, true, 99.0, 21.1, Some(&curve), Some(4.4));
5058 assert!((v - 798.6).abs() < 1e-9);
5059 }
5060
5061 #[test]
5062 fn curve_falls_back_to_ambient_and_clamps() {
5063 let curve = [(4.4, 798.6), (37.8, 841.2)];
5064 let v = resolve_powder_adjusted_velocity(823.0, -40.0, true, 1.0, 21.1, Some(&curve), None);
5066 assert!((v - 798.6).abs() < 1e-9);
5067 let v_hot = resolve_powder_adjusted_velocity(823.0, 60.0, true, 1.0, 21.1, Some(&curve), None);
5068 assert!((v_hot - 841.2).abs() < 1e-9);
5069 }
5070
5071 #[test]
5072 fn empty_curve_suppresses_linear_fallback() {
5073 let v = resolve_powder_adjusted_velocity(823.0, 40.0, true, 0.5486, 21.1, Some(&[]), None);
5075 assert_eq!(v, 823.0);
5076 }
5077
5078 #[test]
5079 fn sweep_huge_range_errors_instead_of_overflowing() {
5080 assert!(parse_powder_sweep("0:1e20:1").is_err());
5083 assert!(parse_powder_sweep("0:1e308:1e-3").is_err());
5084 }
5085
5086 #[test]
5087 fn sweep_fractional_step_keeps_end_row() {
5088 let rows = parse_powder_sweep("0:0.3:0.1").unwrap();
5090 assert_eq!(rows.len(), 4);
5091 assert!((rows[3] - 0.3).abs() < 1e-9);
5092 }
5093
5094 #[test]
5095 fn solver_and_helper_agree_on_linear_model() {
5096 let inputs = BallisticInputs {
5098 use_powder_sensitivity: true,
5099 powder_temp_sensitivity: 0.5486,
5100 powder_temp: 21.1,
5101 temperature: 4.4,
5102 ..Default::default()
5103 };
5104 let expected = resolve_powder_adjusted_velocity(
5105 inputs.muzzle_velocity,
5106 inputs.temperature,
5107 true,
5108 0.5486,
5109 21.1,
5110 None,
5111 None,
5112 );
5113 let solver = TrajectorySolver::new(
5114 inputs,
5115 WindConditions::default(),
5116 AtmosphericConditions::default(),
5117 );
5118 assert!((solver.inputs.muzzle_velocity - expected).abs() < 1e-12);
5119 }
5120}
5121
5122#[cfg(test)]
5123mod mba1302_solver_seam_tests {
5124 use super::*;
5125 use crate::wind::WindSegment;
5126
5127 #[test]
5128 fn authoritative_station_atmosphere_preserves_explicit_standard_values_at_altitude() {
5129 let atmosphere = AtmosphericConditions {
5130 temperature: 15.0,
5131 pressure: 1013.25,
5132 humidity: 50.0,
5133 altitude: 2_000.0,
5134 };
5135 let legacy = TrajectorySolver::new(
5136 BallisticInputs::default(),
5137 WindConditions::default(),
5138 atmosphere.clone(),
5139 );
5140 let authoritative = TrajectorySolver::new_with_resolved_station_atmosphere(
5141 BallisticInputs::default(),
5142 WindConditions::default(),
5143 atmosphere,
5144 );
5145
5146 let (legacy_density, _, legacy_temp_c, legacy_pressure_hpa) = legacy.resolved_atmosphere();
5147 let (authoritative_density, _, authoritative_temp_c, authoritative_pressure_hpa) =
5148 authoritative.resolved_atmosphere();
5149 let (icao_temp_k, icao_pressure_pa) =
5150 crate::atmosphere::calculate_icao_standard_atmosphere(2_000.0);
5151 let (expected_authoritative_density, _) =
5152 crate::atmosphere::calculate_atmosphere(2_000.0, Some(15.0), Some(1013.25), 50.0);
5153
5154 assert!((legacy_temp_c - (icao_temp_k - 273.15)).abs() < 1e-12);
5155 assert!((legacy_pressure_hpa - icao_pressure_pa / 100.0).abs() < 1e-12);
5156 assert_eq!(authoritative_temp_c.to_bits(), 15.0_f64.to_bits());
5157 assert_eq!(authoritative_pressure_hpa.to_bits(), 1013.25_f64.to_bits());
5158 assert_eq!(
5159 authoritative_density.to_bits(),
5160 expected_authoritative_density.to_bits()
5161 );
5162 assert!(
5163 (authoritative_density - legacy_density).abs() > 0.1,
5164 "explicit standard values at altitude must differ from ICAO-at-altitude: explicit={authoritative_density}, ICAO={legacy_density}"
5165 );
5166 }
5167
5168 #[test]
5177 fn precomputed_absolute_resolution_via_authoritative_matches_legacy_new() {
5178 for (temperature, pressure, altitude) in [
5179 (15.0, 1013.25, 0.0), (15.0, 1013.25, 2000.0), (-5.0, 850.0, 2000.0), (22.0, 950.0, 500.0),
5183 ] {
5184 let atmosphere = AtmosphericConditions {
5185 temperature,
5186 pressure,
5187 humidity: 50.0,
5188 altitude,
5189 };
5190 let legacy = TrajectorySolver::new(
5191 BallisticInputs::default(),
5192 WindConditions::default(),
5193 atmosphere.clone(),
5194 );
5195
5196 let (resolved_temp_c, resolved_pressure_hpa) =
5197 crate::atmosphere::resolve_station_conditions_with_pressure_mode(
5198 temperature,
5199 pressure,
5200 altitude,
5201 crate::atmosphere::PressureReferenceMode::Absolute,
5202 );
5203 let precomputed_atmosphere = AtmosphericConditions {
5204 temperature: resolved_temp_c,
5205 pressure: resolved_pressure_hpa,
5206 humidity: 50.0,
5207 altitude,
5208 };
5209 let precomputed = TrajectorySolver::new_with_resolved_station_atmosphere(
5210 BallisticInputs::default(),
5211 WindConditions::default(),
5212 precomputed_atmosphere,
5213 );
5214
5215 let (legacy_density, legacy_sos, legacy_temp_c, legacy_pressure_hpa) =
5216 legacy.resolved_atmosphere();
5217 let (pre_density, pre_sos, pre_temp_c, pre_pressure_hpa) =
5218 precomputed.resolved_atmosphere();
5219
5220 assert_eq!(
5221 legacy_temp_c.to_bits(),
5222 pre_temp_c.to_bits(),
5223 "temperature=({temperature}, {pressure}, {altitude})"
5224 );
5225 assert_eq!(
5226 legacy_pressure_hpa.to_bits(),
5227 pre_pressure_hpa.to_bits(),
5228 "pressure=({temperature}, {pressure}, {altitude})"
5229 );
5230 assert_eq!(legacy_density.to_bits(), pre_density.to_bits());
5231 assert_eq!(legacy_sos.to_bits(), pre_sos.to_bits());
5232 }
5233 }
5234
5235 fn configured_euler_zero(vertical_wind_mps: f64, time_step_s: f64) -> TrajectorySolver {
5236 let inputs = BallisticInputs {
5237 muzzle_velocity: 800.0,
5238 bc_value: 0.5,
5239 bc_type: DragModel::G7,
5240 bullet_mass: 0.0109,
5241 bullet_diameter: 0.00782,
5242 bullet_length: 0.0309,
5243 sight_height: 0.05,
5244 ground_threshold: -100.0,
5245 use_rk4: false,
5246 use_adaptive_rk45: false,
5247 ..BallisticInputs::default()
5248 };
5249 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(
5250 inputs,
5251 WindConditions::default(),
5252 AtmosphericConditions::default(),
5253 );
5254 solver.set_max_range(300.0);
5255 solver.set_time_step(time_step_s);
5256 if vertical_wind_mps != 0.0 {
5257 solver.set_wind_segments(vec![WindSegment {
5258 speed_kmh: 0.0,
5259 angle_deg: 0.0,
5260 until_m: 400.0,
5261 vertical_mps: vertical_wind_mps,
5262 }]);
5263 }
5264 solver
5265 }
5266
5267 #[test]
5268 fn inclined_shot_zeroes_like_a_level_rifle() {
5269 const ZERO_DISTANCE_M: f64 = 91.44; const SIGHT_HEIGHT_M: f64 = 0.0381; let inputs = BallisticInputs {
5278 bc_value: 0.5,
5279 bullet_mass: 150.0 * 0.06479891 / 1000.0,
5280 muzzle_velocity: 2700.0 * 0.3048,
5281 sight_height: SIGHT_HEIGHT_M,
5282 ..Default::default()
5283 };
5284
5285 let mut level = inputs.clone();
5286 level.shooting_angle = 0.0;
5287 let level_angle = TrajectorySolver::new(level, Default::default(), Default::default())
5288 .find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
5289 .expect("level zero must solve");
5290
5291 let mut inclined = inputs;
5292 inclined.shooting_angle = 5.71_f64.to_radians();
5293 let inclined_angle =
5294 TrajectorySolver::new(inclined, Default::default(), Default::default())
5295 .find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
5296 .expect("MBA-1412: a 5.71 deg incline at a 100 yd zero must be solvable");
5297
5298 assert!(
5299 (inclined_angle - level_angle).abs() < 1e-9,
5300 "zeroing is level-rifle sight geometry; incline must not move the solved zero: \
5301 level={level_angle}, inclined={inclined_angle}"
5302 );
5303 }
5304
5305 #[test]
5306 fn configured_zero_keeps_segments_method_and_time_step_then_sets_base_angle() {
5307 const TARGET_DISTANCE_M: f64 = 150.0;
5308 const TARGET_HEIGHT_M: f64 = 0.05;
5309
5310 let mut segmented = configured_euler_zero(-10.0, 0.02);
5313 let coarse_height = segmented
5314 .zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
5315 .expect("coarse configured trial")
5316 .expect("coarse trial reaches target");
5317 let mut fine = segmented.clone();
5318 fine.set_time_step(0.001);
5319 let fine_height = fine
5320 .zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
5321 .expect("fine configured trial")
5322 .expect("fine trial reaches target");
5323 assert!(
5324 (coarse_height - fine_height).abs() > 1e-5,
5325 "configured Euler step must affect zero trials: coarse={coarse_height}, fine={fine_height}"
5326 );
5327
5328 let segmented_angle = segmented
5329 .calculate_and_set_zero_angle(TARGET_DISTANCE_M, TARGET_HEIGHT_M, ZeroTargetFrame::SightLine)
5330 .expect("segmented zero");
5331 assert_eq!(
5332 segmented.inputs.muzzle_angle.to_bits(),
5333 segmented_angle.to_bits(),
5334 "successful zero must install its angle on the configured solver"
5335 );
5336 assert_eq!(segmented.time_step.to_bits(), 0.02_f64.to_bits());
5337 assert_eq!(segmented.max_range.to_bits(), 300.0_f64.to_bits());
5338 assert!(segmented.wind_sock.is_some());
5339 assert_eq!(
5340 segmented.station_atmosphere_resolution,
5341 StationAtmosphereResolution::Authoritative
5342 );
5343 let zero_height = segmented
5344 .zero_trial_height_at(segmented_angle, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
5345 .expect("verify segmented zero")
5346 .expect("zeroed trial reaches target");
5347 assert!(
5348 (zero_height - TARGET_HEIGHT_M).abs() < 0.0001,
5349 "configured zero missed target: height={zero_height}"
5350 );
5351
5352 let mut calm = configured_euler_zero(0.0, 0.02);
5353 let calm_angle = calm
5354 .calculate_and_set_zero_angle(TARGET_DISTANCE_M, TARGET_HEIGHT_M, ZeroTargetFrame::SightLine)
5355 .expect("calm zero");
5356 assert!(
5357 (segmented_angle - calm_angle).abs() > 1e-5,
5358 "segmented vertical wind must participate in zero trials: segmented={segmented_angle}, calm={calm_angle}"
5359 );
5360 }
5361}
5362
5363#[cfg(test)]
5364mod result_sanity_tests {
5365 use super::*;
5366
5367 fn default_solver() -> TrajectorySolver {
5368 TrajectorySolver::new(
5369 BallisticInputs::default(),
5370 WindConditions::default(),
5371 AtmosphericConditions::default(),
5372 )
5373 }
5374
5375 fn minimal_result() -> TrajectoryResult {
5376 TrajectoryResult {
5377 max_range: 100.0,
5378 max_height: 1.0,
5379 time_of_flight: 0.5,
5380 impact_velocity: 700.0,
5381 impact_energy: 2450.0,
5382 projectile_mass_kg: 0.01,
5383 line_of_sight_height_m: 1.5,
5384 station_speed_of_sound_mps: 340.0,
5385 termination: TrajectoryTermination::MaxRange,
5386 points: vec![],
5387 sampled_points: None,
5388 min_pitch_damping: None,
5389 transonic_mach: None,
5390 angular_state: None,
5391 max_yaw_angle: None,
5392 max_precession_angle: None,
5393 aerodynamic_jump: None,
5394 mach_1_2_distance_m: None,
5395 mach_1_0_distance_m: None,
5396 mach_0_9_distance_m: None,
5397 }
5398 }
5399
5400 #[test]
5401 fn mba1293_negative_scalars_fail_the_result_postcondition() {
5402 let solver = default_solver();
5403 solver
5404 .validate_result_sanity(&minimal_result())
5405 .expect("a sane result must pass");
5406
5407 for (name, mutate) in [
5408 ("max_range", (|r| r.max_range = -50.588) as fn(&mut TrajectoryResult)),
5409 ("time_of_flight", |r| r.time_of_flight = -1.0),
5410 ("impact_velocity", |r| r.impact_velocity = -700.0),
5411 ("impact_energy", |r| r.impact_energy = -1.0),
5412 ] {
5413 let mut result = minimal_result();
5414 mutate(&mut result);
5415 let error = solver
5416 .validate_result_sanity(&result)
5417 .expect_err("negative scalar must fail");
5418 assert!(
5419 error.to_string().contains(name),
5420 "error for {name} did not name the field: {error}"
5421 );
5422 }
5423 }
5424
5425 #[test]
5426 fn mba1293_speed_budget_bounds_legitimate_states_and_rejects_divergence() {
5427 let solver = default_solver();
5428 let mv = solver.inputs.muzzle_velocity;
5429
5430 let position = Vector3::new(10.0, 0.0, 0.0);
5432 solver
5433 .validate_integration_state(&position, &Vector3::new(mv, 0.0, 0.0), 0.01)
5434 .expect("muzzle-speed state must pass");
5435
5436 let error = solver
5438 .validate_integration_state(&position, &Vector3::new(-13.0 * mv, 0.0, 0.0), 0.01)
5439 .expect_err("13x muzzle speed must fail the budget");
5440 assert!(error.to_string().contains("diverged"), "{error}");
5441
5442 let after_fall = mv + crate::constants::G_ACCEL_MPS2 * 60.0;
5444 solver
5445 .validate_integration_state(&position, &Vector3::new(0.0, -after_fall, 0.0), 60.0)
5446 .expect("gravity-accelerated speed within g*t must pass");
5447 }
5448}
5449
5450#[cfg(test)]
5451mod trajectory_point_budget_tests {
5452 use super::*;
5453 use crate::MAX_TRAJECTORY_SAMPLES;
5454
5455 fn solver_with_budget(
5456 use_rk4: bool,
5457 use_adaptive_rk45: bool,
5458 point_budget: usize,
5459 max_range: f64,
5460 ) -> TrajectorySolver {
5461 let inputs = BallisticInputs {
5462 use_rk4,
5463 use_adaptive_rk45,
5464 ground_threshold: f64::NEG_INFINITY,
5465 ..BallisticInputs::default()
5466 };
5467 let mut solver = TrajectorySolver::new(
5468 inputs,
5469 WindConditions::default(),
5470 AtmosphericConditions::default(),
5471 );
5472 solver.max_trajectory_points = point_budget;
5473 solver.set_max_range(max_range);
5474 solver.set_time_step(0.001);
5475 solver
5476 }
5477
5478 #[test]
5479 fn mba1283_every_solver_errors_instead_of_exceeding_point_budget() {
5480 for (mode, use_rk4, use_adaptive_rk45) in [
5481 ("Euler", false, false),
5482 ("RK4", true, false),
5483 ("RK45", true, true),
5484 ] {
5485 let error = solver_with_budget(use_rk4, use_adaptive_rk45, 3, 10.0)
5486 .solve()
5487 .expect_err("a solve requiring more than three points must fail");
5488 assert!(
5489 error.to_string().contains("point limit of 3"),
5490 "unexpected {mode} point-budget error: {error}"
5491 );
5492 }
5493 }
5494
5495 #[test]
5496 fn mba1283_interpolated_endpoint_counts_toward_point_budget() {
5497 for (mode, use_rk4, use_adaptive_rk45) in [
5498 ("Euler", false, false),
5499 ("RK4", true, false),
5500 ("RK45", true, true),
5501 ] {
5502 let result = solver_with_budget(use_rk4, use_adaptive_rk45, 2, 0.1)
5503 .solve()
5504 .expect("the initial point plus exact endpoint fit a two-point budget");
5505 assert_eq!(result.points.len(), 2, "unexpected {mode} point count");
5506
5507 let error = solver_with_budget(use_rk4, use_adaptive_rk45, 1, 0.1)
5508 .solve()
5509 .expect_err("the exact endpoint must not exceed a one-point budget");
5510 assert!(
5511 error.to_string().contains("point limit of 1"),
5512 "unexpected {mode} endpoint-budget error: {error}"
5513 );
5514 }
5515 }
5516
5517 #[test]
5518 fn mba1299_every_solver_preflights_the_sample_budget() {
5519 for (mode, use_rk4, use_adaptive_rk45) in [
5520 ("Euler", false, false),
5521 ("RK4", true, false),
5522 ("RK45", true, true),
5523 ] {
5524 let inputs = BallisticInputs {
5525 use_rk4,
5526 use_adaptive_rk45,
5527 enable_trajectory_sampling: true,
5528 sample_interval: 1.0,
5529 ground_threshold: f64::NEG_INFINITY,
5530 ..BallisticInputs::default()
5531 };
5532 let mut solver = TrajectorySolver::new(
5533 inputs,
5534 WindConditions::default(),
5535 AtmosphericConditions::default(),
5536 );
5537 solver.set_max_range(MAX_TRAJECTORY_SAMPLES as f64);
5538 solver.max_trajectory_points = 0;
5541
5542 let error = solver
5543 .solve()
5544 .expect_err("an over-limit sample grid must fail before integration");
5545 assert!(
5546 error
5547 .to_string()
5548 .contains("trajectory sample limit of 250000 exceeded"),
5549 "unexpected {mode} sample-budget error: {error}"
5550 );
5551 }
5552 }
5553
5554 #[test]
5555 fn mba1299_normal_sampling_does_not_change_solver_results() {
5556 for (mode, use_rk4, use_adaptive_rk45) in [
5557 ("Euler", false, false),
5558 ("RK4", true, false),
5559 ("RK45", true, true),
5560 ] {
5561 let solve = |enable_trajectory_sampling| {
5562 let inputs = BallisticInputs {
5563 use_rk4,
5564 use_adaptive_rk45,
5565 enable_trajectory_sampling,
5566 sample_interval: 0.5,
5567 ground_threshold: f64::NEG_INFINITY,
5568 ..BallisticInputs::default()
5569 };
5570 let mut solver = TrajectorySolver::new(
5571 inputs,
5572 WindConditions::default(),
5573 AtmosphericConditions::default(),
5574 );
5575 solver.set_max_range(2.0);
5576 solver.solve().expect("normal short-range solve")
5577 };
5578
5579 let baseline = solve(false);
5580 let sampled = solve(true);
5581 for (field, left, right) in [
5582 ("max_range", baseline.max_range, sampled.max_range),
5583 ("max_height", baseline.max_height, sampled.max_height),
5584 (
5585 "time_of_flight",
5586 baseline.time_of_flight,
5587 sampled.time_of_flight,
5588 ),
5589 (
5590 "impact_velocity",
5591 baseline.impact_velocity,
5592 sampled.impact_velocity,
5593 ),
5594 (
5595 "impact_energy",
5596 baseline.impact_energy,
5597 sampled.impact_energy,
5598 ),
5599 ] {
5600 assert_eq!(
5601 left.to_bits(),
5602 right.to_bits(),
5603 "{mode} sampling changed {field}"
5604 );
5605 }
5606 assert_eq!(baseline.points.len(), sampled.points.len());
5607 for (index, (left, right)) in baseline
5608 .points
5609 .iter()
5610 .zip(&sampled.points)
5611 .enumerate()
5612 {
5613 assert_eq!(left.time.to_bits(), right.time.to_bits(), "{mode} point {index}");
5614 assert_eq!(
5615 left.position.map(f64::to_bits),
5616 right.position.map(f64::to_bits),
5617 "{mode} point {index} position"
5618 );
5619 assert_eq!(
5620 left.velocity_magnitude.to_bits(),
5621 right.velocity_magnitude.to_bits(),
5622 "{mode} point {index} velocity"
5623 );
5624 assert_eq!(
5625 left.kinetic_energy.to_bits(),
5626 right.kinetic_energy.to_bits(),
5627 "{mode} point {index} energy"
5628 );
5629 }
5630 assert!(baseline.sampled_points.is_none());
5631 let samples = sampled
5632 .sampled_points
5633 .expect("sampling-enabled solve should return observations");
5634 assert_eq!(
5635 samples
5636 .iter()
5637 .map(|sample| sample.distance_m)
5638 .collect::<Vec<_>>(),
5639 vec![0.0, 0.5, 1.0, 1.5, 2.0],
5640 "{mode} normal sampling grid changed"
5641 );
5642 }
5643 }
5644}
5645
5646#[cfg(test)]
5647mod monte_carlo_result_tests {
5648 use super::*;
5649
5650 fn make_results(impact_positions: Vec<Vector3<f64>>) -> MonteCarloResults {
5651 let count = impact_positions.len();
5652 MonteCarloResults {
5653 ranges: vec![500.0; count],
5654 impact_velocities: vec![300.0; count],
5655 impact_positions,
5656 }
5657 }
5658
5659 #[test]
5660 fn target_plane_cep_excludes_shortfall_markers() {
5661 let mut positions: Vec<Vector3<f64>> = (1..=5)
5662 .map(|radius| Vector3::new(0.0, radius as f64, 0.0))
5663 .collect();
5664 positions.extend(
5665 (0..5).map(|_| Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)),
5666 );
5667 let results = make_results(positions);
5668
5669 assert_eq!(results.target_arrival_count(), 5);
5670 assert_eq!(results.target_shortfall_fraction(), 0.5);
5671 assert_eq!(results.target_plane_cep(), Some(3.0));
5672
5673 let one_shortfall = make_results(vec![
5674 Vector3::new(0.0, 1.0, 0.0),
5675 Vector3::new(0.0, 2.0, 0.0),
5676 Vector3::new(0.0, 3.0, 0.0),
5677 Vector3::new(0.0, 4.0, 0.0),
5678 Vector3::new(0.0, 5.0, 0.0),
5679 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5680 ]);
5681 assert_eq!(one_shortfall.target_plane_cep(), Some(3.0));
5682 }
5683
5684 #[test]
5685 fn all_shortfalls_have_no_cep_but_still_count_as_misses() {
5686 let all_shortfalls = make_results(vec![
5687 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5688 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5689 ]);
5690 assert_eq!(all_shortfalls.target_arrival_count(), 0);
5691 assert_eq!(all_shortfalls.target_shortfall_fraction(), 1.0);
5692 assert_eq!(all_shortfalls.target_plane_cep(), None);
5693 assert_eq!(all_shortfalls.hit_probability(0.3), 0.0);
5694
5695 let one_hit_one_shortfall = make_results(vec![
5696 Vector3::new(0.0, 0.1, 0.0),
5697 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5698 ]);
5699 assert_eq!(one_hit_one_shortfall.hit_probability(0.3), 0.5);
5700 }
5701
5702 #[test]
5704 fn rect_hit_probability_checks_independent_axis_halves() {
5705 let results = make_results(vec![
5706 Vector3::new(0.0, 0.1, 0.1),
5708 Vector3::new(0.0, 0.0, 0.2),
5710 Vector3::new(0.0, 0.0, 0.201),
5712 Vector3::new(0.0, 0.301, 0.0),
5714 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5716 ]);
5717 assert!((results.rect_hit_probability(0.4, 0.6) - 0.4).abs() < 1e-12);
5719 }
5720
5721 #[test]
5722 fn rect_hit_probability_matches_circular_hit_probability_for_a_centered_hit() {
5723 let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
5724 assert_eq!(results.rect_hit_probability(0.5, 0.5), 1.0);
5725 assert_eq!(results.hit_probability(0.3), 1.0);
5726 }
5727
5728 #[test]
5729 fn rect_hit_probability_is_zero_for_empty_or_nonpositive_dimensions() {
5730 let empty = make_results(vec![]);
5731 assert_eq!(empty.rect_hit_probability(1.0, 1.0), 0.0);
5732
5733 let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
5734 assert_eq!(results.rect_hit_probability(0.0, 1.0), 0.0);
5735 assert_eq!(results.rect_hit_probability(1.0, 0.0), 0.0);
5736 assert_eq!(results.rect_hit_probability(-1.0, 1.0), 0.0);
5737 }
5738}
5739
5740#[cfg(test)]
5741mod monte_carlo_seeded_tests {
5742 use super::*;
5743
5744 fn seeded_test_fixture() -> (BallisticInputs, WindConditions) {
5753 (
5754 BallisticInputs {
5755 muzzle_velocity: 800.0,
5756 ..BallisticInputs::default()
5757 },
5758 WindConditions::default(),
5759 )
5760 }
5761
5762 fn loose_params() -> MonteCarloParams {
5772 MonteCarloParams {
5773 num_simulations: 1, velocity_std_dev: 3.0,
5775 angle_std_dev: 3.5e-4,
5776 bc_std_dev: 0.01,
5777 wind_speed_std_dev: 1.0,
5778 target_distance: Some(300.0),
5779 base_wind_speed: 0.0,
5780 base_wind_direction: 0.0,
5781 azimuth_std_dev: 3.5e-4,
5782 }
5783 }
5784
5785 fn mixed_arrival_params() -> MonteCarloParams {
5793 MonteCarloParams {
5794 num_simulations: 1,
5795 target_distance: Some(1920.0),
5796 ..MonteCarloParams::default()
5797 }
5798 }
5799
5800 #[test]
5816 fn legacy_seeded_estimates_are_pinned_bit_for_bit() {
5817 let (inputs, wind) = seeded_test_fixture();
5818 let params = MonteCarloParams {
5819 num_simulations: 200,
5820 target_distance: Some(500.0),
5821 ..MonteCarloParams::default()
5822 };
5823
5824 let results = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5825 inputs,
5826 wind,
5827 params,
5828 0.01,
5829 0x1352_5EED,
5830 )
5831 .expect("seeded legacy run");
5832
5833 assert_eq!(results.ranges.len(), 200, "ranges length");
5838 assert_eq!(results.impact_velocities.len(), 200, "impact_velocities length");
5839 assert_eq!(results.impact_positions.len(), 200, "impact_positions length");
5840
5841 assert_eq!(
5843 results.hit_probability(DEFAULT_HIT_RADIUS_M).to_bits(),
5844 0.14_f64.to_bits(),
5845 "hit_probability = {:?}",
5846 results.hit_probability(DEFAULT_HIT_RADIUS_M)
5847 );
5848
5849 let expected_ranges: [f64; 3] =
5850 [1907.972891143359, 1936.408435469319, 1912.8150447617645];
5851 let expected_velocities: [f64; 3] =
5852 [238.6187151542299, 239.91923651600106, 243.14112455427164];
5853 let expected_positions: [(f64, f64, f64); 3] = [
5854 (0.0, -0.0643556039548101, 0.7344970252014579),
5855 (0.0, 0.5769422971539162, 0.27227201756386726),
5856 (0.0, -0.7440425792472842, 0.1541804446282822),
5857 ];
5858
5859 for (i, expected) in expected_ranges.iter().enumerate() {
5860 assert_eq!(
5861 results.ranges[i].to_bits(),
5862 expected.to_bits(),
5863 "ranges[{i}] = {:?}, pinned {expected:?}",
5864 results.ranges[i]
5865 );
5866 }
5867 for (i, expected) in expected_velocities.iter().enumerate() {
5868 assert_eq!(
5869 results.impact_velocities[i].to_bits(),
5870 expected.to_bits(),
5871 "impact_velocities[{i}] = {:?}, pinned {expected:?}",
5872 results.impact_velocities[i]
5873 );
5874 }
5875 for (i, (x, y, z)) in expected_positions.iter().enumerate() {
5876 let actual = results.impact_positions[i];
5877 assert_eq!(actual.x.to_bits(), x.to_bits(), "impact_positions[{i}].x = {:?}", actual.x);
5878 assert_eq!(actual.y.to_bits(), y.to_bits(), "impact_positions[{i}].y = {:?}", actual.y);
5879 assert_eq!(actual.z.to_bits(), z.to_bits(), "impact_positions[{i}].z = {:?}", actual.z);
5880 }
5881 }
5882
5883 #[test]
5884 fn seeded_runs_are_deterministic_and_match_the_using_rng_path() {
5885 let inputs = BallisticInputs {
5886 muzzle_velocity: 800.0,
5887 ..BallisticInputs::default()
5888 };
5889 let params = MonteCarloParams {
5890 num_simulations: 64,
5891 target_distance: Some(200.0),
5892 ..MonteCarloParams::default()
5893 };
5894
5895 let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5896 inputs.clone(),
5897 WindConditions::default(),
5898 params.clone(),
5899 0.01,
5900 42,
5901 )
5902 .expect("seeded run a");
5903 let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5904 inputs,
5905 WindConditions::default(),
5906 params,
5907 0.01,
5908 42,
5909 )
5910 .expect("seeded run b");
5911
5912 assert_eq!(a.ranges.len(), b.ranges.len());
5913 for (ra, rb) in a.ranges.iter().zip(b.ranges.iter()) {
5914 assert_eq!(ra.to_bits(), rb.to_bits());
5915 }
5916 for (pa, pb) in a.impact_positions.iter().zip(b.impact_positions.iter()) {
5917 assert_eq!(pa.x.to_bits(), pb.x.to_bits());
5918 assert_eq!(pa.y.to_bits(), pb.y.to_bits());
5919 assert_eq!(pa.z.to_bits(), pb.z.to_bits());
5920 }
5921 }
5922
5923 #[test]
5924 fn different_seeds_generally_produce_different_draws() {
5925 let inputs = BallisticInputs {
5926 muzzle_velocity: 800.0,
5927 ..BallisticInputs::default()
5928 };
5929 let params = MonteCarloParams {
5930 num_simulations: 32,
5931 velocity_std_dev: 5.0,
5932 target_distance: Some(200.0),
5933 ..MonteCarloParams::default()
5934 };
5935
5936 let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5937 inputs.clone(),
5938 WindConditions::default(),
5939 params.clone(),
5940 0.0,
5941 1,
5942 )
5943 .expect("seeded run a");
5944 let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5945 inputs,
5946 WindConditions::default(),
5947 params,
5948 0.0,
5949 2,
5950 )
5951 .expect("seeded run b");
5952
5953 assert_ne!(a.impact_velocities, b.impact_velocities);
5954 }
5955
5956 #[test]
5957 fn adaptive_stops_at_target_half_width_on_an_easy_case() {
5958 let (inputs, wind) = seeded_test_fixture();
5959 let conv = McConvergence {
5960 target_half_width: 0.05,
5961 ..Default::default()
5962 };
5963 let r = run_monte_carlo_adaptive_seeded(
5964 &inputs,
5965 &wind,
5966 &loose_params(),
5967 &conv,
5968 DEFAULT_HIT_RADIUS_M,
5969 0x1352_ADA9,
5970 )
5971 .unwrap();
5972
5973 assert_eq!(r.stop_reason, McStopReason::TargetHalfWidthMet);
5974 assert!(
5975 (r.ci_high - r.ci_low) / 2.0 <= 0.05 + 1e-12,
5976 "half-width {} exceeds the requested 0.05",
5977 (r.ci_high - r.ci_low) / 2.0
5978 );
5979 assert!(r.samples >= conv.min_samples, "stopped below min_samples");
5980 assert!(r.samples < conv.max_samples, "did not actually stop early");
5981 assert!(r.samples.is_multiple_of(conv.batch_size) || r.samples == conv.min_samples);
5982 assert!(r.ci_low <= r.hit_probability && r.hit_probability <= r.ci_high);
5983
5984 assert_eq!(r.hit_radius_m, DEFAULT_HIT_RADIUS_M);
5986 assert_eq!(r.target_distance_m, 300.0);
5987 assert_eq!(r.confidence_percent, 95);
5988 assert!(
5990 r.samples > 1,
5991 "params.num_simulations must be ignored by the adaptive driver"
5992 );
5993 assert!(
5996 r.mean_impact_velocity_mps > 0.0,
5997 "no impact velocity accumulated"
5998 );
5999 assert!(
6000 r.std_drop_at_target_m > 0.0 && r.std_wind_drift_at_target_m > 0.0,
6001 "dispersion collapsed: drop sd {} drift sd {}",
6002 r.std_drop_at_target_m,
6003 r.std_wind_drift_at_target_m
6004 );
6005 }
6006
6007 #[test]
6008 fn adaptive_caps_at_max_samples_on_an_impossible_target() {
6009 let (inputs, wind) = seeded_test_fixture();
6010 let conv = McConvergence {
6011 target_half_width: 1e-6,
6012 max_samples: 3_000,
6013 batch_size: 500,
6014 min_samples: 1_000,
6015 level: ConfidenceLevel::P95,
6016 };
6017 let r = run_monte_carlo_adaptive_seeded(
6018 &inputs,
6019 &wind,
6020 &loose_params(),
6021 &conv,
6022 DEFAULT_HIT_RADIUS_M,
6023 7,
6024 )
6025 .unwrap();
6026
6027 assert_eq!(r.stop_reason, McStopReason::MaxSamplesReached);
6028 assert_eq!(r.samples, 3_000);
6029 }
6030
6031 #[test]
6040 fn adaptive_stops_between_the_floor_and_the_ceiling() {
6041 let (inputs, wind) = seeded_test_fixture();
6042 let conv = McConvergence {
6043 level: ConfidenceLevel::P95,
6044 target_half_width: 0.03,
6045 min_samples: 0,
6046 max_samples: 10_000,
6047 batch_size: 100,
6048 };
6049 let r = run_monte_carlo_adaptive_seeded(
6050 &inputs,
6051 &wind,
6052 &loose_params(),
6053 &conv,
6054 DEFAULT_HIT_RADIUS_M,
6055 0x1352_5A1D,
6056 )
6057 .unwrap();
6058
6059 assert_eq!(r.stop_reason, McStopReason::TargetHalfWidthMet);
6060 assert!(r.samples > 0);
6061 assert!(
6062 r.samples.is_multiple_of(conv.batch_size),
6063 "samples {} is not a whole number of batches",
6064 r.samples
6065 );
6066 assert!(
6067 r.samples > conv.min_samples,
6068 "stopped on the floor, not on the data"
6069 );
6070 assert!(
6071 r.samples < conv.max_samples,
6072 "ran to the ceiling, so nothing adaptive was exercised"
6073 );
6074 assert!(
6075 r.samples > conv.batch_size,
6076 "stopped on the very first look ({} samples); the multi-batch path is untested",
6077 r.samples
6078 );
6079 assert!((r.ci_high - r.ci_low) / 2.0 <= 0.03 + 1e-12);
6080 assert!(r.ci_low <= r.hit_probability && r.hit_probability <= r.ci_high);
6081 assert_eq!(r.attempts, r.samples, "no trial should have been dropped");
6082 }
6083
6084 #[test]
6092 fn adaptive_runs_a_truncated_final_batch_up_to_max_samples() {
6093 let (inputs, wind) = seeded_test_fixture();
6094 let conv = McConvergence {
6095 level: ConfidenceLevel::P95,
6096 target_half_width: 1e-6,
6097 min_samples: 0,
6098 max_samples: 750,
6099 batch_size: 500,
6100 };
6101 let r = run_monte_carlo_adaptive_seeded(
6102 &inputs,
6103 &wind,
6104 &loose_params(),
6105 &conv,
6106 DEFAULT_HIT_RADIUS_M,
6107 0x1352_7B10,
6108 )
6109 .unwrap();
6110
6111 assert_eq!(r.stop_reason, McStopReason::MaxSamplesReached);
6112 assert_eq!(
6113 r.samples, 750,
6114 "the 250-trial final batch did not run, or was not truncated"
6115 );
6116 assert_eq!(r.attempts, 750);
6117 assert!(!r.samples.is_multiple_of(conv.batch_size));
6118 }
6119
6120 #[test]
6121 fn adaptive_is_deterministic_for_a_seed() {
6122 let (inputs, wind) = seeded_test_fixture();
6123 let conv = McConvergence::default();
6124 let a = run_monte_carlo_adaptive_seeded(
6125 &inputs,
6126 &wind,
6127 &loose_params(),
6128 &conv,
6129 DEFAULT_HIT_RADIUS_M,
6130 99,
6131 )
6132 .unwrap();
6133 let b = run_monte_carlo_adaptive_seeded(
6134 &inputs,
6135 &wind,
6136 &loose_params(),
6137 &conv,
6138 DEFAULT_HIT_RADIUS_M,
6139 99,
6140 )
6141 .unwrap();
6142
6143 assert_eq!(a.hit_probability.to_bits(), b.hit_probability.to_bits());
6144 assert_eq!(a.samples, b.samples);
6145 assert_eq!(a.ci_low.to_bits(), b.ci_low.to_bits());
6146 assert_eq!(a.ci_high.to_bits(), b.ci_high.to_bits());
6147 assert_eq!(
6150 a.mean_impact_velocity_mps.to_bits(),
6151 b.mean_impact_velocity_mps.to_bits()
6152 );
6153 assert_eq!(
6154 a.std_drop_at_target_m.to_bits(),
6155 b.std_drop_at_target_m.to_bits()
6156 );
6157 }
6158
6159 #[test]
6160 fn adaptive_report_carries_schema_method_and_all_four_assumptions() {
6161 let (inputs, wind) = seeded_test_fixture();
6162 let conv = McConvergence {
6166 min_samples: 0,
6167 max_samples: 50,
6168 batch_size: 50,
6169 target_half_width: 1.0,
6170 level: ConfidenceLevel::P90,
6171 };
6172 let r = run_monte_carlo_adaptive_seeded(
6173 &inputs,
6174 &wind,
6175 &mixed_arrival_params(),
6176 &conv,
6177 DEFAULT_HIT_RADIUS_M,
6178 0x1352_D0C5,
6179 )
6180 .unwrap();
6181
6182 assert_eq!(r.schema_version, MC_ADAPTIVE_SCHEMA_VERSION_V1);
6183 assert_eq!(r.schema_version, 1);
6184 assert_eq!(r.method, "anytime_beta_binomial_mixture_cs_v1");
6185 assert_eq!(r.confidence_percent, 90);
6186
6187 assert_eq!(r.attempts, 50, "one full batch was drawn");
6193 assert_eq!(r.samples, 50, "no trial was dropped by the solver");
6194 assert!(
6195 r.arrivals > 0 && r.arrivals < r.samples,
6196 "fixture must split the run: arrivals {} of samples {}",
6197 r.arrivals,
6198 r.samples
6199 );
6200 assert!(r.arrivals >= 2, "arrivals {} too few for a sample sd", r.arrivals);
6203 assert!(r.std_drop_at_target_m > 0.0 && r.std_impact_velocity_mps > 0.0);
6204 assert!(r.attempts >= r.samples && r.samples >= r.arrivals);
6206
6207 assert_eq!(r.assumptions.len(), 4, "exactly four assumptions expected");
6213 assert_eq!(
6214 r.assumptions[0],
6215 "Sampling uncertainty only: intervals cover Monte Carlo sampling error, not model error in the trajectory solver or its inputs."
6216 );
6217 assert_eq!(
6218 r.assumptions[1],
6219 "Anytime-valid stopping: the beta-binomial mixture confidence sequence keeps its coverage guarantee despite stopping the moment the target half-width is met."
6220 );
6221 assert_eq!(
6222 r.assumptions[2],
6223 "Input dispersions are the independent normal distributions declared in MonteCarloParams; correlations between inputs are not modeled."
6224 );
6225 assert_eq!(
6226 r.assumptions[3],
6227 "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."
6228 );
6229
6230 assert_eq!(
6232 serde_json::to_string(&McStopReason::TargetHalfWidthMet).unwrap(),
6233 "\"target_half_width_met\""
6234 );
6235 assert_eq!(
6236 serde_json::to_string(&McStopReason::MaxSamplesReached).unwrap(),
6237 "\"max_samples_reached\""
6238 );
6239 }
6240
6241 #[test]
6242 fn adaptive_rejects_nonsense_convergence() {
6243 let (inputs, wind) = seeded_test_fixture();
6244 let run = |conv: McConvergence| {
6245 run_monte_carlo_adaptive_seeded(
6246 &inputs,
6247 &wind,
6248 &loose_params(),
6249 &conv,
6250 DEFAULT_HIT_RADIUS_M,
6251 1,
6252 )
6253 .unwrap_err()
6254 };
6255
6256 for bad_width in [0.0, -0.01, f64::NAN] {
6257 let err = run(McConvergence {
6258 target_half_width: bad_width,
6259 ..Default::default()
6260 });
6261 assert!(
6262 err.contains("target_half_width"),
6263 "error must name the field, got: {err}"
6264 );
6265 }
6266
6267 let err = run(McConvergence {
6268 batch_size: 0,
6269 ..Default::default()
6270 });
6271 assert!(err.contains("batch_size"), "got: {err}");
6272
6273 let err = run(McConvergence {
6274 min_samples: 5_000,
6275 max_samples: 1_000,
6276 ..Default::default()
6277 });
6278 assert!(err.contains("max_samples"), "got: {err}");
6279 assert!(err.contains("min_samples"), "got: {err}");
6280
6281 let err = run(McConvergence {
6282 max_samples: 0,
6283 min_samples: 0,
6284 ..Default::default()
6285 });
6286 assert!(err.contains("max_samples"), "got: {err}");
6287
6288 let err = McConvergence {
6290 batch_size: 0,
6291 ..Default::default()
6292 }
6293 .validate()
6294 .unwrap_err();
6295 assert!(err.contains("batch_size"));
6296 }
6297
6298 #[test]
6299 fn wilson_companion_matches_hit_probability_and_wilson_interval() {
6300 let (inputs, wind) = seeded_test_fixture();
6301 let params = MonteCarloParams {
6302 num_simulations: 128,
6303 target_distance: Some(500.0),
6304 ..MonteCarloParams::default()
6305 };
6306 let results = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
6307 inputs,
6308 wind,
6309 params,
6310 0.01,
6311 0x1352_C0DE,
6312 )
6313 .expect("seeded legacy run");
6314
6315 for level in [
6316 ConfidenceLevel::P90,
6317 ConfidenceLevel::P95,
6318 ConfidenceLevel::P99,
6319 ] {
6320 let (p_hat, (lo, hi), n) =
6321 results.hit_probability_wilson(DEFAULT_HIT_RADIUS_M, level);
6322
6323 assert_eq!(
6326 p_hat.to_bits(),
6327 results.hit_probability(DEFAULT_HIT_RADIUS_M).to_bits()
6328 );
6329 assert_eq!(n, results.impact_positions.len() as u64);
6330
6331 let hits = results
6335 .impact_positions
6336 .iter()
6337 .filter(|p| MonteCarloResults::position_is_hit(p, DEFAULT_HIT_RADIUS_M))
6338 .count() as u64;
6339 let (want_lo, want_hi) = wilson_interval(hits, n, level);
6340 assert_eq!(lo.to_bits(), want_lo.to_bits());
6341 assert_eq!(hi.to_bits(), want_hi.to_bits());
6342 assert!(lo <= p_hat && p_hat <= hi, "interval excludes p_hat");
6343 }
6344
6345 let empty = MonteCarloResults {
6347 ranges: Vec::new(),
6348 impact_velocities: Vec::new(),
6349 impact_positions: Vec::new(),
6350 };
6351 assert_eq!(
6352 empty.hit_probability_wilson(DEFAULT_HIT_RADIUS_M, ConfidenceLevel::P95),
6353 (0.0, (0.0, 1.0), 0)
6354 );
6355 }
6356}
6357
6358#[cfg(test)]
6359mod monte_carlo_powder_curve_tests {
6360 use super::*;
6361 use rand::{rngs::StdRng, SeedableRng};
6362
6363 #[test]
6364 fn powder_curve_preserves_sampled_muzzle_velocity_dispersion() {
6365 let inputs = BallisticInputs {
6366 muzzle_velocity: 700.0,
6367 powder_temp_curve: Some(vec![(15.0, 800.0)]),
6368 powder_curve_temp_c: Some(15.0),
6369 ..BallisticInputs::default()
6370 };
6371 let params = MonteCarloParams {
6372 num_simulations: 16,
6373 velocity_std_dev: 20.0,
6374 angle_std_dev: 1e-12,
6375 bc_std_dev: 1e-12,
6376 wind_speed_std_dev: 1e-12,
6377 target_distance: Some(100.0),
6378 azimuth_std_dev: 1e-12,
6379 ..MonteCarloParams::default()
6380 };
6381
6382 let mut rng = StdRng::seed_from_u64(0x5EED_1176);
6383 let results = run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
6384 inputs,
6385 WindConditions::default(),
6386 params,
6387 0.0,
6388 &mut rng,
6389 )
6390 .expect("Monte Carlo solve");
6391 let min_velocity = results
6392 .impact_velocities
6393 .iter()
6394 .copied()
6395 .fold(f64::INFINITY, f64::min);
6396 let max_velocity = results
6397 .impact_velocities
6398 .iter()
6399 .copied()
6400 .fold(f64::NEG_INFINITY, f64::max);
6401
6402 assert!(
6403 max_velocity - min_velocity > 1.0,
6404 "20 m/s muzzle spread collapsed after curve resolution: impact-velocity span={} m/s",
6405 max_velocity - min_velocity
6406 );
6407 }
6408}
6409
6410#[cfg(test)]
6411mod monte_carlo_wind_sampling_tests {
6412 use super::*;
6413 use rand::{rngs::StdRng, SeedableRng};
6414
6415 #[test]
6416 fn wind_speed_sigma_does_not_change_seeded_direction_draws() {
6417 let base_wind = WindConditions {
6418 speed: 100.0,
6419 direction: 0.37,
6420 vertical_speed: 0.0,
6421 };
6422 let narrow_speed = MonteCarloWindSampler::new(&base_wind, 0.5, 0.2).unwrap();
6423 let wide_speed = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
6424 let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
6425 let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
6426 let mut speed_changed = false;
6427
6428 for _ in 0..32 {
6429 let narrow = narrow_speed.sample(&mut narrow_rng);
6430 let wide = wide_speed.sample(&mut wide_rng);
6431 assert!(narrow.speed > 0.0 && wide.speed > 0.0);
6432 assert_eq!(narrow.direction.to_bits(), wide.direction.to_bits());
6433 speed_changed |= narrow.speed.to_bits() != wide.speed.to_bits();
6434 }
6435 assert!(
6436 speed_changed,
6437 "different speed sigmas must still vary speed draws"
6438 );
6439 }
6440
6441 #[test]
6442 fn zero_direction_sigma_has_no_angular_jitter() {
6443 let base_wind = WindConditions {
6444 speed: 100.0,
6445 direction: 0.37,
6446 vertical_speed: 0.0,
6447 };
6448 let sampler = MonteCarloWindSampler::new(&base_wind, 4.0, 0.0).unwrap();
6449 let mut rng = StdRng::seed_from_u64(0x5EED_1223);
6450 let mut speed_changed = false;
6451
6452 for _ in 0..32 {
6453 let wind = sampler.sample(&mut rng);
6454 speed_changed |= wind.speed.to_bits() != base_wind.speed.to_bits();
6455 assert_eq!(wind.direction.to_bits(), base_wind.direction.to_bits());
6456 }
6457 assert!(speed_changed, "speed uncertainty should remain active");
6458 }
6459
6460 #[test]
6461 fn direction_sigma_controls_seeded_angular_spread_in_radians() {
6462 let base_wind = WindConditions {
6463 speed: 100.0,
6464 direction: 0.37,
6465 vertical_speed: 0.0,
6466 };
6467 let narrow = MonteCarloWindSampler::new(&base_wind, 4.0, 0.1).unwrap();
6468 let wide = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
6469 let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
6470 let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
6471 let mut nonzero_direction_draw = false;
6472
6473 for _ in 0..32 {
6474 let narrow_wind = narrow.sample(&mut narrow_rng);
6475 let wide_wind = wide.sample(&mut wide_rng);
6476 assert_eq!(narrow_wind.speed.to_bits(), wide_wind.speed.to_bits());
6477
6478 let narrow_delta = narrow_wind.direction - base_wind.direction;
6479 let wide_delta = wide_wind.direction - base_wind.direction;
6480 assert!((wide_delta - 2.0 * narrow_delta).abs() < 1e-12);
6481 nonzero_direction_draw |= narrow_delta.abs() > 1e-6;
6482 }
6483 assert!(
6484 nonzero_direction_draw,
6485 "positive radians sigma must vary direction"
6486 );
6487 }
6488
6489 #[test]
6490 fn direction_sigma_rejects_negative_or_nonfinite_values() {
6491 let base_wind = WindConditions::default();
6492 for sigma in [-0.1, f64::NAN, f64::INFINITY] {
6493 assert!(MonteCarloWindSampler::new(&base_wind, 1.0, sigma).is_err());
6494 }
6495 }
6496
6497 #[test]
6498 fn base_vertical_wind_rides_into_every_mc_sample() {
6499 use rand::SeedableRng;
6503 let base_wind = WindConditions { vertical_speed: 4.2, ..Default::default() };
6504 let sampler = MonteCarloWindSampler::new(&base_wind, 1.0, 0.2).unwrap();
6505 let mut rng = rand::rngs::StdRng::seed_from_u64(7);
6506 for _ in 0..32 {
6507 let w = sampler.sample(&mut rng);
6508 assert_eq!(w.vertical_speed, 4.2);
6509 }
6510 }
6511
6512 #[test]
6513 fn negative_speed_sample_reverses_wind_direction() {
6514 let direction = 0.25;
6515 let signed_speed = -2.5;
6516 let wind = wind_from_signed_speed_sample(signed_speed, direction, 0.0);
6517 let positive_wind = wind_from_signed_speed_sample(2.5, direction, 0.0);
6518
6519 assert_eq!(wind.speed, 2.5);
6520 assert!(
6521 (wind.direction - (direction + std::f64::consts::PI)).abs() < f64::EPSILON,
6522 "negative speed must reverse direction by pi: got {}",
6523 wind.direction
6524 );
6525 assert_eq!(positive_wind.speed, 2.5);
6526 assert_eq!(positive_wind.direction, direction);
6527
6528 let normalized_x = -wind.speed * wind.direction.cos();
6529 let normalized_z = -wind.speed * wind.direction.sin();
6530 let signed_x = -signed_speed * direction.cos();
6531 let signed_z = -signed_speed * direction.sin();
6532 assert!((normalized_x - signed_x).abs() < 1e-12);
6533 assert!((normalized_z - signed_z).abs() < 1e-12);
6534 }
6535}
6536
6537#[cfg(test)]
6538mod bc_fit_objective_tests {
6539 use super::*;
6540
6541 fn velocity_point(range_m: f64, velocity_mps: f64) -> TrajectoryPoint {
6542 TrajectoryPoint {
6543 time: 0.0,
6544 position: Vector3::new(range_m, 0.0, 0.0),
6545 velocity_magnitude: velocity_mps,
6546 kinetic_energy: 0.0,
6547 drag_coefficient: None,
6548 }
6549 }
6550
6551 #[test]
6552 fn candidate_that_misses_an_observation_has_no_score() {
6553 let trajectory = vec![velocity_point(0.0, 800.0), velocity_point(100.0, 700.0)];
6554 let observations = vec![(50.0, 750.0), (150.0, 600.0)];
6555
6556 assert!(
6557 fit_residual_sse(&trajectory, &observations, BcFitMode::Velocity, 0.0).is_none(),
6558 "a candidate that reaches only one of two observations must not compete on partial SSE"
6559 );
6560
6561 let complete_observations = vec![(50.0, 740.0), (100.0, 680.0)];
6562 assert_eq!(
6563 fit_residual_sse(
6564 &trajectory,
6565 &complete_observations,
6566 BcFitMode::Velocity,
6567 0.0,
6568 ),
6569 Some(500.0)
6570 );
6571 }
6572}
6573
6574#[cfg(test)]
6575mod cluster_bc_reference_space_tests {
6576 use super::*;
6577
6578 fn acceleration_at_1100_fps(inputs: BallisticInputs) -> Vector3<f64> {
6579 let solver = TrajectorySolver::new(
6580 inputs,
6581 WindConditions::default(),
6582 AtmosphericConditions::default(),
6583 );
6584 let position = Vector3::zeros();
6585 let velocity = Vector3::new(1100.0 / 3.28084, 0.0, 0.0);
6586 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6587 solver.calculate_acceleration(
6588 &position,
6589 &velocity,
6590 &Vector3::zeros(),
6591 (temp_c, pressure_hpa, density / 1.225),
6592 )
6593 }
6594
6595 #[test]
6596 fn solver_passes_g7_reference_model_to_cluster_classifier() {
6597 let inputs = BallisticInputs {
6598 bc_value: 0.190,
6599 bc_type: DragModel::G7,
6600 bullet_mass: 77.0 * crate::constants::GRAINS_TO_KG,
6601 bullet_diameter: 0.224 * 0.0254,
6602 use_cluster_bc: true,
6603 ..BallisticInputs::default()
6604 };
6605
6606 let solver = TrajectorySolver::new(
6607 inputs,
6608 WindConditions::default(),
6609 AtmosphericConditions::default(),
6610 );
6611 let corrected = solver.apply_cluster_bc_correction(0.190, 2800.0);
6612
6613 assert!(
6614 (corrected / 0.190 - 1.004).abs() < 1e-12,
6615 "solver selected the wrong G7 cluster multiplier: {}",
6616 corrected / 0.190
6617 );
6618 }
6619
6620 #[test]
6621 fn velocity_bc_segments_are_not_cluster_corrected_twice() {
6622 let segmented_clustered = BallisticInputs {
6623 bc_value: 0.5,
6624 bc_type: DragModel::G7,
6625 use_bc_segments: true,
6626 bc_segments_data: Some(vec![
6627 crate::BCSegmentData {
6628 velocity_min: 0.0,
6629 velocity_max: 1_600.0,
6630 bc_value: 0.4,
6631 },
6632 crate::BCSegmentData {
6633 velocity_min: 1_600.0,
6634 velocity_max: 5_000.0,
6635 bc_value: 0.45,
6636 },
6637 ]),
6638 use_cluster_bc: true,
6639 ..BallisticInputs::default()
6640 };
6641 let mut segmented_only = segmented_clustered.clone();
6642 segmented_only.use_cluster_bc = false;
6643 let mut constant_clustered = segmented_clustered.clone();
6644 constant_clustered.bc_value = 0.4;
6645 constant_clustered.bc_segments_data = None;
6646
6647 let stacked = acceleration_at_1100_fps(segmented_clustered);
6648 let segment_only = acceleration_at_1100_fps(segmented_only);
6649 let cluster_only = acceleration_at_1100_fps(constant_clustered);
6650
6651 assert!(
6652 (stacked.x - segment_only.x).abs() < 1e-12,
6653 "segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
6654 stacked.x,
6655 segment_only.x
6656 );
6657 assert!(
6658 (cluster_only.x - segment_only.x).abs() > 1e-6,
6659 "cluster correction must remain active for a constant BC"
6660 );
6661 }
6662
6663 #[test]
6664 fn mach_bc_segments_are_not_cluster_corrected_twice() {
6665 let mach_segmented_clustered = BallisticInputs {
6666 bc_value: 0.5,
6667 bc_type: DragModel::G7,
6668 use_bc_segments: false,
6669 bc_segments: Some(vec![(0.5, 0.3), (1.5, 0.5)]),
6670 use_cluster_bc: true,
6671 ..BallisticInputs::default()
6672 };
6673 let mut mach_segmented_only = mach_segmented_clustered.clone();
6674 mach_segmented_only.use_cluster_bc = false;
6675
6676 let stacked = acceleration_at_1100_fps(mach_segmented_clustered);
6677 let segment_only = acceleration_at_1100_fps(mach_segmented_only);
6678
6679 assert!(
6680 (stacked.x - segment_only.x).abs() < 1e-12,
6681 "Mach segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
6682 stacked.x,
6683 segment_only.x
6684 );
6685 }
6686}
6687
6688#[cfg(test)]
6689mod velocity_bc_flag_tests {
6690 use super::*;
6691
6692 fn acceleration_at_600_mps(inputs: BallisticInputs) -> Vector3<f64> {
6693 let solver = TrajectorySolver::new(
6694 inputs,
6695 WindConditions::default(),
6696 AtmosphericConditions::default(),
6697 );
6698 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6699 solver.calculate_acceleration(
6700 &Vector3::zeros(),
6701 &Vector3::new(600.0, 0.0, 0.0),
6702 &Vector3::zeros(),
6703 (temp_c, pressure_hpa, density / 1.225),
6704 )
6705 }
6706
6707 #[test]
6708 fn velocity_bc_data_requires_opt_in_in_trajectory_solver() {
6709 let scalar_inputs = BallisticInputs {
6710 bc_value: 0.5,
6711 bc_type: DragModel::G7,
6712 ..BallisticInputs::default()
6713 };
6714 let mut disabled_inputs = scalar_inputs.clone();
6715 disabled_inputs.bc_segments_data = Some(vec![crate::BCSegmentData {
6716 velocity_min: 0.0,
6717 velocity_max: 4_000.0,
6718 bc_value: 0.46,
6719 }]);
6720 disabled_inputs.use_bc_segments = false;
6721 let mut enabled_inputs = disabled_inputs.clone();
6722 enabled_inputs.use_bc_segments = true;
6723 let mut mach_only_inputs = scalar_inputs.clone();
6724 mach_only_inputs.bc_segments = Some(vec![(0.0, 0.4), (3.0, 0.4)]);
6725 let mut disabled_with_both = mach_only_inputs.clone();
6726 disabled_with_both.bc_segments_data = disabled_inputs.bc_segments_data.clone();
6727
6728 let scalar = acceleration_at_600_mps(scalar_inputs);
6729 let disabled = acceleration_at_600_mps(disabled_inputs);
6730 let enabled = acceleration_at_600_mps(enabled_inputs);
6731 let mach_only = acceleration_at_600_mps(mach_only_inputs);
6732 let disabled_with_both = acceleration_at_600_mps(disabled_with_both);
6733
6734 assert_eq!(
6735 disabled.x.to_bits(),
6736 scalar.x.to_bits(),
6737 "a populated velocity table must not change drag while use_bc_segments is false"
6738 );
6739 assert!(
6740 enabled.x < disabled.x - 1.0,
6741 "enabling the lower BC table must increase drag: disabled ax={} enabled ax={}",
6742 disabled.x,
6743 enabled.x
6744 );
6745 assert_eq!(
6746 disabled_with_both.x.to_bits(),
6747 mach_only.x.to_bits(),
6748 "disabling velocity data must fall through to an explicit Mach table"
6749 );
6750 }
6751}
6752
6753#[cfg(test)]
6754mod mach_bc_segment_tests {
6755 use super::*;
6756
6757 #[test]
6758 fn trajectory_solver_interpolates_explicit_mach_bc_segments() {
6759 let segmented_inputs = BallisticInputs {
6760 bc_value: 0.8,
6761 use_bc_segments: false,
6762 bc_segments: Some(vec![(1.0, 0.2), (2.0, 0.4)]),
6763 bc_segments_data: None,
6764 ..BallisticInputs::default()
6765 };
6766
6767 let mut expected_inputs = segmented_inputs.clone();
6768 expected_inputs.bc_value = 0.3;
6769 expected_inputs.bc_segments = None;
6770
6771 let atmosphere = AtmosphericConditions::default();
6772 let segmented_solver = TrajectorySolver::new(
6773 segmented_inputs,
6774 WindConditions::default(),
6775 atmosphere.clone(),
6776 );
6777 let expected_solver = TrajectorySolver::new(
6778 expected_inputs,
6779 WindConditions::default(),
6780 atmosphere,
6781 );
6782 let position = Vector3::zeros();
6783 let (density, _, temp_c, pressure_hpa) = segmented_solver.resolved_atmosphere();
6784 let (_, local_speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
6785 segmented_solver.atmosphere.altitude,
6786 segmented_solver.atmosphere.altitude,
6787 temp_c,
6788 pressure_hpa,
6789 density / 1.225,
6790 segmented_solver.atmosphere.humidity,
6791 );
6792 let velocity = Vector3::new(1.5 * local_speed_of_sound, 0.0, 0.0);
6793 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
6794
6795 let segmented_acceleration = segmented_solver.calculate_acceleration(
6796 &position,
6797 &velocity,
6798 &Vector3::zeros(),
6799 resolved_atmo,
6800 );
6801 let expected_acceleration = expected_solver.calculate_acceleration(
6802 &position,
6803 &velocity,
6804 &Vector3::zeros(),
6805 resolved_atmo,
6806 );
6807
6808 assert!(
6809 (segmented_acceleration.x - expected_acceleration.x).abs() < 1e-12,
6810 "Mach 1.5 must interpolate BC 0.3: segmented ax={} expected ax={}",
6811 segmented_acceleration.x,
6812 expected_acceleration.x
6813 );
6814 }
6815}
6816
6817#[cfg(test)]
6818mod custom_drag_table_validation_tests {
6819 use super::*;
6820
6821 #[test]
6822 fn solve_accepts_zero_bc_when_custom_table_present() {
6823 let inputs = BallisticInputs {
6824 bc_value: 0.0, bullet_mass: 0.0106,
6826 bullet_diameter: 0.00782,
6827 muzzle_velocity: 850.0,
6828 custom_drag_table: Some(crate::drag::DragTable::new(
6829 vec![0.5, 1.0, 2.0, 3.0],
6830 vec![0.23, 0.40, 0.30, 0.26],
6831 )),
6832 ..BallisticInputs::default()
6833 };
6834 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
6835 assert!(solver.solve().is_ok());
6837 }
6838
6839 #[test]
6840 fn solve_still_requires_bc_without_table() {
6841 let inputs = BallisticInputs {
6842 bc_value: 0.0,
6843 bullet_mass: 0.0106,
6844 bullet_diameter: 0.00782,
6845 muzzle_velocity: 850.0,
6846 ..BallisticInputs::default()
6847 };
6848 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
6849 assert!(solver.solve().is_err());
6850 }
6851}
6852
6853#[cfg(test)]
6855mod cd_scale_tests {
6856 use super::*;
6857
6858 fn deck() -> crate::drag::DragTable {
6859 crate::drag::DragTable::new(vec![0.5, 1.0, 2.0, 3.0], vec![0.23, 0.40, 0.30, 0.26])
6860 }
6861
6862 fn deck_inputs(cd_scale: f64) -> BallisticInputs {
6863 BallisticInputs {
6864 bullet_mass: 0.0106,
6865 bullet_diameter: 0.00782,
6866 muzzle_velocity: 850.0,
6867 custom_drag_table: Some(deck()),
6868 cd_scale,
6869 ..BallisticInputs::default()
6870 }
6871 }
6872
6873 #[test]
6874 fn default_cd_scale_is_one() {
6875 assert_eq!(BallisticInputs::default().cd_scale, 1.0);
6876 }
6877
6878 #[test]
6881 fn cd_scale_absent_is_byte_identical_to_explicit_one() {
6882 let omitted = BallisticInputs {
6883 bullet_mass: 0.0106,
6884 bullet_diameter: 0.00782,
6885 muzzle_velocity: 850.0,
6886 custom_drag_table: Some(deck()),
6887 ..BallisticInputs::default()
6888 };
6889 let explicit = BallisticInputs {
6890 cd_scale: 1.0,
6891 ..omitted.clone()
6892 };
6893
6894 let solver_omitted =
6895 TrajectorySolver::new(omitted, WindConditions::default(), AtmosphericConditions::default());
6896 let solver_explicit =
6897 TrajectorySolver::new(explicit, WindConditions::default(), AtmosphericConditions::default());
6898
6899 let cd_omitted = solver_omitted.calculate_drag_coefficient(700.0, 340.0);
6900 let cd_explicit = solver_explicit.calculate_drag_coefficient(700.0, 340.0);
6901 assert_eq!(
6902 cd_omitted.to_bits(),
6903 cd_explicit.to_bits(),
6904 "default cd_scale must be bit-identical to an explicit 1.0"
6905 );
6906
6907 let result = solver_omitted.solve();
6910 assert!(result.is_ok(), "existing custom-deck solves must pass unchanged");
6911 }
6912
6913 #[test]
6915 fn cd_scale_multiplies_the_interpolated_cd_exactly() {
6916 let velocity = 700.0;
6917 let speed_of_sound = 340.0;
6918 let mach = velocity / speed_of_sound;
6919 let expected_unscaled = deck().interpolate(mach);
6920
6921 for &scale in &[0.90, 1.0, 1.10, 1.5] {
6922 let solver = TrajectorySolver::new(
6923 deck_inputs(scale),
6924 WindConditions::default(),
6925 AtmosphericConditions::default(),
6926 );
6927 let cd = solver.calculate_drag_coefficient(velocity, speed_of_sound);
6928 assert!(
6929 (cd - expected_unscaled * scale).abs() < 1e-12,
6930 "scale={scale}: cd={cd} expected={}",
6931 expected_unscaled * scale
6932 );
6933 }
6934 }
6935
6936 #[test]
6940 fn cd_scale_direction_on_cli_api_solver() {
6941 let solve = |scale: f64| {
6942 TrajectorySolver::new(
6943 deck_inputs(scale),
6944 WindConditions::default(),
6945 AtmosphericConditions::default(),
6946 )
6947 .solve()
6948 .expect("custom-deck solve should succeed")
6949 };
6950
6951 let baseline = solve(1.0);
6952 let scaled_up = solve(1.10);
6953 let scaled_down = solve(0.90);
6954
6955 assert!(
6956 scaled_up.impact_velocity < baseline.impact_velocity,
6957 "cd_scale=1.10 must increase drag -> lower impact velocity: base={} up={}",
6958 baseline.impact_velocity,
6959 scaled_up.impact_velocity
6960 );
6961 assert!(
6962 scaled_down.impact_velocity > baseline.impact_velocity,
6963 "cd_scale=0.90 must decrease drag -> higher impact velocity: base={} down={}",
6964 baseline.impact_velocity,
6965 scaled_down.impact_velocity
6966 );
6967 }
6968
6969 #[test]
6971 fn validate_for_solve_rejects_invalid_cd_scale() {
6972 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
6973 let solver = TrajectorySolver::new(
6974 deck_inputs(bad),
6975 WindConditions::default(),
6976 AtmosphericConditions::default(),
6977 );
6978 assert!(
6979 solver.solve().is_err(),
6980 "cd_scale={bad} must be rejected by validate_for_solve"
6981 );
6982 }
6983 }
6984
6985 #[test]
6992 fn validate_for_solve_rejects_invalid_cd_scale_without_a_custom_drag_table() {
6993 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
6994 let inputs = BallisticInputs {
6995 bc_value: 0.5,
6996 bc_type: crate::DragModel::G1,
6997 bullet_mass: 0.0106,
6998 bullet_diameter: 0.00782,
6999 muzzle_velocity: 850.0,
7000 cd_scale: bad,
7001 ..BallisticInputs::default()
7002 };
7003 assert!(inputs.custom_drag_table.is_none(), "precondition: no custom deck");
7004 let solver = TrajectorySolver::new(
7005 inputs,
7006 WindConditions::default(),
7007 AtmosphericConditions::default(),
7008 );
7009 assert!(
7010 solver.solve().is_err(),
7011 "cd_scale={bad} must be rejected by validate_for_solve even without a custom \
7012 drag table"
7013 );
7014 }
7015 }
7016
7017 #[test]
7021 fn cd_scale_is_inert_without_a_custom_drag_table() {
7022 let make = |cd_scale: f64| BallisticInputs {
7023 bc_value: 0.5,
7024 bc_type: crate::DragModel::G1,
7025 bullet_mass: 0.0106,
7026 bullet_diameter: 0.00782,
7027 muzzle_velocity: 850.0,
7028 cd_scale,
7029 ..BallisticInputs::default()
7030 };
7031 let solver_neutral = TrajectorySolver::new(
7032 make(1.0),
7033 WindConditions::default(),
7034 AtmosphericConditions::default(),
7035 );
7036 let solver_far = TrajectorySolver::new(
7037 make(1.5),
7038 WindConditions::default(),
7039 AtmosphericConditions::default(),
7040 );
7041 let cd_neutral = solver_neutral.calculate_drag_coefficient(700.0, 340.0);
7042 let cd_far = solver_far.calculate_drag_coefficient(700.0, 340.0);
7043 assert_eq!(
7044 cd_neutral.to_bits(),
7045 cd_far.to_bits(),
7046 "cd_scale must not affect the G-model/BC drag path"
7047 );
7048 }
7049
7050 #[test]
7053 fn cd_scale_shifts_all_three_solver_paths_in_the_same_direction() {
7054 let cli_solve = |scale: f64| {
7056 TrajectorySolver::new(
7057 deck_inputs(scale),
7058 WindConditions::default(),
7059 AtmosphericConditions::default(),
7060 )
7061 .solve()
7062 .expect("cli_api custom-deck solve should succeed")
7063 };
7064 let cli_baseline = cli_solve(1.0);
7065 let cli_scaled = cli_solve(1.10);
7066 assert!(
7067 cli_scaled.impact_velocity < cli_baseline.impact_velocity,
7068 "cli_api: cd_scale=1.10 must lower impact velocity"
7069 );
7070
7071 let derivatives_accel_x = |scale: f64| {
7073 let inputs = deck_inputs(scale);
7074 crate::derivatives::compute_derivatives(
7075 nalgebra::Vector3::zeros(),
7076 nalgebra::Vector3::new(700.0, 0.0, 0.0),
7077 &inputs,
7078 nalgebra::Vector3::zeros(),
7079 (1.225, 340.0, 0.0, 0.0),
7080 inputs.bc_value,
7081 None,
7082 0.0,
7083 None,
7084 )[3]
7085 };
7086 let deriv_baseline = derivatives_accel_x(1.0);
7087 let deriv_scaled = derivatives_accel_x(1.10);
7088 assert!(
7089 deriv_scaled < deriv_baseline,
7090 "derivatives: cd_scale=1.10 must make x-acceleration more negative (more drag): \
7091 base={deriv_baseline} scaled={deriv_scaled}"
7092 );
7093
7094 let fast_final_speed = |scale: f64| {
7096 let inputs = deck_inputs(scale);
7097 let wind_sock = crate::wind::WindSock::new(vec![]);
7098 let params = crate::fast_trajectory::FastIntegrationParams {
7099 horiz: 500.0,
7100 vert: 0.0,
7101 initial_state: [0.0, 0.0, 0.0, 850.0, 0.0, 0.0],
7102 t_span: (0.0, 5.0),
7103 atmo_params: (0.0, 15.0, 1013.25, 1.0),
7104 atmo_sock: None,
7105 };
7106 let solution = crate::fast_trajectory::fast_integrate(&inputs, &wind_sock, params);
7107 assert!(solution.success, "fast_integrate must succeed for scale={scale}");
7108 let last = solution.t.len() - 1;
7109 let (vx, vy, vz) = (
7110 solution.y[3][last],
7111 solution.y[4][last],
7112 solution.y[5][last],
7113 );
7114 (vx * vx + vy * vy + vz * vz).sqrt()
7115 };
7116 let fast_baseline = fast_final_speed(1.0);
7117 let fast_scaled = fast_final_speed(1.10);
7118 assert!(
7119 fast_scaled < fast_baseline,
7120 "fast_trajectory: cd_scale=1.10 must lower final speed: base={fast_baseline} scaled={fast_scaled}"
7121 );
7122 }
7123}
7124
7125#[cfg(test)]
7126mod humid_local_mach_tests {
7127 use super::*;
7128
7129 fn solver_with_station_humidity(humidity_percent: f64) -> TrajectorySolver {
7130 let inputs = BallisticInputs {
7131 custom_drag_table: Some(crate::drag::DragTable::new(vec![0.5, 1.5], vec![0.1, 1.1])),
7132 ..BallisticInputs::default()
7133 };
7134 TrajectorySolver::new(
7135 inputs,
7136 WindConditions::default(),
7137 AtmosphericConditions {
7138 temperature: 30.0,
7139 pressure: 1013.25,
7140 humidity: humidity_percent,
7141 altitude: 0.0,
7142 },
7143 )
7144 }
7145
7146 fn acceleration(solver: &TrajectorySolver, base_ratio: f64) -> Vector3<f64> {
7147 solver.calculate_acceleration(
7148 &Vector3::zeros(),
7149 &Vector3::new(350.0, 0.0, 0.0),
7150 &Vector3::zeros(),
7151 (30.0, 1013.25, base_ratio),
7152 )
7153 }
7154
7155 #[test]
7156 fn local_mach_uses_station_humidity_when_density_is_held_constant() {
7157 let dry = acceleration(&solver_with_station_humidity(0.0), 1.0);
7158 let humid = acceleration(&solver_with_station_humidity(100.0), 1.0);
7159
7160 assert!(
7161 humid.x > dry.x,
7162 "humid sound speed should lower Mach and drag on the rising test curve: dry ax={} humid ax={}",
7163 dry.x,
7164 humid.x
7165 );
7166 }
7167
7168 #[test]
7169 fn active_atmosphere_zone_uses_zone_humidity_instead_of_station_humidity() {
7170 let zone_humidity = 80.0;
7171 let zone_ratio =
7172 crate::atmosphere::calculate_air_density_cimp(30.0, 1013.25, zone_humidity) / 1.225;
7173 let station_solver = solver_with_station_humidity(zone_humidity);
7174 let mut zoned_solver = solver_with_station_humidity(0.0);
7175 zoned_solver.set_atmo_segments(vec![(30.0, 1013.25, zone_humidity, 1_000.0)]);
7176
7177 let station = acceleration(&station_solver, zone_ratio);
7178 let zoned = acceleration(&zoned_solver, zone_ratio);
7179
7180 assert!(
7181 (zoned - station).norm() < 1e-12,
7182 "active zone T/P/RH should override the station atmosphere: station={station:?} zoned={zoned:?}"
7183 );
7184 }
7185}
7186
7187#[cfg(test)]
7188mod inclined_atmosphere_frame_tests {
7189 use super::*;
7190
7191 fn expected_shot_frame_vector(level: Vector3<f64>, angle: f64) -> Vector3<f64> {
7192 let (sin_angle, cos_angle) = angle.sin_cos();
7193 Vector3::new(
7194 level.x * cos_angle + level.y * sin_angle,
7195 -level.x * sin_angle + level.y * cos_angle,
7196 level.z,
7197 )
7198 }
7199
7200 #[test]
7201 fn inclined_positions_at_same_world_altitude_have_same_solver_acceleration() {
7202 let angle = std::f64::consts::FRAC_PI_6;
7203 let inputs = BallisticInputs {
7204 shooting_angle: angle,
7205 ..BallisticInputs::default()
7206 };
7207 let atmosphere = AtmosphericConditions {
7208 altitude: 100.0,
7209 ..AtmosphericConditions::default()
7210 };
7211 let solver = TrajectorySolver::new(inputs, WindConditions::default(), atmosphere);
7212 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7213 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
7214 let velocity = Vector3::new(600.0, 0.0, 0.0);
7215 let along_slant = Vector3::new(1_000.0, 0.0, 0.0);
7216 let across_slant = Vector3::new(0.0, 500.0 / angle.cos(), 0.0);
7217
7218 let a = solver.calculate_acceleration(
7219 &along_slant,
7220 &velocity,
7221 &Vector3::zeros(),
7222 resolved_atmo,
7223 );
7224 let b = solver.calculate_acceleration(
7225 &across_slant,
7226 &velocity,
7227 &Vector3::zeros(),
7228 resolved_atmo,
7229 );
7230
7231 assert!(
7232 (a - b).norm() < 1e-10,
7233 "solver acceleration differs at equal world altitude: {a:?} vs {b:?}"
7234 );
7235 }
7236
7237 #[test]
7238 fn inclined_headwind_is_rotated_into_solver_frame() {
7239 let angle = std::f64::consts::FRAC_PI_6;
7240 let inputs = BallisticInputs {
7241 shooting_angle: angle,
7242 ..BallisticInputs::default()
7243 };
7244 let solver = TrajectorySolver::new(
7245 inputs,
7246 WindConditions::default(),
7247 AtmosphericConditions::default(),
7248 );
7249 let level_headwind = Vector3::new(-100.0, 0.0, 0.0);
7250 let velocity = expected_shot_frame_vector(level_headwind, angle);
7251 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7252 let actual = solver.calculate_acceleration(
7253 &Vector3::zeros(),
7254 &velocity,
7255 &level_headwind,
7256 (temp_c, pressure_hpa, density / 1.225),
7257 );
7258
7259 assert!(
7260 (actual - solver.gravity_acceleration()).norm() < 1e-12,
7261 "co-moving horizontal wind must leave only shot-frame gravity: {actual:?}"
7262 );
7263 }
7264
7265 #[test]
7266 fn inclined_coriolis_is_rotated_into_solver_frame() {
7267 let angle = std::f64::consts::FRAC_PI_6;
7268 let latitude_deg = 45.0_f64;
7269 let shot_azimuth = 0.4_f64;
7270 let velocity = Vector3::new(600.0, 20.0, 5.0);
7271 let base_inputs = BallisticInputs {
7272 shooting_angle: angle,
7273 latitude: Some(latitude_deg),
7274 shot_azimuth,
7275 ..BallisticInputs::default()
7276 };
7277 let acceleration = |enable_coriolis| {
7278 let mut inputs = base_inputs.clone();
7279 inputs.enable_coriolis = enable_coriolis;
7280 let solver = TrajectorySolver::new(
7281 inputs,
7282 WindConditions::default(),
7283 AtmosphericConditions::default(),
7284 );
7285 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7286 solver.calculate_acceleration(
7287 &Vector3::zeros(),
7288 &velocity,
7289 &Vector3::zeros(),
7290 (temp_c, pressure_hpa, density / 1.225),
7291 )
7292 };
7293
7294 let omega_earth = 7.2921159e-5_f64;
7295 let latitude = latitude_deg.to_radians();
7296 let level_omega = Vector3::new(
7297 omega_earth * latitude.cos() * shot_azimuth.cos(),
7298 omega_earth * latitude.sin(),
7299 -omega_earth * latitude.cos() * shot_azimuth.sin(),
7300 );
7301 let expected = -2.0 * expected_shot_frame_vector(level_omega, angle).cross(&velocity);
7302 let actual = acceleration(true) - acceleration(false);
7303
7304 assert!(
7305 (actual - expected).norm() < 1e-12,
7306 "inclined Coriolis mismatch: actual={actual:?}, expected={expected:?}"
7307 );
7308 }
7309}
7310
7311#[cfg(test)]
7312mod terminal_range_interpolation_tests {
7313 use super::*;
7314
7315 #[test]
7316 fn terminal_finalizer_selects_the_earliest_crossed_boundary() {
7317 let inputs = BallisticInputs {
7318 ground_threshold: 0.0,
7319 ..BallisticInputs::default()
7320 };
7321 let mut solver = TrajectorySolver::new(
7322 inputs,
7323 WindConditions::default(),
7324 AtmosphericConditions::default(),
7325 );
7326 solver.set_max_range(120.0);
7327
7328 let previous_speed = 700.0;
7329 let mut points = vec![TrajectoryPoint {
7330 time: 99.0,
7331 position: Vector3::new(90.0, 1.0, -1.0),
7332 velocity_magnitude: previous_speed,
7333 kinetic_energy: 0.5 * solver.inputs.bullet_mass * previous_speed.powi(2),
7334 drag_coefficient: None,
7335 }];
7336 let mut max_height = 1.0;
7337 let termination = solver
7338 .append_terminal_endpoint(
7339 &mut points,
7340 Vector3::new(130.0, -3.0, 3.0),
7341 Vector3::new(600.0, 0.0, 0.0),
7342 101.0,
7343 &mut max_height,
7344 )
7345 .expect("the final step brackets supported boundaries");
7346
7347 assert_eq!(termination, TrajectoryTermination::GroundThreshold);
7348 assert_eq!(points.len(), 2);
7349 let terminal = points.last().expect("terminal point");
7350 assert_eq!(terminal.time, 99.5);
7351 assert_eq!(terminal.position, Vector3::new(100.0, 0.0, 0.0));
7352 assert_eq!(terminal.velocity_magnitude, 675.0);
7353 assert_eq!(
7354 terminal.kinetic_energy,
7355 0.5 * solver.inputs.bullet_mass * 675.0_f64.powi(2)
7356 );
7357
7358 solver.set_max_range(100.0);
7360 let mut tied_points = vec![points[0].clone()];
7361 assert_eq!(
7362 solver
7363 .append_terminal_endpoint(
7364 &mut tied_points,
7365 Vector3::new(130.0, -3.0, 3.0),
7366 Vector3::new(600.0, 0.0, 0.0),
7367 101.0,
7368 &mut max_height,
7369 )
7370 .expect("tied boundaries remain a valid terminal"),
7371 TrajectoryTermination::GroundThreshold
7372 );
7373 }
7374
7375 #[test]
7376 fn sub_ulp_terminal_crossing_replaces_instead_of_duplicating_range() {
7377 let ground_threshold = f64::from_bits(1.0_f64.to_bits() - 1);
7378 let inputs = BallisticInputs {
7379 ground_threshold,
7380 ..BallisticInputs::default()
7381 };
7382 let mut solver = TrajectorySolver::new(
7383 inputs,
7384 WindConditions::default(),
7385 AtmosphericConditions::default(),
7386 );
7387 solver.set_max_range(1_000.0);
7388
7389 let speed = 700.0;
7390 let mut points = vec![TrajectoryPoint {
7391 time: 0.0,
7392 position: Vector3::new(100.0, 1.0, 0.0),
7393 velocity_magnitude: speed,
7394 kinetic_energy: 0.5 * solver.inputs.bullet_mass * speed.powi(2),
7395 drag_coefficient: None,
7396 }];
7397 let mut max_height = 1.0;
7398 let termination = solver
7399 .append_terminal_endpoint(
7400 &mut points,
7401 Vector3::new(101.0, 0.0, 0.0),
7402 Vector3::new(699.0, 0.0, 0.0),
7403 1.0,
7404 &mut max_height,
7405 )
7406 .expect("sub-ULP ground crossing remains representable as one terminal state");
7407
7408 assert_eq!(termination, TrajectoryTermination::GroundThreshold);
7409 assert_eq!(points.len(), 1);
7410 assert_eq!(points[0].position.x, 100.0);
7411 assert_eq!(points[0].position.y.to_bits(), ground_threshold.to_bits());
7412 assert!(points[0].time > 0.0);
7413 }
7414
7415 #[test]
7416 fn every_solver_appends_an_exact_max_range_endpoint() {
7417 let target_range = 0.1;
7418 let modes = [
7419 ("Euler", false, false),
7420 ("RK4", true, false),
7421 ("RK45", true, true),
7422 ];
7423
7424 for (name, use_rk4, use_adaptive_rk45) in modes {
7425 let inputs = BallisticInputs {
7426 use_rk4,
7427 use_adaptive_rk45,
7428 ground_threshold: f64::NEG_INFINITY,
7429 enable_trajectory_sampling: true,
7430 sample_interval: target_range,
7431 ..BallisticInputs::default()
7432 };
7433 let mut solver = TrajectorySolver::new(
7434 inputs,
7435 WindConditions::default(),
7436 AtmosphericConditions::default(),
7437 );
7438 solver.set_max_range(target_range);
7439
7440 let result = solver.solve().expect("short-range solve should succeed");
7441 let terminal = result.points.last().expect("terminal point is missing");
7442 let muzzle = result.points.first().expect("muzzle point is missing");
7443
7444 assert_eq!(result.termination, TrajectoryTermination::MaxRange);
7445 assert_eq!(
7446 terminal.position.x.to_bits(),
7447 target_range.to_bits(),
7448 "{name} did not terminate exactly at max_range"
7449 );
7450 assert_eq!(result.max_range.to_bits(), target_range.to_bits());
7451 assert!(
7452 result.time_of_flight > 0.0 && result.time_of_flight < solver.time_step,
7453 "{name} terminal time was not interpolated within the crossing step: {}",
7454 result.time_of_flight
7455 );
7456 assert_eq!(result.time_of_flight.to_bits(), terminal.time.to_bits());
7457 assert_eq!(
7458 result.impact_velocity.to_bits(),
7459 terminal.velocity_magnitude.to_bits()
7460 );
7461 assert_eq!(
7462 result.impact_energy.to_bits(),
7463 terminal.kinetic_energy.to_bits()
7464 );
7465 let expected_energy = 0.5 * solver.inputs.bullet_mass * result.impact_velocity.powi(2);
7466 assert!((result.impact_energy - expected_energy).abs() < 1e-12);
7467 assert!(terminal.velocity_magnitude < muzzle.velocity_magnitude);
7468 assert!(terminal.kinetic_energy < muzzle.kinetic_energy);
7469
7470 let terminal_sample = result
7471 .sampled_points
7472 .as_ref()
7473 .and_then(|samples| samples.last())
7474 .expect("terminal trajectory sample is missing");
7475 assert_eq!(
7476 terminal_sample.distance_m.to_bits(),
7477 target_range.to_bits(),
7478 "{name} sampling did not include max_range"
7479 );
7480 assert_eq!(
7481 terminal_sample.time_s.to_bits(),
7482 result.time_of_flight.to_bits()
7483 );
7484 assert_eq!(
7485 terminal_sample.velocity_mps.to_bits(),
7486 result.impact_velocity.to_bits()
7487 );
7488 assert!((terminal_sample.energy_j - result.impact_energy).abs() < 1e-12);
7489 }
7490 }
7491}
7492
7493#[cfg(test)]
7494mod precession_inertia_wiring_tests {
7495 use super::*;
7496
7497 #[test]
7498 fn solver_uses_projectile_specific_moments_of_inertia() {
7499 let mass_kg = 55.0 * crate::constants::GRAINS_TO_KG;
7500 let caliber_m = 0.224 * 0.0254;
7501 let length_m = 0.75 * 0.0254;
7502 let inputs = BallisticInputs {
7503 bullet_mass: mass_kg,
7504 bullet_diameter: caliber_m,
7505 bullet_length: length_m,
7506 muzzle_velocity: 800.0,
7507 twist_rate: 7.0,
7508 enable_precession_nutation: true,
7509 use_rk4: false,
7510 use_adaptive_rk45: false,
7511 ..BallisticInputs::default()
7512 };
7513 let mut solver = TrajectorySolver::new(
7514 inputs,
7515 WindConditions::default(),
7516 AtmosphericConditions::default(),
7517 );
7518 solver.set_max_range(0.1);
7519
7520 let (air_density, speed_of_sound, _, _) = solver.resolved_atmosphere();
7521 let velocity_mps = solver.inputs.muzzle_velocity;
7522 let velocity_fps = velocity_mps * 3.28084;
7523 let twist_rate_ft = solver.inputs.twist_rate / 12.0;
7524 let spin_rate_rad_s = (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI;
7525 let initial_state = AngularState {
7526 pitch_angle: 0.001,
7527 yaw_angle: 0.001,
7528 pitch_rate: 0.0,
7529 yaw_rate: 0.0,
7530 precession_angle: 0.0,
7531 nutation_phase: 0.0,
7532 };
7533 let params = PrecessionNutationParams {
7534 mass_kg,
7535 caliber_m,
7536 length_m,
7537 spin_rate_rad_s,
7538 spin_inertia: crate::spin_decay::calculate_moment_of_inertia(
7539 mass_kg, caliber_m, length_m, "ogive",
7540 ),
7541 transverse_inertia: crate::pitch_damping::calculate_transverse_moment_of_inertia(
7542 mass_kg, caliber_m, length_m, "ogive",
7543 ),
7544 velocity_mps,
7545 air_density_kg_m3: air_density,
7546 mach: velocity_mps / speed_of_sound,
7547 pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
7548 nutation_damping_factor: 0.05,
7549 };
7550 let expected = calculate_combined_angular_motion(
7551 ¶ms,
7552 &initial_state,
7553 0.0,
7554 solver.time_step,
7555 0.001,
7556 );
7557 let actual = solver
7558 .solve()
7559 .expect("one-step solve should succeed")
7560 .angular_state
7561 .expect("precession/nutation was enabled");
7562
7563 assert!(
7564 (actual.precession_angle - expected.precession_angle).abs() < 1e-15,
7565 "precession phase used the wrong inertia: actual={}, expected={}",
7566 actual.precession_angle,
7567 expected.precession_angle
7568 );
7569 assert!(
7570 (actual.nutation_phase - expected.nutation_phase).abs() < 1e-15,
7571 "nutation phase used the wrong inertia: actual={}, expected={}",
7572 actual.nutation_phase,
7573 expected.nutation_phase
7574 );
7575 }
7576}
7577
7578#[cfg(test)]
7579mod form_factor_drag_tests {
7580 use super::*;
7581
7582 fn acceleration_with_form_factor_flag(enabled: bool) -> Vector3<f64> {
7583 let inputs = BallisticInputs {
7584 bc_value: 0.462,
7585 bc_type: DragModel::G1,
7586 bullet_model: Some("168gr SMK Match".to_string()),
7587 use_form_factor: enabled,
7588 ..BallisticInputs::default()
7589 };
7590 let solver = TrajectorySolver::new(
7591 inputs,
7592 WindConditions::default(),
7593 AtmosphericConditions::default(),
7594 );
7595 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7596 solver.calculate_acceleration(
7597 &Vector3::zeros(),
7598 &Vector3::new(600.0, 0.0, 0.0),
7599 &Vector3::zeros(),
7600 (temp_c, pressure_hpa, density / 1.225),
7601 )
7602 }
7603
7604 #[test]
7605 fn measured_bc_drag_does_not_apply_name_based_form_factor_again() {
7606 let baseline = acceleration_with_form_factor_flag(false);
7607 let flagged = acceleration_with_form_factor_flag(true);
7608
7609 assert!(
7610 (flagged - baseline).norm() < 1e-12,
7611 "published BC already encodes form factor: baseline={baseline:?} flagged={flagged:?}"
7612 );
7613 }
7614}
7615
7616#[cfg(test)]
7617mod rk45_adaptivity_tests {
7618 use super::*;
7619
7620 #[test]
7621 fn cli_rk45_error_norm_scales_components_independently() {
7622 let position = Vector3::new(1.0e9, 0.0, 0.0);
7623 let velocity = Vector3::new(800.0, 0.0, 0.0);
7624 let fifth_position = position;
7625 let fifth_velocity = velocity;
7626 let fourth_position = position;
7627 let fourth_velocity = Vector3::new(800.0, 1.0e-3, 0.0);
7628
7629 let error = cli_rk45_error_norm(
7630 &position,
7631 &velocity,
7632 &fifth_position,
7633 &fifth_velocity,
7634 &fourth_position,
7635 &fourth_velocity,
7636 );
7637 let expected = 1.0e-3 / 6.0_f64.sqrt();
7638
7639 assert!(
7640 (error - expected).abs() <= 1e-15,
7641 "large downrange position masked a velocity-component error: {error}"
7642 );
7643 }
7644
7645 fn discontinuous_wind_solver() -> TrajectorySolver {
7646 let inputs = BallisticInputs::default();
7647 let mut solver = TrajectorySolver::new(
7648 inputs,
7649 WindConditions::default(),
7650 AtmosphericConditions::default(),
7651 );
7652 solver.set_wind_segments(vec![
7653 crate::wind::WindSegment::new(0.0, 90.0, 4.0),
7654 crate::wind::WindSegment::new(1_000.0, 90.0, 10_000.0),
7655 ]);
7656 solver
7657 }
7658
7659 #[test]
7660 fn rk45_retries_discontinuous_trial_before_advancing() {
7661 let solver = discontinuous_wind_solver();
7662 let position = Vector3::new(0.0, solver.inputs.muzzle_height, 0.0);
7663 let velocity = Vector3::new(solver.inputs.muzzle_velocity, 0.0, 0.0);
7664 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7665 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
7666 let dt = 0.01;
7667
7668 let rejected_trial = solver.rk45_step(
7669 &position,
7670 &velocity,
7671 dt,
7672 &Vector3::zeros(),
7673 RK45_TOLERANCE,
7674 resolved_atmo,
7675 );
7676 assert!(
7677 rejected_trial.error > RK45_TOLERANCE,
7678 "discontinuous full step must exceed tolerance, got {}",
7679 rejected_trial.error
7680 );
7681
7682 let accepted = solver.adaptive_rk45_step(
7683 &position,
7684 &velocity,
7685 dt,
7686 &Vector3::zeros(),
7687 resolved_atmo,
7688 );
7689 assert!(accepted.used_dt < dt, "oversized trial was not retried");
7690 assert!(
7691 accepted.error <= RK45_TOLERANCE || accepted.used_dt <= RK45_MIN_DT,
7692 "accepted error {} exceeds tolerance at dt {}",
7693 accepted.error,
7694 accepted.used_dt
7695 );
7696
7697 let accepted_trial = solver.rk45_step(
7698 &position,
7699 &velocity,
7700 accepted.used_dt,
7701 &Vector3::zeros(),
7702 RK45_TOLERANCE,
7703 resolved_atmo,
7704 );
7705 assert_eq!(accepted.position, accepted_trial.position);
7706 assert_eq!(accepted.velocity, accepted_trial.velocity);
7707 assert!((RK45_MIN_DT..=RK45_MAX_DT).contains(&accepted.next_dt));
7708 }
7709}
7710
7711#[cfg(test)]
7712mod ground_termination_tests {
7713 use super::*;
7714 use crate::trajectory_observation::TrajectoryObservationFlag;
7715
7716 #[test]
7717 fn every_solver_reports_one_exact_early_ground_endpoint() {
7718 for (name, use_rk4, use_adaptive_rk45) in [
7719 ("Euler", false, false),
7720 ("RK4", true, false),
7721 ("RK45", true, true),
7722 ] {
7723 let inputs = BallisticInputs {
7724 muzzle_height: 1.0,
7725 muzzle_angle: -0.2,
7726 ground_threshold: 0.0,
7727 use_rk4,
7728 use_adaptive_rk45,
7729 ..BallisticInputs::default()
7730 };
7731 let mut solver = TrajectorySolver::new(
7732 inputs,
7733 WindConditions::default(),
7734 AtmosphericConditions::default(),
7735 );
7736 solver.set_max_range(1_000.0);
7737
7738 let result = solver.solve().expect("early-ground solve should succeed");
7739 let terminal = result.points.last().expect("terminal point is missing");
7740
7741 assert_eq!(result.termination, TrajectoryTermination::GroundThreshold);
7742 assert_eq!(terminal.position.y.to_bits(), 0.0_f64.to_bits());
7743 assert!(
7744 terminal.position.x < 1_000.0,
7745 "{name} incorrectly reached max range"
7746 );
7747 assert_eq!(result.max_range.to_bits(), terminal.position.x.to_bits());
7748 assert_eq!(
7749 result
7750 .points
7751 .iter()
7752 .filter(|point| point.position.y == 0.0)
7753 .count(),
7754 1,
7755 "{name} did not retain exactly one ground endpoint"
7756 );
7757
7758 let observations = result
7759 .sample_observations(1.0, 100)
7760 .expect("checked early-ground sampling should succeed");
7761 assert!(observations[..observations.len() - 1]
7762 .iter()
7763 .all(|observation| observation.distance_m < terminal.position.x));
7764 let terminal_observation = observations.last().expect("terminal observation");
7765 assert_eq!(
7766 terminal_observation.distance_m.to_bits(),
7767 terminal.position.x.to_bits()
7768 );
7769 assert!(terminal_observation
7770 .flags
7771 .contains(&TrajectoryObservationFlag::Terminal));
7772 assert!(terminal_observation
7773 .flags
7774 .contains(&TrajectoryObservationFlag::GroundThreshold));
7775 assert_eq!(
7776 observations
7777 .iter()
7778 .filter(|observation| observation
7779 .flags
7780 .contains(&TrajectoryObservationFlag::Terminal))
7781 .count(),
7782 1,
7783 "{name} repeated the terminal observation"
7784 );
7785 }
7786 }
7787
7788 #[test]
7793 fn rk4_and_rk45_descend_to_ground_threshold() {
7794 for adaptive in [false, true] {
7795 let inputs = BallisticInputs {
7796 muzzle_angle: 0.1, use_rk4: true,
7798 use_adaptive_rk45: adaptive,
7799 ..BallisticInputs::default()
7800 };
7801 assert_eq!(
7802 inputs.ground_threshold, -100.0,
7803 "default ground_threshold is -100 m"
7804 );
7805
7806 let mut solver = TrajectorySolver::new(
7807 inputs,
7808 WindConditions::default(),
7809 AtmosphericConditions::default(),
7810 );
7811 solver.set_max_range(1.0e7);
7813
7814 let result = solver.solve().expect("solve should succeed");
7815 let final_y = result
7816 .points
7817 .last()
7818 .expect("trajectory has points")
7819 .position
7820 .y;
7821 assert!(
7822 final_y < -1.0,
7823 "adaptive_rk45={adaptive}: final y = {final_y} m; a lofted shot should descend \
7824 past launch level toward the ground_threshold floor, not stop at y = 0"
7825 );
7826 }
7827 }
7828}
7829
7830#[cfg(test)]
7831mod magnus_stability_tests {
7832 use super::*;
7833
7834 #[test]
7835 fn yaw_of_repose_magnus_force_is_vertical_and_twist_signed() {
7836 let acceleration = |enable_magnus, is_twist_right| {
7837 let inputs = BallisticInputs {
7838 muzzle_velocity: 822.96,
7839 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
7840 bullet_diameter: 0.308 * 0.0254,
7841 bullet_length: 1.215 * 0.0254,
7842 twist_rate: 10.0,
7843 is_twist_right,
7844 enable_magnus,
7845 ..BallisticInputs::default()
7846 };
7847 let solver = TrajectorySolver::new(
7848 inputs,
7849 WindConditions::default(),
7850 AtmosphericConditions::default(),
7851 );
7852 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7853 solver.calculate_acceleration(
7854 &Vector3::zeros(),
7855 &Vector3::new(822.96, 0.0, 0.0),
7856 &Vector3::zeros(),
7857 (temp_c, pressure_hpa, density / 1.225),
7858 )
7859 };
7860
7861 let baseline = acceleration(false, true);
7862 let right_twist = acceleration(true, true) - baseline;
7863 let left_twist = acceleration(true, false) - baseline;
7864
7865 assert!(
7866 right_twist.y < 0.0,
7867 "right-hand Magnus must point down, got {right_twist:?}"
7868 );
7869 assert!(
7870 left_twist.y > 0.0,
7871 "left-hand Magnus must point up, got {left_twist:?}"
7872 );
7873 assert!((right_twist.y + left_twist.y).abs() < 1e-12);
7874 assert!(right_twist.x.abs() < 1e-12 && right_twist.z.abs() < 1e-12);
7875 assert!(left_twist.x.abs() < 1e-12 && left_twist.z.abs() < 1e-12);
7876 }
7877
7878 #[test]
7879 fn magnus_uses_velocity_corrected_muzzle_stability_gate() {
7880 let muzzle_velocity = 1_400.0 / 3.28084;
7881 let inputs = BallisticInputs {
7882 muzzle_velocity,
7883 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
7884 bullet_diameter: 0.308 * 0.0254,
7885 bullet_length: 1.215 * 0.0254,
7886 twist_rate: 15.0,
7887 enable_magnus: true,
7888 ..BallisticInputs::default()
7889 };
7890 let solver = TrajectorySolver::new(
7891 inputs.clone(),
7892 WindConditions::default(),
7893 AtmosphericConditions::default(),
7894 );
7895
7896 let bare_sg = crate::spin_drift::miller_stability(0.308, 168.0, 15.0, 1.215);
7897 let canonical_sg = solver.effective_spin_drift_sg();
7898 assert!(bare_sg > 1.0, "test requires bare Sg above the Magnus gate");
7899 assert!(
7900 canonical_sg < 1.0,
7901 "velocity-corrected Sg must be below the gate, got {canonical_sg}"
7902 );
7903
7904 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7905 let acceleration = solver.calculate_acceleration(
7906 &Vector3::zeros(),
7907 &Vector3::new(muzzle_velocity, 0.0, 0.0),
7908 &Vector3::zeros(),
7909 (temp_c, pressure_hpa, density / 1.225),
7910 );
7911 let mut baseline_inputs = inputs;
7912 baseline_inputs.enable_magnus = false;
7913 let baseline_solver = TrajectorySolver::new(
7914 baseline_inputs,
7915 WindConditions::default(),
7916 AtmosphericConditions::default(),
7917 );
7918 let baseline = baseline_solver.calculate_acceleration(
7919 &Vector3::zeros(),
7920 &Vector3::new(muzzle_velocity, 0.0, 0.0),
7921 &Vector3::zeros(),
7922 (temp_c, pressure_hpa, density / 1.225),
7923 );
7924
7925 assert_eq!(
7926 acceleration, baseline,
7927 "canonical Sg below 1 must suppress every Magnus acceleration component"
7928 );
7929 }
7930
7931 #[test]
7932 fn magnus_force_grows_as_fixed_spin_projectile_slows() {
7933 let inputs = BallisticInputs {
7934 muzzle_velocity: 800.0,
7935 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
7936 bullet_diameter: 0.308 * 0.0254,
7937 bullet_length: 1.215 * 0.0254,
7938 twist_rate: 12.0,
7939 enable_magnus: true,
7940 ..BallisticInputs::default()
7941 };
7942
7943 let magnus_acceleration = |speed_mps| {
7944 let evaluate = |enable_magnus| {
7945 let mut run_inputs = inputs.clone();
7946 run_inputs.enable_magnus = enable_magnus;
7947 let solver = TrajectorySolver::new(
7948 run_inputs,
7949 WindConditions::default(),
7950 AtmosphericConditions::default(),
7951 );
7952 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7953 solver
7954 .calculate_acceleration(
7955 &Vector3::zeros(),
7956 &Vector3::new(speed_mps, 0.0, 0.0),
7957 &Vector3::zeros(),
7958 (temp_c, pressure_hpa, density / 1.225),
7959 )
7960 .y
7961 };
7962 (evaluate(true) - evaluate(false)).abs()
7963 };
7964
7965 let fast = magnus_acceleration(200.0);
7966 let slow = magnus_acceleration(100.0);
7967 let ratio = slow / fast;
7968 let expected_ratio = 2.0_f64.powf(5.0 / 3.0);
7969
7970 assert!(fast > 0.0 && slow > 0.0, "fast={fast}, slow={slow}");
7971 assert!(
7972 (ratio - expected_ratio).abs() < 1e-3,
7973 "fixed-spin Magnus acceleration must grow downrange; slow/fast={ratio}, \
7974 expected={expected_ratio}"
7975 );
7976 }
7977}
7978
7979#[cfg(test)]
7980mod coriolis_direction_tests {
7981 use super::*;
7982 use std::f64::consts::FRAC_PI_2;
7983
7984 #[test]
7985 fn supersonic_crossing_flags_a_positive_range_sample() {
7986 use crate::trajectory_sampling::TrajectoryFlag;
7990
7991 for (solver_name, use_rk4, use_adaptive_rk45) in [
7992 ("Euler", false, false),
7993 ("RK4", true, false),
7994 ("RK45", true, true),
7995 ] {
7996 let inputs = BallisticInputs {
7997 muzzle_velocity: 850.0,
7998 bc_value: 0.2,
7999 bc_type: DragModel::G7,
8000 muzzle_angle: 0.03,
8001 enable_trajectory_sampling: true,
8002 sample_interval: 50.0,
8003 use_rk4,
8004 use_adaptive_rk45,
8005 ..BallisticInputs::default()
8006 };
8007 let mut solver = TrajectorySolver::new(
8008 inputs,
8009 WindConditions::default(),
8010 AtmosphericConditions::default(),
8011 );
8012 solver.set_max_range(2000.0);
8013 let samples = solver
8014 .solve()
8015 .expect("supersonic solve should succeed")
8016 .sampled_points
8017 .expect("sampling was enabled");
8018 let flagged_distances: Vec<_> = samples
8019 .iter()
8020 .filter(|sample| sample.flags.contains(&TrajectoryFlag::MachTransition))
8021 .map(|sample| sample.distance_m)
8022 .collect();
8023
8024 assert!(
8025 !flagged_distances.is_empty()
8026 && flagged_distances.iter().all(|distance| *distance > 0.0),
8027 "{solver_name} must flag genuine crossings only at positive range: {flagged_distances:?}"
8028 );
8029 }
8030 }
8031
8032 #[test]
8033 fn subsonic_launch_does_not_flag_a_muzzle_transition() {
8034 use crate::trajectory_sampling::TrajectoryFlag;
8035
8036 for (solver_name, use_rk4, use_adaptive_rk45) in [
8037 ("Euler", false, false),
8038 ("RK4", true, false),
8039 ("RK45", true, true),
8040 ] {
8041 let inputs = BallisticInputs {
8042 muzzle_velocity: 250.0,
8043 muzzle_angle: 0.02,
8044 enable_trajectory_sampling: true,
8045 sample_interval: 25.0,
8046 use_rk4,
8047 use_adaptive_rk45,
8048 ..BallisticInputs::default()
8049 };
8050 let mut solver = TrajectorySolver::new(
8051 inputs,
8052 WindConditions::default(),
8053 AtmosphericConditions::default(),
8054 );
8055 solver.set_max_range(300.0);
8056 let samples = solver
8057 .solve()
8058 .expect("subsonic solve should succeed")
8059 .sampled_points
8060 .expect("sampling was enabled");
8061
8062 assert!(
8063 samples
8064 .iter()
8065 .all(|sample| !sample.flags.contains(&TrajectoryFlag::MachTransition)),
8066 "{solver_name} marked a Mach transition for a launch already below Mach 1"
8067 );
8068 }
8069 }
8070
8071 #[test]
8072 fn mach_transition_tracker_requires_a_downward_crossing() {
8073 fn record(mach_values: &[f64]) -> Vec<f64> {
8074 let mut tracker = MachTransitionTracker::default();
8075 let mut distances = Vec::new();
8076 for (index, mach) in mach_values.iter().copied().enumerate() {
8077 tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
8078 }
8079 distances
8080 }
8081
8082 assert!(record(&[0.9, 0.8, 0.7]).is_empty());
8083 assert_eq!(record(&[1.1, 1.05, 0.99]), vec![20.0]);
8084 assert_eq!(record(&[1.2, 1.19, 1.0, 0.99]), vec![10.0, 30.0]);
8085 assert_eq!(record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]), vec![20.0, 30.0]);
8086 assert!(record(&[1.3, f64::NAN, 1.1]).is_empty());
8087 }
8088
8089 #[test]
8090 fn mach_transition_tracker_labels_0_9_without_touching_the_flat_vec() {
8091 fn record(mach_values: &[f64]) -> (Vec<f64>, MachTransitionTracker) {
8097 let mut tracker = MachTransitionTracker::default();
8098 let mut distances = Vec::new();
8099 for (index, mach) in mach_values.iter().copied().enumerate() {
8100 tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
8101 }
8102 (distances, tracker)
8103 }
8104
8105 let (distances, tracker) = record(&[0.9, 0.8, 0.7]);
8108 assert!(distances.is_empty()); assert_eq!(tracker.mach_1_2_distance_m, None);
8110 assert_eq!(tracker.mach_1_0_distance_m, None);
8111 assert_eq!(tracker.mach_0_9_distance_m, Some(10.0));
8112
8113 let (distances, tracker) = record(&[1.1, 1.05, 0.99]);
8115 assert_eq!(distances, vec![20.0]);
8116 assert_eq!(tracker.mach_1_2_distance_m, None);
8117 assert_eq!(tracker.mach_1_0_distance_m, Some(20.0));
8118 assert_eq!(tracker.mach_0_9_distance_m, None);
8119
8120 let (distances, tracker) = record(&[1.2, 1.19, 1.0, 0.99]);
8122 assert_eq!(distances, vec![10.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(10.0));
8124 assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
8125 assert_eq!(tracker.mach_0_9_distance_m, None);
8126
8127 let (distances, tracker) = record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]);
8130 assert_eq!(distances, vec![20.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(20.0));
8132 assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
8133 assert_eq!(tracker.mach_0_9_distance_m, Some(50.0));
8134 assert!(
8135 tracker.mach_1_2_distance_m < tracker.mach_1_0_distance_m
8136 && tracker.mach_1_0_distance_m < tracker.mach_0_9_distance_m,
8137 "labeled crossings must be strictly increasing downrange"
8138 );
8139
8140 let (distances, tracker) = record(&[1.3, f64::NAN, 1.1]);
8142 assert!(distances.is_empty());
8143 assert_eq!(tracker.mach_1_2_distance_m, None);
8144 assert_eq!(tracker.mach_1_0_distance_m, None);
8145 assert_eq!(tracker.mach_0_9_distance_m, None);
8146 }
8147
8148 #[test]
8149 fn humidity_percent_converts_and_clamps() {
8150 let mut i = BallisticInputs {
8152 humidity: 0.5,
8153 ..BallisticInputs::default()
8154 };
8155 assert!((i.humidity_percent() - 50.0).abs() < 1e-9, "0.5 -> 50%");
8156 i.humidity = 0.0;
8157 assert_eq!(i.humidity_percent(), 0.0);
8158 i.humidity = 1.0;
8159 assert_eq!(i.humidity_percent(), 100.0);
8160 i.humidity = 1.5; assert_eq!(i.humidity_percent(), 100.0);
8162 }
8163
8164 fn vertical_at(shot_azimuth: f64, range_m: f64) -> f64 {
8167 let inputs = BallisticInputs {
8168 muzzle_velocity: 800.0,
8169 bc_value: 0.5,
8170 bc_type: DragModel::G7,
8171 muzzle_angle: 0.02, enable_coriolis: true,
8173 latitude: Some(45.0),
8174 shot_azimuth,
8175 ground_threshold: f64::NEG_INFINITY, ..BallisticInputs::default()
8177 };
8178 let mut solver = TrajectorySolver::new(
8179 inputs,
8180 WindConditions::default(),
8181 AtmosphericConditions::default(),
8182 );
8183 solver.set_max_range(range_m + 50.0);
8184 let r = solver.solve().expect("solve");
8185 let pts = &r.points;
8186 for i in 1..pts.len() {
8187 if pts[i].position.x >= range_m {
8188 let p1 = &pts[i - 1];
8189 let p2 = &pts[i];
8190 let t = (range_m - p1.position.x) / (p2.position.x - p1.position.x);
8191 return p1.position.y + t * (p2.position.y - p1.position.y);
8192 }
8193 }
8194 panic!("range {range_m} not reached");
8195 }
8196
8197 #[test]
8202 fn eotvos_east_higher_than_west() {
8203 let range = 600.0;
8204 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!(
8208 east > west,
8209 "east ({east:.5}) must be higher than west ({west:.5}) at {range} m (Eötvös)"
8210 );
8211 assert!(
8212 east > north && north > west,
8213 "north ({north:.5}) must lie between east ({east:.5}) and west ({west:.5})"
8214 );
8215 assert!(
8216 (east - west) > 1e-3,
8217 "E-W vertical separation ({:.6} m) should be physically meaningful, not FP noise",
8218 east - west
8219 );
8220 }
8221
8222 #[test]
8230 fn labeled_mach_crossings_match_pinned_pre_change_flat_vec_across_solvers() {
8231 let cases = [
8233 ("Euler", false, false, 670.9878683238721_f64, 805.5274119916264_f64),
8234 ("RK4", true, false, 671.7257336844475_f64, 805.933409072171_f64),
8235 ("RK45", true, true, 672.4905711917901_f64, 806.5709746782849_f64),
8236 ];
8237
8238 for (solver_name, use_rk4, use_adaptive_rk45, expected_1_2, expected_1_0) in cases {
8239 let inputs = BallisticInputs {
8240 muzzle_velocity: 850.0,
8241 bc_value: 0.2,
8242 bc_type: DragModel::G7,
8243 muzzle_angle: 0.03,
8244 use_rk4,
8245 use_adaptive_rk45,
8246 ..BallisticInputs::default()
8247 };
8248 let mut solver = TrajectorySolver::new(
8249 inputs,
8250 WindConditions::default(),
8251 AtmosphericConditions::default(),
8252 );
8253 solver.set_max_range(2000.0);
8254 let result = solver.solve().expect("solve should succeed");
8255
8256 assert_eq!(
8257 result.mach_1_2_distance_m,
8258 Some(expected_1_2),
8259 "{solver_name}: mach_1_2_distance_m must match the pinned pre-change flat-Vec value"
8260 );
8261 assert_eq!(
8262 result.mach_1_0_distance_m,
8263 Some(expected_1_0),
8264 "{solver_name}: mach_1_0_distance_m must match the pinned pre-change flat-Vec value"
8265 );
8266
8267 let mach_1_2 = result.mach_1_2_distance_m.expect("crosses 1.2");
8268 let mach_1_0 = result.mach_1_0_distance_m.expect("crosses 1.0");
8269 let mach_0_9 = result
8270 .mach_0_9_distance_m
8271 .expect("this trajectory also goes past 0.9 within 2000 m");
8272 assert!(
8273 mach_1_2 < mach_1_0 && mach_1_0 < mach_0_9,
8274 "{solver_name}: labeled crossings must be strictly increasing downrange \
8275 (1.2={mach_1_2}, 1.0={mach_1_0}, 0.9={mach_0_9})"
8276 );
8277 }
8278 }
8279
8280 #[test]
8283 fn labeled_mach_crossings_are_none_for_a_fully_supersonic_trajectory() {
8284 for (solver_name, use_rk4, use_adaptive_rk45) in [
8285 ("Euler", false, false),
8286 ("RK4", true, false),
8287 ("RK45", true, true),
8288 ] {
8289 let inputs = BallisticInputs {
8290 muzzle_velocity: 850.0,
8291 bc_value: 0.2,
8292 bc_type: DragModel::G7,
8293 muzzle_angle: 0.03,
8294 use_rk4,
8295 use_adaptive_rk45,
8296 ..BallisticInputs::default()
8297 };
8298 let mut solver = TrajectorySolver::new(
8299 inputs,
8300 WindConditions::default(),
8301 AtmosphericConditions::default(),
8302 );
8303 solver.set_max_range(200.0);
8305 let result = solver.solve().expect("solve should succeed");
8306
8307 assert_eq!(
8308 result.mach_1_2_distance_m, None,
8309 "{solver_name}: must not report a 1.2 crossing that never happens"
8310 );
8311 assert_eq!(
8312 result.mach_1_0_distance_m, None,
8313 "{solver_name}: must not report a 1.0 crossing that never happens"
8314 );
8315 assert_eq!(
8316 result.mach_0_9_distance_m, None,
8317 "{solver_name}: must not report a 0.9 crossing that never happens"
8318 );
8319 }
8320 }
8321}
8322
8323#[cfg(test)]
8324mod cant_tests {
8325 use super::*;
8326
8327 fn base_inputs() -> BallisticInputs {
8328 BallisticInputs {
8329 muzzle_velocity: 800.0,
8330 bc_value: 0.5,
8331 bc_type: DragModel::G7,
8332 bullet_mass: 0.0109,
8333 bullet_diameter: 0.00782,
8334 bullet_length: 0.0309,
8335 sight_height: 0.05,
8336 twist_rate: 10.0,
8337 use_rk4: true,
8338 ..BallisticInputs::default()
8339 }
8340 }
8341
8342 fn solve_with(inputs: BallisticInputs, max_range: f64) -> TrajectoryResult {
8343 let mut s = TrajectorySolver::new(
8344 inputs,
8345 WindConditions::default(),
8346 AtmosphericConditions::default(),
8347 );
8348 s.set_max_range(max_range);
8349 s.solve().expect("solve")
8350 }
8351
8352 fn yz_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
8354 let pts = &result.points;
8355 for i in 1..pts.len() {
8356 if pts[i].position.x >= x {
8357 let (p1, p2) = (&pts[i - 1], &pts[i]);
8358 let dx = p2.position.x - p1.position.x;
8359 let t = if dx.abs() < 1e-12 { 0.0 } else { (x - p1.position.x) / dx };
8360 return (
8361 p1.position.y + t * (p2.position.y - p1.position.y),
8362 p1.position.z + t * (p2.position.z - p1.position.z),
8363 );
8364 }
8365 }
8366 panic!("trajectory never reached {x} m");
8367 }
8368
8369 #[test]
8370 fn cant_sign_clockwise_up_offset_goes_right_and_low() {
8371 let mut level = base_inputs();
8373 level.muzzle_angle = 0.003; let mut canted = level.clone();
8375 canted.cant_angle = 10f64.to_radians();
8376
8377 let (y0, z0) = yz_at(&solve_with(level, 400.0), 300.0);
8378 let (y1, z1) = yz_at(&solve_with(canted, 400.0), 300.0);
8379 assert!(z1 > z0 + 0.01, "clockwise cant must move POI right: z0={z0} z1={z1}");
8380 assert!(y1 < y0 - 0.001, "clockwise cant must move POI low: y0={y0} y1={y1}");
8381 }
8382
8383 #[test]
8384 fn pure_cant_shows_bore_offset_near_range() {
8385 let mut i = base_inputs();
8388 i.muzzle_angle = 0.0;
8389 i.cant_angle = 10f64.to_radians();
8390 let sh = i.sight_height;
8391 let r = solve_with(i, 60.0);
8392 let first = &r.points[1]; let expected = -sh * 10f64.to_radians().sin();
8394 assert!(
8395 (first.position.z - expected).abs() < 0.005,
8396 "near-muzzle lateral {} should be ~bore offset {expected}",
8397 first.position.z
8398 );
8399 }
8400
8401 #[test]
8402 fn zero_angle_is_independent_of_cant() {
8403 let a = base_inputs();
8404 let mut b = base_inputs();
8405 b.cant_angle = 15f64.to_radians();
8406 let za = calculate_zero_angle(a.clone(), 100.0, 0.0).expect("zero a");
8407 let zb = calculate_zero_angle(b.clone(), 100.0, 0.0).expect("zero b");
8408 assert_eq!(za.to_bits(), zb.to_bits(), "zeroing must ignore cant: {za} vs {zb}");
8409 let _ = (a.cant_angle, b.cant_angle);
8411 }
8412
8413 #[test]
8414 fn nonfinite_cant_is_rejected() {
8415 let mut i = base_inputs();
8416 i.cant_angle = f64::NAN;
8417 let s = TrajectorySolver::new(i, WindConditions::default(), AtmosphericConditions::default());
8418 assert!(s.solve().is_err());
8419 }
8420
8421 #[test]
8422 fn incline_and_cant_compose_without_breaking() {
8423 let mut flat = base_inputs();
8425 flat.muzzle_angle = 0.003;
8426 flat.shooting_angle = 15f64.to_radians();
8427 let mut canted = flat.clone();
8428 canted.cant_angle = 10f64.to_radians();
8429 let (_, z_flat) = yz_at(&solve_with(flat, 400.0), 300.0);
8430 let (_, z_cant) = yz_at(&solve_with(canted, 400.0), 300.0);
8431 assert!(z_cant > z_flat, "cant must still deflect right on an incline");
8432 }
8433}
8434
8435#[cfg(test)]
8436mod vertical_wind_tests {
8437 use super::*;
8438
8439 fn base_inputs() -> BallisticInputs {
8440 BallisticInputs {
8441 muzzle_velocity: 800.0,
8442 bc_value: 0.5,
8443 bc_type: DragModel::G7,
8444 bullet_mass: 0.0109,
8445 bullet_diameter: 0.00782,
8446 bullet_length: 0.0309,
8447 sight_height: 0.05,
8448 twist_rate: 10.0,
8449 use_rk4: true,
8450 ..BallisticInputs::default()
8451 }
8452 }
8453
8454 fn y_at(result: &TrajectoryResult, x: f64) -> f64 {
8456 let pts = &result.points;
8457 for i in 1..pts.len() {
8458 if pts[i].position.x >= x {
8459 let (p1, p2) = (&pts[i - 1], &pts[i]);
8460 let dx = p2.position.x - p1.position.x;
8461 let t = if dx.abs() < 1e-12 { 0.0 } else { (x - p1.position.x) / dx };
8462 return p1.position.y + t * (p2.position.y - p1.position.y);
8463 }
8464 }
8465 panic!("trajectory never reached {x} m");
8466 }
8467
8468 fn solve_with(inputs: BallisticInputs, wind: WindConditions, max_range: f64) -> TrajectoryResult {
8469 let mut s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
8470 s.set_max_range(max_range);
8471 s.solve().expect("solve")
8472 }
8473
8474 #[test]
8475 fn updraft_raises_poi_downrange() {
8476 let calm_inputs = base_inputs();
8479 let calm_wind = WindConditions::default();
8480 let updraft = WindConditions {
8481 vertical_speed: 5.0,
8482 ..Default::default()
8483 };
8484
8485 let calm = solve_with(calm_inputs.clone(), calm_wind, 500.0);
8486 let updraft_result = solve_with(calm_inputs, updraft, 500.0);
8487
8488 let y_calm = y_at(&calm, 400.0);
8489 let y_updraft = y_at(&updraft_result, 400.0);
8490 assert!(
8491 y_updraft > y_calm,
8492 "5 m/s updraft must raise POI at 400m: calm={y_calm}, updraft={y_updraft}"
8493 );
8494 }
8495
8496 #[test]
8497 fn zero_vertical_is_default_and_finite_required() {
8498 assert_eq!(WindConditions::default().vertical_speed, 0.0);
8499
8500 let inputs = base_inputs();
8501 let wind = WindConditions {
8502 vertical_speed: f64::NAN,
8503 ..Default::default()
8504 };
8505 let s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
8506 assert!(
8507 s.solve().is_err(),
8508 "NaN wind.vertical_speed must be rejected by validate_for_solve"
8509 );
8510 }
8511}
8512
8513#[cfg(test)]
8515mod bc_reference_standard_tests {
8516 use super::*;
8517
8518 fn base_inputs() -> BallisticInputs {
8519 BallisticInputs {
8520 muzzle_velocity: 800.0,
8521 bc_value: 0.5,
8522 bc_type: DragModel::G7,
8523 bullet_mass: 0.0109,
8524 bullet_diameter: 0.00782,
8525 bullet_length: 0.0309,
8526 sight_height: 0.05,
8527 twist_rate: 10.0,
8528 use_rk4: true,
8529 ..BallisticInputs::default()
8530 }
8531 }
8532
8533 fn y_and_speed_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
8535 let pts = &result.points;
8536 for i in 1..pts.len() {
8537 if pts[i].position.x >= x {
8538 let (p1, p2) = (&pts[i - 1], &pts[i]);
8539 let dx = p2.position.x - p1.position.x;
8540 let t = if dx.abs() < 1e-12 {
8541 0.0
8542 } else {
8543 (x - p1.position.x) / dx
8544 };
8545 return (
8546 p1.position.y + t * (p2.position.y - p1.position.y),
8547 p1.velocity_magnitude + t * (p2.velocity_magnitude - p1.velocity_magnitude),
8548 );
8549 }
8550 }
8551 panic!("trajectory never reached {x} m");
8552 }
8553
8554 #[test]
8557 fn asm_to_icao_ratio_matches_documented_value() {
8558 assert!(
8559 (crate::constants::ASM_TO_ICAO_BC - 0.98237).abs() < 1e-5,
8560 "ASM_TO_ICAO_BC = {} must equal 0.98237 to 5 decimal places",
8561 crate::constants::ASM_TO_ICAO_BC
8562 );
8563 assert_eq!(
8568 crate::constants::ASM_TO_ICAO_BC,
8569 crate::constants::ASM_DENSITY_LB_FT3 / crate::constants::ICAO_DENSITY_LB_FT3
8570 );
8571 }
8572
8573 #[test]
8576 fn default_bc_reference_standard_is_icao() {
8577 assert_eq!(
8578 BallisticInputs::default().bc_reference_standard,
8579 BcReferenceStandard::Icao
8580 );
8581 }
8582
8583 #[test]
8588 fn icao_reference_leaves_bc_value_bit_identical() {
8589 let raw_bc: f64 = 0.4372911; let inputs = BallisticInputs {
8591 bc_value: raw_bc,
8592 bc_reference_standard: BcReferenceStandard::Icao,
8593 ..base_inputs()
8594 };
8595 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8596 assert_eq!(solver.inputs.bc_value.to_bits(), raw_bc.to_bits());
8597 }
8598
8599 #[test]
8600 fn default_inputs_solve_is_unaffected_by_the_new_field_existing() {
8601 let a = TrajectorySolver::new(base_inputs(), WindConditions::default(), AtmosphericConditions::default())
8605 .solve()
8606 .expect("solve a");
8607 let b = TrajectorySolver::new(
8608 BallisticInputs { ..base_inputs() },
8609 WindConditions::default(),
8610 AtmosphericConditions::default(),
8611 )
8612 .solve()
8613 .expect("solve b");
8614 assert_eq!(a.impact_velocity.to_bits(), b.impact_velocity.to_bits());
8615 assert_eq!(a.max_range.to_bits(), b.max_range.to_bits());
8616 }
8617
8618 #[test]
8621 fn army_standard_metro_scales_bc_value_by_exactly_the_derived_ratio() {
8622 let raw_bc = 0.5;
8623 let inputs = BallisticInputs {
8624 bc_value: raw_bc,
8625 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8626 ..base_inputs()
8627 };
8628 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8629 assert_eq!(
8630 solver.inputs.bc_value,
8631 raw_bc * crate::constants::ASM_TO_ICAO_BC
8632 );
8633 }
8634
8635 #[test]
8636 fn army_standard_metro_scales_mach_keyed_bc_segments() {
8637 let inputs = BallisticInputs {
8638 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8639 bc_segments: Some(vec![(0.5, 0.40), (1.5, 0.30)]),
8640 ..base_inputs()
8641 };
8642 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8643 let segments = solver.inputs.bc_segments.as_ref().expect("segments");
8644 assert_eq!(segments[0], (0.5, 0.40 * crate::constants::ASM_TO_ICAO_BC));
8645 assert_eq!(segments[1], (1.5, 0.30 * crate::constants::ASM_TO_ICAO_BC));
8646 }
8647
8648 #[test]
8649 fn army_standard_metro_scales_velocity_keyed_bc_segments_data() {
8650 let inputs = BallisticInputs {
8651 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8652 bc_segments_data: Some(vec![
8653 crate::BCSegmentData {
8654 velocity_min: 0.0,
8655 velocity_max: 500.0,
8656 bc_value: 0.40,
8657 },
8658 crate::BCSegmentData {
8659 velocity_min: 500.0,
8660 velocity_max: 900.0,
8661 bc_value: 0.45,
8662 },
8663 ]),
8664 ..base_inputs()
8665 };
8666 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8667 let segments = solver.inputs.bc_segments_data.as_ref().expect("segments");
8668 assert_eq!(segments[0].bc_value, 0.40 * crate::constants::ASM_TO_ICAO_BC);
8669 assert_eq!(segments[1].bc_value, 0.45 * crate::constants::ASM_TO_ICAO_BC);
8670 assert_eq!(segments[0].velocity_min, 0.0);
8672 assert_eq!(segments[1].velocity_max, 900.0);
8673 }
8674
8675 #[test]
8680 fn army_standard_metro_moves_impact_in_the_more_drag_direction() {
8681 let solve_at = |standard: BcReferenceStandard| {
8682 let inputs = BallisticInputs {
8683 bc_value: 0.475,
8684 bc_reference_standard: standard,
8685 ..base_inputs()
8686 };
8687 let mut solver =
8688 TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8689 solver.set_max_range(500.0);
8690 solver.solve().expect("solve")
8691 };
8692
8693 let icao = solve_at(BcReferenceStandard::Icao);
8694 let asm = solve_at(BcReferenceStandard::ArmyStandardMetro);
8695
8696 let (y_icao, v_icao) = y_and_speed_at(&icao, 400.0);
8697 let (y_asm, v_asm) = y_and_speed_at(&asm, 400.0);
8698
8699 assert!(
8700 y_asm < y_icao,
8701 "ArmyStandardMetro must drop MORE (lower y) at 400m than Icao for the same raw \
8702 bc_value: icao_y={y_icao}, asm_y={y_asm}"
8703 );
8704 assert!(
8705 v_asm < v_icao,
8706 "ArmyStandardMetro must retain LESS velocity at 400m than Icao for the same raw \
8707 bc_value: icao_v={v_icao}, asm_v={v_asm}"
8708 );
8709 }
8710
8711 #[test]
8718 fn monte_carlo_inherits_the_normalized_bc_reference() {
8719 let base_inputs_asm = BallisticInputs {
8720 bc_value: 0.475,
8721 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8722 ..base_inputs()
8723 };
8724 let wind = WindConditions::default();
8725
8726 let mut direct_solver =
8732 TrajectorySolver::new(base_inputs_asm.clone(), wind.clone(), AtmosphericConditions::default());
8733 direct_solver.set_max_range(base_inputs_asm.target_distance.max(1000.0) * 2.0);
8734 let direct = direct_solver.solve().expect("direct solve");
8735
8736 let mc_params = MonteCarloParams {
8737 num_simulations: 1,
8738 velocity_std_dev: 0.0,
8739 angle_std_dev: 0.0,
8740 bc_std_dev: 0.0,
8741 wind_speed_std_dev: 0.0,
8742 target_distance: None,
8743 base_wind_speed: 0.0,
8744 base_wind_direction: 0.0,
8745 azimuth_std_dev: 0.0,
8746 };
8747 let mc = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
8748 base_inputs_asm,
8749 wind,
8750 mc_params,
8751 0.0,
8752 42,
8753 )
8754 .expect("monte carlo");
8755
8756 assert_eq!(mc.ranges.len(), 1);
8757 assert_eq!(
8758 mc.ranges[0].to_bits(),
8759 direct.max_range.to_bits(),
8760 "a zero-dispersion single MC sample must match a plain solve of the same \
8761 ASM-referenced inputs bit-for-bit"
8762 );
8763 assert_eq!(
8764 mc.impact_velocities[0].to_bits(),
8765 direct.impact_velocity.to_bits()
8766 );
8767 }
8768
8769 #[test]
8777 fn estimate_bc_fit_recovers_an_icao_referenced_bc() {
8778 let known_bc = 0.475;
8779 let velocity = 800.0;
8780 let mass = 0.0109;
8781 let diameter = 0.00782;
8782 let atmosphere = AtmosphericConditions::default();
8783
8784 let synth_inputs = BallisticInputs {
8785 muzzle_velocity: velocity,
8786 bc_value: known_bc,
8787 bc_type: DragModel::G7,
8788 bullet_mass: mass,
8789 bullet_diameter: diameter,
8790 bullet_length: 0.0309,
8791 sight_height: 0.05,
8792 twist_rate: 10.0,
8793 use_rk4: true,
8794 bc_reference_standard: BcReferenceStandard::Icao,
8795 ..BallisticInputs::default()
8796 };
8797 let mut solver = TrajectorySolver::new(synth_inputs, WindConditions::default(), atmosphere.clone());
8798 solver.set_max_range(500.0);
8799 let trajectory = solver.solve().expect("synthetic solve");
8800
8801 let points: Vec<(f64, f64)> = [100.0, 200.0, 300.0, 400.0]
8802 .iter()
8803 .map(|&d| {
8804 let (y, _) = {
8805 let pts = &trajectory.points;
8806 let mut found = None;
8807 for i in 1..pts.len() {
8808 if pts[i].position.x >= d {
8809 let (p1, p2) = (&pts[i - 1], &pts[i]);
8810 let dx = p2.position.x - p1.position.x;
8811 let t = if dx.abs() < 1e-12 {
8812 0.0
8813 } else {
8814 (d - p1.position.x) / dx
8815 };
8816 found = Some((
8817 p1.position.y + t * (p2.position.y - p1.position.y),
8818 0.0,
8819 ));
8820 break;
8821 }
8822 }
8823 found.expect("trajectory reached observation distance")
8824 };
8825 (d, -y) })
8827 .collect();
8828
8829 let estimate = estimate_bc_fit(
8830 velocity,
8831 mass,
8832 diameter,
8833 &points,
8834 DragModel::G7,
8835 BcFitMode::Drop,
8836 atmosphere,
8837 None,
8838 0.05,
8839 )
8840 .expect("fit should converge");
8841
8842 assert!(
8843 (estimate.bc - known_bc).abs() < 0.02,
8844 "fit should recover the known ICAO-referenced bc={known_bc}, got {}",
8845 estimate.bc
8846 );
8847 }
8848
8849 #[test]
8852 fn custom_drag_table_makes_bc_reference_standard_numerically_inert() {
8853 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])
8854 .expect("valid table");
8855
8856 let solve_with = |standard: BcReferenceStandard| {
8857 let inputs = BallisticInputs {
8858 bc_value: 0.5, bc_reference_standard: standard,
8860 custom_drag_table: Some(table.clone()),
8861 ..base_inputs()
8862 };
8863 let mut solver =
8864 TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
8865 solver.set_max_range(500.0);
8866 solver.solve().expect("solve")
8867 };
8868
8869 let icao = solve_with(BcReferenceStandard::Icao);
8870 let asm = solve_with(BcReferenceStandard::ArmyStandardMetro);
8871
8872 assert_eq!(
8873 icao.impact_velocity.to_bits(),
8874 asm.impact_velocity.to_bits(),
8875 "a custom drag table must make bc_reference_standard fully inert"
8876 );
8877 assert_eq!(icao.max_range.to_bits(), asm.max_range.to_bits());
8878 }
8879
8880 #[test]
8881 fn custom_drag_table_inert_warning_fires_only_for_army_standard_metro_with_a_table() {
8882 let table = crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0], vec![0.3, 0.4, 0.3])
8883 .expect("valid table");
8884
8885 let no_table_icao = base_inputs();
8887 assert!(no_table_icao.bc_reference_standard_inert_warning().is_none());
8888 let no_table_asm = BallisticInputs {
8889 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8890 ..base_inputs()
8891 };
8892 assert!(no_table_asm.bc_reference_standard_inert_warning().is_none());
8893
8894 let table_icao = BallisticInputs {
8896 custom_drag_table: Some(table.clone()),
8897 ..base_inputs()
8898 };
8899 assert!(table_icao.bc_reference_standard_inert_warning().is_none());
8900
8901 let table_asm = BallisticInputs {
8903 custom_drag_table: Some(table),
8904 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8905 ..base_inputs()
8906 };
8907 let warning = table_asm
8908 .bc_reference_standard_inert_warning()
8909 .expect("must warn");
8910 assert!(warning.contains("--bc-reference"));
8911 assert!(warning.contains("--drag-table"));
8912 }
8913}
8914
8915#[cfg(test)]
8921mod effective_drag_coefficient_tests {
8922 use super::*;
8923
8924 fn inputs_175gr_g7() -> BallisticInputs {
8925 let mut inputs = BallisticInputs {
8926 bc_value: 0.243,
8927 bc_type: DragModel::G7,
8928 muzzle_velocity: 823.0,
8929 ..Default::default()
8930 };
8931 inputs.bullet_mass = 175.0 * crate::constants::GRAINS_TO_KG;
8935 inputs.bullet_diameter = 0.308 * 0.0254;
8936 inputs.weight_grains = 175.0;
8937 inputs.caliber_inches = 0.308;
8938 inputs
8939 }
8940
8941 fn solver(inputs: BallisticInputs) -> TrajectorySolver {
8942 TrajectorySolver::new(
8943 inputs,
8944 WindConditions::default(),
8945 AtmosphericConditions::default(),
8946 )
8947 }
8948
8949 #[test]
8953 fn reports_the_projectiles_own_cd_not_the_reference_tables() {
8954 let inputs = inputs_175gr_g7();
8955 let sd = inputs.sectional_density_lb_in2().expect("SD");
8956 let solver = solver(inputs);
8957
8958 let sos = 340.0;
8959 let velocity = 800.0;
8960 let mach = velocity / sos;
8961
8962 let reference = crate::drag::get_drag_coefficient(mach, &DragModel::G7);
8963 let reported = solver
8964 .effective_drag_coefficient(velocity, sos)
8965 .expect("mass and diameter are set");
8966
8967 let expected = reference * sd / 0.243;
8968 assert!(
8969 (reported - expected).abs() < 1e-12,
8970 "reported {reported} != Cd_ref * SD / BC {expected}"
8971 );
8972 assert!(
8975 (reported - reference).abs() > 1e-6,
8976 "form factor collapsed to 1; this fixture no longer distinguishes the two values"
8977 );
8978 }
8979
8980 #[test]
8983 fn a_custom_drag_table_passes_through_unscaled() {
8984 let mut inputs = inputs_175gr_g7();
8985 inputs.custom_drag_table = Some(crate::drag::DragTable::new(
8986 vec![0.5, 3.0],
8987 vec![0.15, 0.40],
8988 ));
8989 let solver = solver(inputs);
8990
8991 let sos = 340.0;
8992 let velocity = 0.9 * sos;
8993 let table_value = solver
8994 .inputs
8995 .custom_drag_table
8996 .as_ref()
8997 .expect("table")
8998 .interpolate(0.9);
8999
9000 let reported = solver
9001 .effective_drag_coefficient(velocity, sos)
9002 .expect("mass and diameter are set");
9003 assert!(
9004 (reported - table_value).abs() < 1e-12,
9005 "custom table Cd {table_value} was rescaled to {reported}"
9006 );
9007 }
9008
9009 #[test]
9012 fn a_velocity_segmented_bc_steps_the_reported_cd() {
9013 let mut inputs = inputs_175gr_g7();
9014 inputs.use_bc_segments = true;
9015 inputs.bc_segments_data = Some(vec![
9016 crate::BCSegmentData { velocity_min: 2400.0, velocity_max: 4000.0, bc_value: 0.243 },
9017 crate::BCSegmentData { velocity_min: 0.0, velocity_max: 2400.0, bc_value: 0.200 },
9018 ]);
9019 let solver = solver(inputs);
9020
9021 let sos = 340.0;
9022 let above = solver.effective_drag_coefficient(2500.0 / 3.28084, sos).expect("cd");
9024 let below = solver.effective_drag_coefficient(2300.0 / 3.28084, sos).expect("cd");
9025
9026 assert!(
9028 below > above,
9029 "expected the 0.200 band to report a higher Cd than the 0.243 band; got {below} vs {above}"
9030 );
9031 }
9032
9033 #[test]
9036 fn is_absent_when_sectional_density_is_unknown() {
9037 let mut inputs = inputs_175gr_g7();
9038 inputs.weight_grains = 0.0;
9039 inputs.bullet_mass = 0.0;
9040 let solver = solver(inputs);
9041 assert!(solver.effective_drag_coefficient(800.0, 340.0).is_none());
9042 }
9043
9044 #[test]
9049 fn the_json_emit_rule_is_flag_gated_and_absent_when_cd_is_unknown() {
9050 let mut point = TrajectoryPoint {
9051 time: 0.0,
9052 position: nalgebra::Vector3::new(0.0, 0.0, 0.0),
9053 velocity_magnitude: 800.0,
9054 kinetic_energy: 3000.0,
9055 drag_coefficient: Some(0.31),
9056 };
9057 assert_eq!(point.drag_coefficient_json_value(true), Some(0.31));
9058 assert_eq!(
9059 point.drag_coefficient_json_value(false),
9060 None,
9061 "without the flag the key must not exist, so default JSON stays byte-identical"
9062 );
9063 point.drag_coefficient = None;
9064 assert_eq!(
9065 point.drag_coefficient_json_value(true),
9066 None,
9067 "unknown sectional density must yield an ABSENT key, not null"
9068 );
9069 }
9070
9071 #[test]
9073 fn every_point_of_a_solved_trajectory_carries_the_value() {
9074 let mut solver = solver(inputs_175gr_g7());
9075 solver.set_max_range(300.0);
9076 let result = solver.solve().expect("solve");
9077 assert!(!result.points.is_empty());
9078 assert!(
9079 result.points.iter().all(|p| p.drag_coefficient.is_some()),
9080 "the post-integration pass missed at least one point"
9081 );
9082 }
9083}