1use crate::cluster_bc::ClusterBCDegradation;
3use crate::mc_stats::{wilson_interval, BernoulliConfidenceSequence, ConfidenceLevel, Welford};
4use crate::pitch_damping::{calculate_pitch_damping_coefficient, PitchDampingCoefficients};
5use crate::precession_nutation::{
6 calculate_combined_angular_motion, projectile_moments_of_inertia, AngularState,
7 PrecessionNutationParams,
8};
9use crate::trajectory_observation::{bracket_param, Bracket, TrajectoryTermination};
10use crate::trajectory_sampling::{
11 projected_sample_count, sample_trajectory, TrajectoryData, TrajectoryOutputs, TrajectorySample,
12};
13use crate::wind_shear::WindShearModel;
14use crate::DragModel;
15use nalgebra::{Vector3, Vector6};
16use std::error::Error;
17use std::fmt;
18
19#[derive(Debug, Clone, Copy, PartialEq)]
29#[cfg_attr(feature = "cli", derive(clap::ValueEnum))]
30pub enum UnitSystem {
31 Metric,
33 Imperial,
35}
36
37#[derive(Debug, Clone, Copy, PartialEq)]
39pub enum OutputFormat {
40 Table,
41 Json,
42 Csv,
43}
44
45#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
64#[cfg_attr(feature = "cli", derive(clap::ValueEnum))]
65pub enum BcReferenceStandard {
66 #[default]
69 Icao,
70 ArmyStandardMetro,
73}
74
75pub const BC_REFERENCE_STANDARD_INERT_WARNING: &str =
79 "warning: --bc-reference army-standard-metro has no effect together with a custom drag \
80 table (--drag-table): the deck's Cd is divided by sectional density, not a BC value, so \
81 no BC-reference conversion applies";
82
83#[derive(Debug)]
85pub struct BallisticsError {
86 message: String,
87}
88
89impl fmt::Display for BallisticsError {
90 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
91 write!(f, "{}", self.message)
92 }
93}
94
95impl Error for BallisticsError {}
96
97impl From<String> for BallisticsError {
98 fn from(msg: String) -> Self {
99 BallisticsError { message: msg }
100 }
101}
102
103impl From<&str> for BallisticsError {
104 fn from(msg: &str) -> Self {
105 BallisticsError {
106 message: msg.to_string(),
107 }
108 }
109}
110
111#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
119pub enum DropsReference {
120 #[default]
122 Los,
123 Target,
126}
127
128#[derive(Debug, Clone)]
132pub struct BallisticInputs {
133 pub bc_value: f64, pub bc_type: DragModel, pub bc_reference_standard: BcReferenceStandard,
141 pub bullet_mass: f64, pub muzzle_velocity: f64, pub bullet_diameter: f64, pub bullet_length: f64, pub muzzle_angle: f64, pub target_distance: f64, pub azimuth_angle: f64, pub shot_azimuth: f64,
155 pub shooting_angle: f64, pub cant_angle: f64,
165 pub sight_height: f64, pub sight_offset_lateral_m: f64,
177 pub muzzle_height: f64, pub target_height: f64, pub zero_poi_vertical_m: f64,
188 pub zero_poi_horizontal_m: f64,
194 pub ground_threshold: f64, pub altitude: f64, pub temperature: f64, pub pressure: f64, pub humidity: f64,
205 pub latitude: Option<f64>, pub wind_speed: f64, pub wind_angle: f64, pub twist_rate: f64, pub is_twist_right: bool, pub caliber_inches: f64, pub weight_grains: f64, pub manufacturer: Option<String>, pub bullet_model: Option<String>, pub bullet_id: Option<String>, pub bullet_cluster: Option<usize>, pub use_rk4: bool, pub use_adaptive_rk45: bool, pub enable_advanced_effects: bool,
227 pub enable_magnus: bool, pub enable_coriolis: bool, pub use_powder_sensitivity: bool,
230 pub powder_temp_sensitivity: f64, pub powder_temp: f64, pub powder_temp_curve: Option<Vec<(f64, f64)>>,
239 pub powder_curve_temp_c: Option<f64>,
243 pub tipoff_yaw: f64, pub tipoff_decay_distance: f64, pub cd_delta2: f64,
251 pub use_bc_segments: bool,
254 pub bc_segments: Option<Vec<(f64, f64)>>, pub bc_segments_data: Option<Vec<crate::BCSegmentData>>, pub use_enhanced_spin_drift: bool,
257 pub use_form_factor: bool,
260 pub enable_wind_shear: bool,
261 pub wind_shear_model: String,
262 pub enable_trajectory_sampling: bool,
263 pub sample_interval: f64, pub drops_reference: DropsReference,
273 pub enable_pitch_damping: bool,
274 pub enable_precession_nutation: bool,
275 pub enable_aerodynamic_jump: bool,
278 pub use_cluster_bc: bool, pub custom_drag_table: Option<crate::drag::DragTable>,
282 pub cd_scale: f64,
290
291 pub bc_type_str: Option<String>,
293}
294
295impl BallisticInputs {
296 pub fn humidity_percent(&self) -> f64 {
301 (self.humidity * 100.0).clamp(0.0, 100.0)
302 }
303
304 pub fn windage_zero_bias_rad(&self, zero_distance_m: f64) -> f64 {
319 if zero_distance_m > 0.0 {
320 (self.zero_poi_horizontal_m + self.sight_offset_lateral_m) / zero_distance_m
321 } else {
322 0.0
323 }
324 }
325
326 pub fn sectional_density_lb_in2(&self) -> Option<f64> {
332 let weight_gr = if self.weight_grains > 0.0 {
333 self.weight_grains
334 } else {
335 self.bullet_mass / crate::constants::GRAINS_TO_KG };
337 let diameter_in = if self.caliber_inches > 0.0 {
338 self.caliber_inches
339 } else {
340 self.bullet_diameter / 0.0254 };
342 if weight_gr > 0.0 && diameter_in > 0.0 {
343 Some(weight_gr / 7000.0 / (diameter_in * diameter_in))
344 } else {
345 None
346 }
347 }
348
349 pub fn custom_drag_denominator(&self, fallback_bc: f64) -> f64 {
361 match self.sectional_density_lb_in2() {
362 Some(sd) => sd,
363 None => {
364 static WARN_ONCE: std::sync::Once = std::sync::Once::new();
365 WARN_ONCE.call_once(|| {
366 eprintln!(
367 "Warning: custom drag table active but bullet mass/diameter are \
368 unavailable; falling back to bc_value for the retardation denominator"
369 );
370 });
371 fallback_bc
372 }
373 }
374 }
375
376 pub fn bc_reference_standard_inert_warning(&self) -> Option<&'static str> {
387 if self.custom_drag_table.is_some()
388 && matches!(
389 self.bc_reference_standard,
390 BcReferenceStandard::ArmyStandardMetro
391 )
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 + t * (points[hi].position.y - points[lo].position.y),
1493 )
1494 }
1495 }
1496 }
1497
1498 let mut inclined = self.clone();
1502 inclined.inputs.enable_trajectory_sampling = false;
1503 let inclined_result = inclined.solve().ok()?;
1504 let los_height = inclined_result.line_of_sight_height_m;
1505 let inclined_drop = los_height - path_y_at(&inclined_result.points, target_range_m)?;
1506 let correction = inclined_drop / target_range_m;
1507 if correction <= 0.0 {
1508 return None;
1509 }
1510
1511 let mut flat = self.clone();
1513 flat.inputs.enable_trajectory_sampling = false;
1514 flat.inputs.shooting_angle = 0.0;
1515 let flat_result = flat.solve().ok()?;
1516 let flat_correction_at = |range_m: f64| -> Option<f64> {
1517 Some((los_height - path_y_at(&flat_result.points, range_m)?) / range_m)
1518 };
1519
1520 let flat_terminal_m = flat_result.points.last().map(|p| p.position.x)?;
1528 let mut low = zero_distance_m.max(1.0);
1529 let mut high = target_range_m.min(flat_terminal_m);
1530 if high <= low {
1531 return None;
1532 }
1533 if flat_correction_at(low)? - correction > 0.0 {
1534 return None; }
1536 if flat_correction_at(high)? - correction < 0.0 {
1537 return None; }
1539 for _ in 0..60 {
1540 let mid = 0.5 * (low + high);
1541 let error = flat_correction_at(mid)? - correction;
1542 if error.abs() == 0.0 {
1543 return Some(mid);
1544 }
1545 if error < 0.0 {
1546 low = mid;
1547 } else {
1548 high = mid;
1549 }
1550 if high - low < 0.01 {
1551 break;
1552 }
1553 }
1554 Some(0.5 * (low + high))
1555 }
1556
1557 fn validate_for_solve(&self) -> Result<(), BallisticsError> {
1563 let require_finite = |name: &str, value: f64| {
1564 if value.is_finite() {
1565 Ok(())
1566 } else {
1567 Err(BallisticsError::from(format!("{name} must be finite")))
1568 }
1569 };
1570 let require_positive = |name: &str, value: f64| {
1571 if value.is_finite() && value > 0.0 {
1572 Ok(())
1573 } else {
1574 Err(BallisticsError::from(format!(
1575 "{name} must be finite and greater than zero"
1576 )))
1577 }
1578 };
1579
1580 if self.inputs.custom_drag_table.is_none() {
1586 require_positive("bc_value", self.inputs.bc_value)?;
1587 }
1588 require_positive("bullet_mass", self.inputs.bullet_mass)?;
1589 require_positive("bullet_diameter", self.inputs.bullet_diameter)?;
1590 require_positive("muzzle_velocity", self.inputs.muzzle_velocity)?;
1591 require_positive("cd_scale", self.inputs.cd_scale)?;
1597
1598 require_finite("muzzle_angle", self.inputs.muzzle_angle)?;
1599 require_finite("azimuth_angle", self.inputs.azimuth_angle)?;
1600 require_finite("shooting_angle", self.inputs.shooting_angle)?;
1601 require_finite("cant_angle", self.inputs.cant_angle)?;
1602 require_finite("muzzle_height", self.inputs.muzzle_height)?;
1603
1604 for (name, value) in [
1608 ("zero_poi_vertical_m", self.inputs.zero_poi_vertical_m),
1609 ("zero_poi_horizontal_m", self.inputs.zero_poi_horizontal_m),
1610 ] {
1611 require_finite(name, value)?;
1612 if value.abs() >= 1.0 {
1613 return Err(BallisticsError::from(format!(
1614 "{name} must be smaller than 1.0 m in magnitude (it is a linear POI \
1615 offset at the zero range, in meters)"
1616 )));
1617 }
1618 }
1619
1620 require_finite("sight_offset_lateral_m", self.inputs.sight_offset_lateral_m)?;
1624 if self.inputs.sight_offset_lateral_m.abs() >= 0.5 {
1625 return Err(BallisticsError::from(
1626 "sight_offset_lateral_m must be smaller than 0.5 m in magnitude (it is \
1627 the lateral sight-to-bore mount offset, in meters)",
1628 ));
1629 }
1630
1631 if !(self.inputs.ground_threshold.is_finite()
1634 || self.inputs.ground_threshold == f64::NEG_INFINITY)
1635 {
1636 return Err(BallisticsError::from(
1637 "ground_threshold must be finite or negative infinity",
1638 ));
1639 }
1640
1641 match &self.wind_sock {
1642 Some(wind_sock) => wind_sock
1643 .validate_segments()
1644 .map_err(BallisticsError::from)?,
1645 None => {
1646 require_finite("wind.speed", self.wind.speed)?;
1647 require_finite("wind.direction", self.wind.direction)?;
1648 require_finite("wind.vertical_speed", self.wind.vertical_speed)?;
1649 }
1650 }
1651
1652 require_finite("atmosphere.temperature", self.atmosphere.temperature)?;
1653 require_finite("atmosphere.pressure", self.atmosphere.pressure)?;
1654 require_finite("atmosphere.humidity", self.atmosphere.humidity)?;
1655 require_finite("atmosphere.altitude", self.atmosphere.altitude)?;
1656
1657 require_positive("max_range", self.max_range)?;
1658 if !self.inputs.use_rk4 || !self.inputs.use_adaptive_rk45 {
1661 require_positive("time_step", self.time_step)?;
1662 }
1663
1664 if self.inputs.enable_trajectory_sampling {
1665 require_finite("sight_height", self.inputs.sight_height)?;
1666 require_positive("sample_interval", self.inputs.sample_interval)?;
1667 projected_sample_count(self.max_range, self.inputs.sample_interval)?;
1668 }
1669
1670 if self.inputs.drops_reference == DropsReference::Target {
1674 require_finite("target_height", self.inputs.target_height)?;
1675 if self.inputs.shooting_angle.cos() <= 1e-9 {
1676 return Err(BallisticsError::from(
1677 "drops reference 'target' is undefined for shooting angles at or beyond 90 degrees",
1678 ));
1679 }
1680 }
1681
1682 if self.inputs.enable_coriolis {
1683 require_finite("shot_azimuth", self.inputs.shot_azimuth)?;
1684 if let Some(latitude) = self.inputs.latitude {
1685 require_finite("latitude", latitude)?;
1686 }
1687 }
1688
1689 Ok(())
1690 }
1691
1692 fn validate_result_sanity(&self, result: &TrajectoryResult) -> Result<(), BallisticsError> {
1699 let require_finite = |name: &str, value: f64| {
1700 if value.is_finite() {
1701 Ok(())
1702 } else {
1703 Err(BallisticsError::from(format!(
1704 "trajectory result contains non-finite {name}"
1705 )))
1706 }
1707 };
1708 let require_non_negative = |name: &str, value: f64| {
1709 if value >= 0.0 {
1710 Ok(())
1711 } else {
1712 Err(BallisticsError::from(format!(
1713 "trajectory result contains non-physical negative {name} ({value})"
1714 )))
1715 }
1716 };
1717 let require_indexed_finite = |collection: &str, index: usize, field: &str, value: f64| {
1718 if value.is_finite() {
1719 Ok(())
1720 } else {
1721 Err(BallisticsError::from(format!(
1722 "trajectory result contains non-finite {collection}[{index}].{field}"
1723 )))
1724 }
1725 };
1726 let require_indexed_non_negative = |collection: &str,
1727 index: usize,
1728 field: &str,
1729 value: f64| {
1730 if value >= 0.0 {
1731 Ok(())
1732 } else {
1733 Err(BallisticsError::from(format!(
1734 "trajectory result contains non-physical negative {collection}[{index}].{field} ({value})"
1735 )))
1736 }
1737 };
1738
1739 require_finite("max_range", result.max_range)?;
1740 require_finite("max_height", result.max_height)?;
1741 require_finite("time_of_flight", result.time_of_flight)?;
1742 require_finite("impact_velocity", result.impact_velocity)?;
1743 require_finite("impact_energy", result.impact_energy)?;
1744 require_finite("projectile_mass_kg", result.projectile_mass_kg)?;
1745 require_finite("line_of_sight_height_m", result.line_of_sight_height_m)?;
1746 require_finite(
1747 "station_speed_of_sound_mps",
1748 result.station_speed_of_sound_mps,
1749 )?;
1750
1751 require_non_negative("max_range", result.max_range)?;
1755 require_non_negative("time_of_flight", result.time_of_flight)?;
1756 require_non_negative("impact_velocity", result.impact_velocity)?;
1757 require_non_negative("impact_energy", result.impact_energy)?;
1758 require_non_negative("projectile_mass_kg", result.projectile_mass_kg)?;
1759 require_non_negative(
1760 "station_speed_of_sound_mps",
1761 result.station_speed_of_sound_mps,
1762 )?;
1763
1764 for (index, point) in result.points.iter().enumerate() {
1765 require_indexed_finite("points", index, "time", point.time)?;
1766 require_indexed_finite("points", index, "position.x", point.position.x)?;
1767 require_indexed_finite("points", index, "position.y", point.position.y)?;
1768 require_indexed_finite("points", index, "position.z", point.position.z)?;
1769 require_indexed_finite(
1770 "points",
1771 index,
1772 "velocity_magnitude",
1773 point.velocity_magnitude,
1774 )?;
1775 require_indexed_finite("points", index, "kinetic_energy", point.kinetic_energy)?;
1776 require_indexed_non_negative("points", index, "time", point.time)?;
1777 require_indexed_non_negative(
1778 "points",
1779 index,
1780 "velocity_magnitude",
1781 point.velocity_magnitude,
1782 )?;
1783 require_indexed_non_negative("points", index, "kinetic_energy", point.kinetic_energy)?;
1784 }
1785
1786 if let Some(samples) = &result.sampled_points {
1787 for (index, sample) in samples.iter().enumerate() {
1788 require_indexed_finite("sampled_points", index, "distance_m", sample.distance_m)?;
1789 require_indexed_finite("sampled_points", index, "drop_m", sample.drop_m)?;
1790 require_indexed_finite(
1791 "sampled_points",
1792 index,
1793 "wind_drift_m",
1794 sample.wind_drift_m,
1795 )?;
1796 require_indexed_finite(
1797 "sampled_points",
1798 index,
1799 "velocity_mps",
1800 sample.velocity_mps,
1801 )?;
1802 require_indexed_finite("sampled_points", index, "energy_j", sample.energy_j)?;
1803 require_indexed_finite("sampled_points", index, "time_s", sample.time_s)?;
1804 }
1805 }
1806
1807 for (name, value) in [
1808 ("min_pitch_damping", result.min_pitch_damping),
1809 ("transonic_mach", result.transonic_mach),
1810 ("max_yaw_angle", result.max_yaw_angle),
1811 ("max_precession_angle", result.max_precession_angle),
1812 ] {
1813 if let Some(value) = value {
1814 require_finite(name, value)?;
1815 }
1816 }
1817
1818 if let Some(state) = result.angular_state {
1819 for (name, value) in [
1820 ("angular_state.pitch_angle", state.pitch_angle),
1821 ("angular_state.yaw_angle", state.yaw_angle),
1822 ("angular_state.pitch_rate", state.pitch_rate),
1823 ("angular_state.yaw_rate", state.yaw_rate),
1824 ("angular_state.precession_angle", state.precession_angle),
1825 ("angular_state.nutation_phase", state.nutation_phase),
1826 ] {
1827 require_finite(name, value)?;
1828 }
1829 }
1830
1831 if let Some(jump) = result.aerodynamic_jump {
1832 for (name, value) in [
1833 ("aerodynamic_jump.vertical_jump_moa", jump.vertical_jump_moa),
1834 (
1835 "aerodynamic_jump.horizontal_jump_moa",
1836 jump.horizontal_jump_moa,
1837 ),
1838 ("aerodynamic_jump.jump_angle_rad", jump.jump_angle_rad),
1839 (
1840 "aerodynamic_jump.magnus_component_moa",
1841 jump.magnus_component_moa,
1842 ),
1843 ("aerodynamic_jump.yaw_component_moa", jump.yaw_component_moa),
1844 (
1845 "aerodynamic_jump.stabilization_factor",
1846 jump.stabilization_factor,
1847 ),
1848 ] {
1849 require_finite(name, value)?;
1850 }
1851 }
1852
1853 Ok(())
1854 }
1855
1856 fn validate_integration_state(
1869 &self,
1870 position: &Vector3<f64>,
1871 velocity: &Vector3<f64>,
1872 time: f64,
1873 ) -> Result<(), BallisticsError> {
1874 if !(position.iter().all(|value| value.is_finite())
1875 && velocity.iter().all(|value| value.is_finite())
1876 && time.is_finite())
1877 {
1878 return Err(BallisticsError::from(
1879 "trajectory integration produced a non-finite state (often from physically \
1880 extreme inputs — e.g. an absurd bore/muzzle height placing the launch far \
1881 from sea level, or a degenerate atmosphere; check those inputs, or set \
1882 --altitude explicitly)",
1883 ));
1884 }
1885
1886 let speed = velocity.magnitude();
1887 let budget = self.speed_budget(time);
1888 if speed > budget {
1889 return Err(BallisticsError::from(format!(
1890 "trajectory integration diverged: speed {speed:.3e} m/s at t={time:.6}s exceeds \
1891 the physical budget of {budget:.3e} m/s"
1892 )));
1893 }
1894 Ok(())
1895 }
1896
1897 fn speed_budget(&self, time: f64) -> f64 {
1902 let scalar_wind = self.wind.speed.abs() + self.wind.vertical_speed.abs();
1903 let wind_bound = match &self.wind_sock {
1904 Some(sock) => scalar_wind.max(sock.max_speed_mps()),
1905 None => scalar_wind,
1906 };
1907 2.0 * (self.inputs.muzzle_velocity + wind_bound + 10.0)
1908 + crate::constants::G_ACCEL_MPS2 * time
1909 }
1910
1911 fn push_trajectory_point(
1913 &self,
1914 points: &mut Vec<TrajectoryPoint>,
1915 point: TrajectoryPoint,
1916 ) -> Result<(), BallisticsError> {
1917 if points.len() >= self.max_trajectory_points {
1918 return Err(BallisticsError::from(format!(
1919 "trajectory point limit of {} exceeded",
1920 self.max_trajectory_points
1921 )));
1922 }
1923 points.push(point);
1924 Ok(())
1925 }
1926
1927 pub fn set_wind_segments(&mut self, segments: Vec<crate::wind::WindSegment>) {
1934 self.wind_sock = if segments.is_empty() {
1935 None
1936 } else {
1937 Some(crate::wind::WindSock::new(segments))
1938 };
1939 }
1940
1941 pub fn set_atmo_segments(&mut self, segments: Vec<crate::atmosphere::AtmoSegment>) {
1950 self.atmo_sock = if segments.is_empty() {
1951 None
1952 } else {
1953 Some(crate::atmosphere::AtmoSock::new(segments))
1954 };
1955 }
1956
1957 fn launch_angles_from(
1966 &self,
1967 aj: Option<&crate::aerodynamic_jump::AerodynamicJumpComponents>,
1968 ) -> (f64, f64) {
1969 let (mut elev, mut azim) = (self.inputs.muzzle_angle, self.inputs.azimuth_angle);
1970 if self.inputs.cant_angle != 0.0 {
1977 let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
1978 let (e0, a0) = (elev, azim);
1979 elev = e0 * cos_c - a0 * sin_c;
1980 azim = a0 * cos_c + e0 * sin_c;
1981 }
1982 match aj {
1983 Some(c) => {
1984 const MOA_PER_RAD: f64 = 3437.7467707849;
1986 (
1987 elev + c.vertical_jump_moa / MOA_PER_RAD,
1988 azim + c.horizontal_jump_moa / MOA_PER_RAD,
1989 )
1990 }
1991 None => (elev, azim),
1992 }
1993 }
1994
1995 fn aerodynamic_jump_components(
2003 &self,
2004 ) -> Option<crate::aerodynamic_jump::AerodynamicJumpComponents> {
2005 if !self.inputs.enable_aerodynamic_jump {
2006 return None;
2007 }
2008 let diameter_m = self.inputs.bullet_diameter;
2012 if !(self.inputs.twist_rate.is_finite()
2013 && self.inputs.twist_rate != 0.0
2014 && diameter_m.is_finite()
2015 && diameter_m > 0.0
2016 && self.inputs.bullet_length.is_finite()
2017 && self.inputs.bullet_length > 0.0
2018 && self.inputs.muzzle_velocity.is_finite())
2019 {
2020 return None;
2021 }
2022
2023 let (_, _, temp_c, pressure_hpa) = self.resolved_atmosphere();
2025 let sg = crate::stability::compute_stability_coefficient(
2026 &self.inputs,
2027 (self.atmosphere.altitude, temp_c, pressure_hpa, 0.0),
2028 );
2029 if !(sg.is_finite() && sg > 0.0) {
2030 return None;
2031 }
2032 let length_calibers = self.inputs.bullet_length / diameter_m;
2033
2034 const MS_TO_MPH: f64 = 2.236_936_292_054_4;
2040 let crosswind_from_right_mps = if let Some(sock) = &self.wind_sock {
2041 -sock.vector_for_range_stateless(0.0)[2]
2042 } else {
2043 self.wind.speed * self.wind.direction.sin()
2044 };
2045 let crosswind_from_right_mph = crosswind_from_right_mps * MS_TO_MPH;
2046
2047 let vertical_jump_moa = crate::aerodynamic_jump::litz_crosswind_jump_moa(
2048 sg,
2049 length_calibers,
2050 crosswind_from_right_mph,
2051 self.inputs.is_twist_right,
2052 );
2053 if !vertical_jump_moa.is_finite() {
2054 return None;
2055 }
2056
2057 const MOA_PER_RAD: f64 = 3437.7467707849;
2058 Some(crate::aerodynamic_jump::AerodynamicJumpComponents {
2059 vertical_jump_moa,
2060 horizontal_jump_moa: 0.0,
2062 jump_angle_rad: vertical_jump_moa.abs() / MOA_PER_RAD,
2063 magnus_component_moa: 0.0,
2064 yaw_component_moa: 0.0,
2065 stabilization_factor: (sg / 1.5).clamp(0.0, 1.0),
2066 })
2067 }
2068
2069 fn resolved_atmosphere(&self) -> (f64, f64, f64, f64) {
2070 let (temp_c, pressure_hpa) = match self.station_atmosphere_resolution {
2071 StationAtmosphereResolution::LegacyDefaultSentinels => {
2072 crate::atmosphere::resolve_station_conditions(
2073 self.atmosphere.temperature,
2074 self.atmosphere.pressure,
2075 self.atmosphere.altitude,
2076 )
2077 }
2078 StationAtmosphereResolution::Authoritative => {
2079 (self.atmosphere.temperature, self.atmosphere.pressure)
2080 }
2081 };
2082 let (density, speed_of_sound) = crate::atmosphere::calculate_atmosphere(
2083 self.atmosphere.altitude,
2084 Some(temp_c),
2085 Some(pressure_hpa),
2086 self.atmosphere.humidity,
2087 );
2088 (density, speed_of_sound, temp_c, pressure_hpa)
2089 }
2090
2091 fn precession_nutation_params(
2092 &self,
2093 velocity_mps: f64,
2094 air_density_kg_m3: f64,
2095 speed_of_sound_mps: f64,
2096 ) -> PrecessionNutationParams {
2097 let (spin_inertia, transverse_inertia) = self.precession_nutation_inertias;
2098 let spin_rate_rad_s = if self.inputs.twist_rate > 0.0 {
2099 let velocity_fps = velocity_mps * 3.28084;
2100 let twist_rate_ft = self.inputs.twist_rate / 12.0;
2101 (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI
2102 } else {
2103 0.0
2104 };
2105
2106 PrecessionNutationParams {
2107 mass_kg: self.inputs.bullet_mass,
2108 caliber_m: self.inputs.bullet_diameter,
2109 length_m: self.inputs.bullet_length,
2110 spin_rate_rad_s,
2111 spin_inertia,
2112 transverse_inertia,
2113 velocity_mps,
2114 air_density_kg_m3,
2115 mach: velocity_mps / speed_of_sound_mps,
2116 pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
2117 nutation_damping_factor: 0.05,
2118 }
2119 }
2120
2121 fn append_terminal_endpoint(
2128 &self,
2129 points: &mut Vec<TrajectoryPoint>,
2130 post_position: Vector3<f64>,
2131 post_velocity: Vector3<f64>,
2132 post_time: f64,
2133 max_height: &mut f64,
2134 ) -> Result<TrajectoryTermination, BallisticsError> {
2135 let previous = points
2136 .last()
2137 .cloned()
2138 .ok_or_else(|| BallisticsError::from("No trajectory points generated"))?;
2139
2140 let mut crossings = Vec::with_capacity(3);
2141 if previous.position.x < self.max_range && post_position.x >= self.max_range {
2142 let span = post_position.x - previous.position.x;
2143 if span.is_finite() && span > 0.0 {
2144 crossings.push((
2145 (self.max_range - previous.position.x) / span,
2146 TrajectoryTermination::MaxRange,
2147 ));
2148 }
2149 }
2150 if self.inputs.ground_threshold.is_finite()
2151 && previous.position.y > self.inputs.ground_threshold
2152 && post_position.y <= self.inputs.ground_threshold
2153 {
2154 let span = post_position.y - previous.position.y;
2155 if span.is_finite() && span < 0.0 {
2156 crossings.push((
2157 (self.inputs.ground_threshold - previous.position.y) / span,
2158 TrajectoryTermination::GroundThreshold,
2159 ));
2160 }
2161 }
2162 if previous.time < TRAJECTORY_TIME_LIMIT_S && post_time >= TRAJECTORY_TIME_LIMIT_S {
2163 let span = post_time - previous.time;
2164 if span.is_finite() && span > 0.0 {
2165 crossings.push((
2166 (TRAJECTORY_TIME_LIMIT_S - previous.time) / span,
2167 TrajectoryTermination::TimeLimit,
2168 ));
2169 }
2170 }
2171
2172 let (fraction, termination) = crossings
2173 .into_iter()
2174 .filter(|(fraction, _)| fraction.is_finite() && (0.0..=1.0).contains(fraction))
2175 .min_by(|left, right| {
2176 let priority = |termination: TrajectoryTermination| match termination {
2177 TrajectoryTermination::GroundThreshold => 0,
2178 TrajectoryTermination::MaxRange => 1,
2179 TrajectoryTermination::TimeLimit => 2,
2180 TrajectoryTermination::VelocityFloor => 3,
2181 };
2182 left.0
2183 .total_cmp(&right.0)
2184 .then_with(|| priority(left.1).cmp(&priority(right.1)))
2185 })
2186 .ok_or_else(|| {
2187 BallisticsError::from(
2188 "trajectory integration stopped without crossing a supported boundary",
2189 )
2190 })?;
2191
2192 let mut position = previous.position + (post_position - previous.position) * fraction;
2193 match termination {
2194 TrajectoryTermination::MaxRange => position.x = self.max_range,
2195 TrajectoryTermination::GroundThreshold => {
2196 position.y = self.inputs.ground_threshold;
2197 }
2198 TrajectoryTermination::TimeLimit | TrajectoryTermination::VelocityFloor => {}
2199 }
2200 let velocity_magnitude = previous.velocity_magnitude
2201 + (post_velocity.magnitude() - previous.velocity_magnitude) * fraction;
2202 let mut time = previous.time + (post_time - previous.time) * fraction;
2203 if termination == TrajectoryTermination::TimeLimit {
2204 time = TRAJECTORY_TIME_LIMIT_S;
2205 }
2206 let kinetic_energy =
2207 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2208
2209 if position.y > *max_height {
2210 *max_height = position.y;
2211 }
2212 let terminal_point = TrajectoryPoint {
2213 time,
2214 position,
2215 velocity_magnitude,
2216 kinetic_energy,
2217 drag_coefficient: None,
2218 };
2219 if terminal_point.position.x < previous.position.x {
2220 return Err(BallisticsError::from(
2221 "trajectory terminal state reversed downrange before the crossed boundary",
2222 ));
2223 }
2224 if terminal_point.position.x == previous.position.x {
2225 let last = points.last_mut().ok_or_else(|| {
2230 BallisticsError::from("trajectory points disappeared during terminal finalization")
2231 })?;
2232 *last = terminal_point;
2233 } else {
2234 self.push_trajectory_point(points, terminal_point)?;
2235 }
2236 Ok(termination)
2237 }
2238
2239 fn gravity_acceleration(&self) -> Vector3<f64> {
2240 let theta = self.inputs.shooting_angle;
2241 Vector3::new(
2242 -crate::constants::G_ACCEL_MPS2 * theta.sin(),
2243 -crate::constants::G_ACCEL_MPS2 * theta.cos(),
2244 0.0,
2245 )
2246 }
2247
2248 fn get_wind_at_altitude(&self, altitude_m: f64) -> Vector3<f64> {
2249 let model = match self.inputs.wind_shear_model.as_str() {
2264 "logarithmic" => WindShearModel::Logarithmic,
2265 "power_law" | "powerlaw" | "exponential" => WindShearModel::PowerLaw,
2266 "ekman_spiral" | "ekman" => WindShearModel::EkmanSpiral,
2267 "custom_layers" | "custom" => WindShearModel::CustomLayers,
2268 _ => WindShearModel::PowerLaw,
2269 };
2270 let speed_ratio = crate::wind_shear::boundary_layer_speed_ratio(altitude_m, model);
2271
2272 crate::wind::wind_vector(self.wind.speed, self.wind.direction, 0.0) * speed_ratio
2279 + Vector3::new(0.0, self.wind.vertical_speed, 0.0)
2280 }
2281
2282 pub fn solve(&self) -> Result<TrajectoryResult, BallisticsError> {
2283 self.validate_for_solve()?;
2284 let mut result = if self.inputs.use_rk4 {
2285 if self.inputs.use_adaptive_rk45 {
2286 self.solve_rk45()?
2287 } else {
2288 self.solve_rk4()?
2289 }
2290 } else {
2291 self.solve_euler()?
2292 };
2293 self.apply_spin_drift(&mut result);
2294 self.validate_result_sanity(&result)?;
2295 Ok(result)
2296 }
2297
2298 fn apply_spin_drift(&self, result: &mut TrajectoryResult) {
2304 if !self.inputs.use_enhanced_spin_drift {
2305 return;
2306 }
2307 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 {
2311 return;
2312 }
2313
2314 let sg = self.effective_spin_drift_sg();
2321
2322 for p in result.points.iter_mut() {
2323 if p.time <= 0.0 {
2324 continue;
2325 }
2326 p.position.z +=
2328 crate::spin_drift::litz_drift_meters(sg, p.time, self.inputs.is_twist_right);
2329 }
2330
2331 if let Some(samples) = result.sampled_points.as_mut() {
2335 for s in samples.iter_mut() {
2336 if s.time_s <= 0.0 {
2337 continue;
2338 }
2339 s.wind_drift_m +=
2340 crate::spin_drift::litz_drift_meters(sg, s.time_s, self.inputs.is_twist_right);
2341 }
2342 }
2343 }
2344
2345 fn effective_spin_drift_sg(&self) -> f64 {
2350 let (_, _, temp_c, press_hpa) = self.resolved_atmosphere();
2351 crate::spin_drift::effective_sg_from_inputs(&self.inputs, temp_c, press_hpa)
2352 }
2353
2354 fn initial_position(&self) -> Vector3<f64> {
2372 if self.inputs.cant_angle == 0.0 && self.inputs.sight_offset_lateral_m == 0.0 {
2373 return Vector3::new(0.0, self.inputs.muzzle_height, 0.0);
2374 }
2375 let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
2376 let sh = self.inputs.sight_height;
2377 let off = self.inputs.sight_offset_lateral_m;
2378 Vector3::new(
2379 0.0,
2380 self.inputs.muzzle_height + sh * (1.0 - cos_c) + off * sin_c,
2381 -sh * sin_c - off * cos_c,
2382 )
2383 }
2384
2385 fn build_sampled_points(
2392 &self,
2393 points: &[TrajectoryPoint],
2394 max_height: f64,
2395 transonic_distances: Vec<f64>,
2396 mach_transitions: &MachTransitionTracker,
2397 ) -> Result<Option<Vec<TrajectorySample>>, BallisticsError> {
2398 if !self.inputs.enable_trajectory_sampling {
2399 return Ok(None);
2400 }
2401
2402 let last_point = points.last().ok_or("No trajectory points generated")?;
2403 let trajectory_data = TrajectoryData {
2404 times: points.iter().map(|p| p.time).collect(),
2405 positions: points.iter().map(|p| p.position).collect(),
2406 velocities: points
2407 .iter()
2408 .map(|p| {
2409 Vector3::new(0.0, 0.0, p.velocity_magnitude)
2411 })
2412 .collect(),
2413 transonic_distances, mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2415 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2416 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2417 };
2418
2419 let sight_position_m = self.inputs.muzzle_height + self.inputs.sight_height;
2424 let target_reference = self.inputs.drops_reference == DropsReference::Target;
2425 let target_vertical_height_m = if target_reference && self.inputs.target_height != 0.0 {
2432 self.inputs.target_height
2433 } else {
2434 sight_position_m
2435 };
2436 let outputs = TrajectoryOutputs {
2437 target_distance_horiz_m: last_point.position.x, target_vertical_height_m,
2439 time_of_flight_s: last_point.time,
2440 max_ord_dist_horiz_m: max_height,
2441 sight_height_m: sight_position_m,
2442 };
2443
2444 let mut samples = sample_trajectory(
2446 &trajectory_data,
2447 &outputs,
2448 self.inputs.sample_interval,
2449 self.inputs.bullet_mass,
2450 )?;
2451
2452 if target_reference {
2459 let cos_theta = self.inputs.shooting_angle.cos();
2460 for sample in &mut samples {
2461 sample.drop_m /= cos_theta;
2462 }
2463 }
2464 Ok(Some(samples))
2465 }
2466
2467 fn solve_euler(&self) -> Result<TrajectoryResult, BallisticsError> {
2468 let mut time = 0.0;
2470 let mut position = self.initial_position();
2474 let aj_components = self.aerodynamic_jump_components();
2480 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2481 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2482 let mut velocity = Vector3::new(
2483 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2487
2488 let mut points = Vec::new();
2489 let mut max_height = position.y;
2490 let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
2496 let mut mach_transitions = MachTransitionTracker::default();
2497
2498 let mut angular_state = if self.inputs.enable_precession_nutation {
2500 Some(AngularState {
2501 pitch_angle: 0.001, yaw_angle: 0.001,
2503 pitch_rate: 0.0,
2504 yaw_rate: 0.0,
2505 precession_angle: 0.0,
2506 nutation_phase: 0.0,
2507 })
2508 } else {
2509 None
2510 };
2511 let mut max_yaw_angle = 0.0;
2512 let mut max_precession_angle = 0.0;
2513
2514 let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2516 self.resolved_atmosphere();
2517 let base_ratio = air_density / 1.225;
2522
2523 let wind_vector = crate::wind::wind_vector(
2529 self.wind.speed,
2530 self.wind.direction,
2531 self.wind.vertical_speed,
2532 );
2533
2534 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2537 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2538 );
2539
2540 while position.x < self.max_range
2542 && position.y > self.inputs.ground_threshold
2543 && time < TRAJECTORY_TIME_LIMIT_S
2544 {
2545 let velocity_magnitude = velocity.magnitude();
2547 let kinetic_energy =
2548 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2549
2550 self.push_trajectory_point(
2551 &mut points,
2552 TrajectoryPoint {
2553 time,
2554 position,
2555 velocity_magnitude,
2556 kinetic_energy,
2557 drag_coefficient: None,
2558 },
2559 )?;
2560
2561 {
2564 let mach_here = if speed_of_sound > 0.0 {
2565 velocity_magnitude / speed_of_sound
2566 } else {
2567 0.0
2568 };
2569 mach_transitions.record_downward_crossings(
2570 mach_here,
2571 position.x,
2572 &mut transonic_distances,
2573 );
2574 }
2575
2576 if position.y > max_height {
2578 max_height = position.y;
2579 }
2580
2581 if self.inputs.enable_pitch_damping {
2583 let mach = velocity_magnitude / speed_of_sound;
2584
2585 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
2587 transonic_mach = Some(mach);
2588 }
2589
2590 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
2592
2593 if pitch_damping < min_pitch_damping {
2595 min_pitch_damping = pitch_damping;
2596 }
2597 }
2598
2599 if self.inputs.enable_precession_nutation {
2601 if let Some(ref mut state) = angular_state {
2602 let velocity_magnitude = velocity.magnitude();
2603 let params = self.precession_nutation_params(
2604 velocity_magnitude,
2605 air_density,
2606 speed_of_sound,
2607 );
2608
2609 *state = calculate_combined_angular_motion(
2611 ¶ms,
2612 state,
2613 time,
2614 self.time_step,
2615 0.001, );
2617
2618 if state.yaw_angle.abs() > max_yaw_angle {
2620 max_yaw_angle = state.yaw_angle.abs();
2621 }
2622 if state.precession_angle.abs() > max_precession_angle {
2623 max_precession_angle = state.precession_angle.abs();
2624 }
2625 }
2626 }
2627
2628 let acceleration = self.calculate_acceleration(
2635 &position,
2636 &velocity,
2637 &wind_vector,
2638 (resolved_temp_c, resolved_press_hpa, base_ratio),
2639 );
2640
2641 velocity += acceleration * self.time_step;
2643 position += velocity * self.time_step;
2644 time += self.time_step;
2645 self.validate_integration_state(&position, &velocity, time)?;
2646 }
2647
2648 let termination =
2649 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
2650
2651 self.annotate_drag_coefficients(&mut points, speed_of_sound);
2656
2657 let last_point = points.last().ok_or("No trajectory points generated")?;
2658
2659 let sampled_points =
2661 self.build_sampled_points(&points, max_height, transonic_distances, &mach_transitions)?;
2662
2663 Ok(TrajectoryResult {
2664 max_range: last_point.position.x, max_height,
2666 time_of_flight: last_point.time,
2667 impact_velocity: last_point.velocity_magnitude,
2668 impact_energy: last_point.kinetic_energy,
2669 projectile_mass_kg: self.inputs.bullet_mass,
2670 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
2671 station_speed_of_sound_mps: speed_of_sound,
2672 termination,
2673 points,
2674 sampled_points,
2675 min_pitch_damping: if self.inputs.enable_pitch_damping {
2676 Some(min_pitch_damping)
2677 } else {
2678 None
2679 },
2680 transonic_mach,
2681 angular_state,
2682 max_yaw_angle: if self.inputs.enable_precession_nutation {
2683 Some(max_yaw_angle)
2684 } else {
2685 None
2686 },
2687 max_precession_angle: if self.inputs.enable_precession_nutation {
2688 Some(max_precession_angle)
2689 } else {
2690 None
2691 },
2692 aerodynamic_jump: aj_components,
2693 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2694 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2695 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2696 })
2697 }
2698
2699 fn solve_rk4(&self) -> Result<TrajectoryResult, BallisticsError> {
2700 let mut time = 0.0;
2702 let mut position = self.initial_position();
2707
2708 let aj_components = self.aerodynamic_jump_components();
2714 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2715 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2716 let mut velocity = Vector3::new(
2717 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2721
2722 let mut points = Vec::new();
2723 let mut max_height = position.y;
2724 let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
2730 let mut mach_transitions = MachTransitionTracker::default();
2731
2732 let mut angular_state = if self.inputs.enable_precession_nutation {
2734 Some(AngularState {
2735 pitch_angle: 0.001, yaw_angle: 0.001,
2737 pitch_rate: 0.0,
2738 yaw_rate: 0.0,
2739 precession_angle: 0.0,
2740 nutation_phase: 0.0,
2741 })
2742 } else {
2743 None
2744 };
2745 let mut max_yaw_angle = 0.0;
2746 let mut max_precession_angle = 0.0;
2747
2748 let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2750 self.resolved_atmosphere();
2751 let base_ratio = air_density / 1.225;
2756
2757 let wind_vector = crate::wind::wind_vector(
2763 self.wind.speed,
2764 self.wind.direction,
2765 self.wind.vertical_speed,
2766 );
2767
2768 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2771 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2772 );
2773
2774 while position.x < self.max_range
2776 && position.y > self.inputs.ground_threshold
2777 && time < TRAJECTORY_TIME_LIMIT_S
2778 {
2779 let velocity_magnitude = velocity.magnitude();
2781 let kinetic_energy =
2782 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2783
2784 self.push_trajectory_point(
2785 &mut points,
2786 TrajectoryPoint {
2787 time,
2788 position,
2789 velocity_magnitude,
2790 kinetic_energy,
2791 drag_coefficient: None,
2792 },
2793 )?;
2794
2795 {
2798 let mach_here = if speed_of_sound > 0.0 {
2799 velocity_magnitude / speed_of_sound
2800 } else {
2801 0.0
2802 };
2803 mach_transitions.record_downward_crossings(
2804 mach_here,
2805 position.x,
2806 &mut transonic_distances,
2807 );
2808 }
2809
2810 if position.y > max_height {
2811 max_height = position.y;
2812 }
2813
2814 if self.inputs.enable_pitch_damping {
2816 let mach = velocity_magnitude / speed_of_sound;
2817
2818 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
2820 transonic_mach = Some(mach);
2821 }
2822
2823 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
2825
2826 if pitch_damping < min_pitch_damping {
2828 min_pitch_damping = pitch_damping;
2829 }
2830 }
2831
2832 if self.inputs.enable_precession_nutation {
2834 if let Some(ref mut state) = angular_state {
2835 let velocity_magnitude = velocity.magnitude();
2836 let params = self.precession_nutation_params(
2837 velocity_magnitude,
2838 air_density,
2839 speed_of_sound,
2840 );
2841
2842 *state = calculate_combined_angular_motion(
2844 ¶ms,
2845 state,
2846 time,
2847 self.time_step,
2848 0.001, );
2850
2851 if state.yaw_angle.abs() > max_yaw_angle {
2853 max_yaw_angle = state.yaw_angle.abs();
2854 }
2855 if state.precession_angle.abs() > max_precession_angle {
2856 max_precession_angle = state.precession_angle.abs();
2857 }
2858 }
2859 }
2860
2861 let dt = self.time_step;
2863
2864 let acc1 = self.calculate_acceleration(
2866 &position,
2867 &velocity,
2868 &wind_vector,
2869 (resolved_temp_c, resolved_press_hpa, base_ratio),
2870 );
2871
2872 let pos2 = position + velocity * (dt * 0.5);
2874 let vel2 = velocity + acc1 * (dt * 0.5);
2875 let acc2 = self.calculate_acceleration(
2876 &pos2,
2877 &vel2,
2878 &wind_vector,
2879 (resolved_temp_c, resolved_press_hpa, base_ratio),
2880 );
2881
2882 let pos3 = position + vel2 * (dt * 0.5);
2884 let vel3 = velocity + acc2 * (dt * 0.5);
2885 let acc3 = self.calculate_acceleration(
2886 &pos3,
2887 &vel3,
2888 &wind_vector,
2889 (resolved_temp_c, resolved_press_hpa, base_ratio),
2890 );
2891
2892 let pos4 = position + vel3 * dt;
2894 let vel4 = velocity + acc3 * dt;
2895 let acc4 = self.calculate_acceleration(
2896 &pos4,
2897 &vel4,
2898 &wind_vector,
2899 (resolved_temp_c, resolved_press_hpa, base_ratio),
2900 );
2901
2902 position += (velocity + vel2 * 2.0 + vel3 * 2.0 + vel4) * (dt / 6.0);
2904 velocity += (acc1 + acc2 * 2.0 + acc3 * 2.0 + acc4) * (dt / 6.0);
2905 time += dt;
2906 self.validate_integration_state(&position, &velocity, time)?;
2907 }
2908
2909 let termination =
2910 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
2911
2912 self.annotate_drag_coefficients(&mut points, speed_of_sound);
2917
2918 let last_point = points.last().ok_or("No trajectory points generated")?;
2919
2920 let sampled_points =
2922 self.build_sampled_points(&points, max_height, transonic_distances, &mach_transitions)?;
2923
2924 Ok(TrajectoryResult {
2925 max_range: last_point.position.x, max_height,
2927 time_of_flight: last_point.time,
2928 impact_velocity: last_point.velocity_magnitude,
2929 impact_energy: last_point.kinetic_energy,
2930 projectile_mass_kg: self.inputs.bullet_mass,
2931 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
2932 station_speed_of_sound_mps: speed_of_sound,
2933 termination,
2934 points,
2935 sampled_points,
2936 min_pitch_damping: if self.inputs.enable_pitch_damping {
2937 Some(min_pitch_damping)
2938 } else {
2939 None
2940 },
2941 transonic_mach,
2942 angular_state,
2943 max_yaw_angle: if self.inputs.enable_precession_nutation {
2944 Some(max_yaw_angle)
2945 } else {
2946 None
2947 },
2948 max_precession_angle: if self.inputs.enable_precession_nutation {
2949 Some(max_precession_angle)
2950 } else {
2951 None
2952 },
2953 aerodynamic_jump: aj_components,
2954 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2955 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2956 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2957 })
2958 }
2959
2960 fn solve_rk45(&self) -> Result<TrajectoryResult, BallisticsError> {
2961 let mut time = 0.0;
2963 let mut position = self.initial_position();
2967
2968 let aj_components = self.aerodynamic_jump_components();
2974 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2975 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2976 let mut velocity = Vector3::new(
2977 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2981
2982 let mut points = Vec::new();
2983 let mut max_height = position.y;
2984 let mut dt = 0.001; let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2989 self.resolved_atmosphere();
2990 let base_ratio = air_density / 1.225;
2995 let wind_vector = crate::wind::wind_vector(
3000 self.wind.speed,
3001 self.wind.direction,
3002 self.wind.vertical_speed,
3003 );
3004
3005 let mut transonic_distances: Vec<f64> = Vec::new();
3007 let mut mach_transitions = MachTransitionTracker::default();
3008
3009 let mut min_pitch_damping = f64::INFINITY;
3014 let mut transonic_mach: Option<f64> = None;
3015 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
3016 self.inputs.bullet_model.as_deref().unwrap_or("default"),
3017 );
3018 let mut angular_state = if self.inputs.enable_precession_nutation {
3019 Some(AngularState {
3020 pitch_angle: 0.001,
3021 yaw_angle: 0.001,
3022 pitch_rate: 0.0,
3023 yaw_rate: 0.0,
3024 precession_angle: 0.0,
3025 nutation_phase: 0.0,
3026 })
3027 } else {
3028 None
3029 };
3030 let mut max_yaw_angle = 0.0;
3031 let mut max_precession_angle = 0.0;
3032
3033 while position.x < self.max_range
3034 && position.y > self.inputs.ground_threshold
3035 && time < TRAJECTORY_TIME_LIMIT_S
3036 {
3037 let velocity_magnitude = velocity.magnitude();
3039 let kinetic_energy = 0.5 * self.inputs.bullet_mass * velocity_magnitude.powi(2);
3040
3041 self.push_trajectory_point(
3042 &mut points,
3043 TrajectoryPoint {
3044 time,
3045 position,
3046 velocity_magnitude,
3047 kinetic_energy,
3048 drag_coefficient: None,
3049 },
3050 )?;
3051
3052 {
3055 let mach_here = if speed_of_sound > 0.0 {
3056 velocity_magnitude / speed_of_sound
3057 } else {
3058 0.0
3059 };
3060 mach_transitions.record_downward_crossings(
3061 mach_here,
3062 position.x,
3063 &mut transonic_distances,
3064 );
3065 }
3066
3067 if position.y > max_height {
3068 max_height = position.y;
3069 }
3070
3071 if self.inputs.enable_pitch_damping {
3074 let mach = velocity_magnitude / speed_of_sound;
3075 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
3076 transonic_mach = Some(mach);
3077 }
3078 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
3079 if pitch_damping < min_pitch_damping {
3080 min_pitch_damping = pitch_damping;
3081 }
3082 }
3083
3084 let accepted_step = self.adaptive_rk45_step(
3087 &position,
3088 &velocity,
3089 dt,
3090 &wind_vector,
3091 (resolved_temp_c, resolved_press_hpa, base_ratio),
3092 );
3093 debug_assert!(
3094 accepted_step.error <= RK45_TOLERANCE || accepted_step.used_dt <= RK45_MIN_DT
3095 );
3096
3097 if self.inputs.enable_precession_nutation {
3101 if let Some(ref mut state) = angular_state {
3102 let params = self.precession_nutation_params(
3103 velocity_magnitude,
3104 air_density,
3105 speed_of_sound,
3106 );
3107
3108 *state = calculate_combined_angular_motion(
3109 ¶ms,
3110 state,
3111 time,
3112 accepted_step.used_dt,
3113 0.001,
3114 );
3115
3116 if state.yaw_angle.abs() > max_yaw_angle {
3117 max_yaw_angle = state.yaw_angle.abs();
3118 }
3119 if state.precession_angle.abs() > max_precession_angle {
3120 max_precession_angle = state.precession_angle.abs();
3121 }
3122 }
3123 }
3124
3125 position = accepted_step.position;
3126 velocity = accepted_step.velocity;
3127 time += accepted_step.used_dt;
3128 self.validate_integration_state(&position, &velocity, time)?;
3129
3130 dt = accepted_step.next_dt;
3132 }
3133
3134 if points.is_empty() {
3136 return Err(BallisticsError::from("No trajectory points calculated"));
3137 }
3138
3139 let termination =
3141 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
3142
3143 self.annotate_drag_coefficients(&mut points, speed_of_sound);
3145
3146 let last_point = points.last().unwrap();
3147
3148 let sampled_points =
3150 self.build_sampled_points(&points, max_height, transonic_distances, &mach_transitions)?;
3151
3152 Ok(TrajectoryResult {
3153 max_range: last_point.position.x, max_height,
3155 time_of_flight: last_point.time,
3156 impact_velocity: last_point.velocity_magnitude,
3157 impact_energy: last_point.kinetic_energy,
3158 projectile_mass_kg: self.inputs.bullet_mass,
3159 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
3160 station_speed_of_sound_mps: speed_of_sound,
3161 termination,
3162 points,
3163 sampled_points,
3164 min_pitch_damping: if self.inputs.enable_pitch_damping {
3165 Some(min_pitch_damping)
3166 } else {
3167 None
3168 },
3169 transonic_mach,
3170 angular_state,
3171 max_yaw_angle: if self.inputs.enable_precession_nutation {
3172 Some(max_yaw_angle)
3173 } else {
3174 None
3175 },
3176 max_precession_angle: if self.inputs.enable_precession_nutation {
3177 Some(max_precession_angle)
3178 } else {
3179 None
3180 },
3181 aerodynamic_jump: aj_components,
3182 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
3183 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
3184 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
3185 })
3186 }
3187
3188 fn adaptive_rk45_step(
3189 &self,
3190 position: &Vector3<f64>,
3191 velocity: &Vector3<f64>,
3192 initial_dt: f64,
3193 wind_vector: &Vector3<f64>,
3194 resolved_atmo: (f64, f64, f64),
3195 ) -> Rk45AcceptedStep {
3196 let mut trial_dt = initial_dt;
3197
3198 loop {
3199 let trial = self.rk45_step(
3200 position,
3201 velocity,
3202 trial_dt,
3203 wind_vector,
3204 RK45_TOLERANCE,
3205 resolved_atmo,
3206 );
3207 let next_dt = if trial.suggested_dt.is_finite() {
3212 (RK45_SAFETY_FACTOR * trial.suggested_dt).clamp(RK45_MIN_DT, RK45_MAX_DT)
3213 } else {
3214 RK45_MIN_DT
3215 };
3216
3217 if trial.error <= RK45_TOLERANCE || trial_dt <= RK45_MIN_DT {
3218 return Rk45AcceptedStep {
3219 position: trial.position,
3220 velocity: trial.velocity,
3221 used_dt: trial_dt,
3222 next_dt,
3223 error: trial.error,
3224 };
3225 }
3226
3227 trial_dt = next_dt;
3228 }
3229 }
3230
3231 fn rk45_step(
3232 &self,
3233 position: &Vector3<f64>,
3234 velocity: &Vector3<f64>,
3235 dt: f64,
3236 wind_vector: &Vector3<f64>,
3237 tolerance: f64,
3238 resolved_atmo: (f64, f64, f64), ) -> Rk45Trial {
3240 const A21: f64 = 1.0 / 5.0;
3242 const A31: f64 = 3.0 / 40.0;
3243 const A32: f64 = 9.0 / 40.0;
3244 const A41: f64 = 44.0 / 45.0;
3245 const A42: f64 = -56.0 / 15.0;
3246 const A43: f64 = 32.0 / 9.0;
3247 const A51: f64 = 19372.0 / 6561.0;
3248 const A52: f64 = -25360.0 / 2187.0;
3249 const A53: f64 = 64448.0 / 6561.0;
3250 const A54: f64 = -212.0 / 729.0;
3251 const A61: f64 = 9017.0 / 3168.0;
3252 const A62: f64 = -355.0 / 33.0;
3253 const A63: f64 = 46732.0 / 5247.0;
3254 const A64: f64 = 49.0 / 176.0;
3255 const A65: f64 = -5103.0 / 18656.0;
3256 const A71: f64 = 35.0 / 384.0;
3257 const A73: f64 = 500.0 / 1113.0;
3258 const A74: f64 = 125.0 / 192.0;
3259 const A75: f64 = -2187.0 / 6784.0;
3260 const A76: f64 = 11.0 / 84.0;
3261
3262 const B1: f64 = 35.0 / 384.0;
3264 const B3: f64 = 500.0 / 1113.0;
3265 const B4: f64 = 125.0 / 192.0;
3266 const B5: f64 = -2187.0 / 6784.0;
3267 const B6: f64 = 11.0 / 84.0;
3268
3269 const B1_ERR: f64 = 5179.0 / 57600.0;
3271 const B3_ERR: f64 = 7571.0 / 16695.0;
3272 const B4_ERR: f64 = 393.0 / 640.0;
3273 const B5_ERR: f64 = -92097.0 / 339200.0;
3274 const B6_ERR: f64 = 187.0 / 2100.0;
3275 const B7_ERR: f64 = 1.0 / 40.0;
3276
3277 let k1_v = self.calculate_acceleration(position, velocity, wind_vector, resolved_atmo);
3279 let k1_p = *velocity;
3280
3281 let p2 = position + dt * A21 * k1_p;
3282 let v2 = velocity + dt * A21 * k1_v;
3283 let k2_v = self.calculate_acceleration(&p2, &v2, wind_vector, resolved_atmo);
3284 let k2_p = v2;
3285
3286 let p3 = position + dt * (A31 * k1_p + A32 * k2_p);
3287 let v3 = velocity + dt * (A31 * k1_v + A32 * k2_v);
3288 let k3_v = self.calculate_acceleration(&p3, &v3, wind_vector, resolved_atmo);
3289 let k3_p = v3;
3290
3291 let p4 = position + dt * (A41 * k1_p + A42 * k2_p + A43 * k3_p);
3292 let v4 = velocity + dt * (A41 * k1_v + A42 * k2_v + A43 * k3_v);
3293 let k4_v = self.calculate_acceleration(&p4, &v4, wind_vector, resolved_atmo);
3294 let k4_p = v4;
3295
3296 let p5 = position + dt * (A51 * k1_p + A52 * k2_p + A53 * k3_p + A54 * k4_p);
3297 let v5 = velocity + dt * (A51 * k1_v + A52 * k2_v + A53 * k3_v + A54 * k4_v);
3298 let k5_v = self.calculate_acceleration(&p5, &v5, wind_vector, resolved_atmo);
3299 let k5_p = v5;
3300
3301 let p6 = position + dt * (A61 * k1_p + A62 * k2_p + A63 * k3_p + A64 * k4_p + A65 * k5_p);
3302 let v6 = velocity + dt * (A61 * k1_v + A62 * k2_v + A63 * k3_v + A64 * k4_v + A65 * k5_v);
3303 let k6_v = self.calculate_acceleration(&p6, &v6, wind_vector, resolved_atmo);
3304 let k6_p = v6;
3305
3306 let p7 = position + dt * (A71 * k1_p + A73 * k3_p + A74 * k4_p + A75 * k5_p + A76 * k6_p);
3307 let v7 = velocity + dt * (A71 * k1_v + A73 * k3_v + A74 * k4_v + A75 * k5_v + A76 * k6_v);
3308 let k7_v = self.calculate_acceleration(&p7, &v7, wind_vector, resolved_atmo);
3309 let k7_p = v7;
3310
3311 let new_pos = position + dt * (B1 * k1_p + B3 * k3_p + B4 * k4_p + B5 * k5_p + B6 * k6_p);
3313 let new_vel = velocity + dt * (B1 * k1_v + B3 * k3_v + B4 * k4_v + B5 * k5_v + B6 * k6_v);
3314
3315 let pos_err = position
3317 + dt * (B1_ERR * k1_p
3318 + B3_ERR * k3_p
3319 + B4_ERR * k4_p
3320 + B5_ERR * k5_p
3321 + B6_ERR * k6_p
3322 + B7_ERR * k7_p);
3323 let vel_err = velocity
3324 + dt * (B1_ERR * k1_v
3325 + B3_ERR * k3_v
3326 + B4_ERR * k4_v
3327 + B5_ERR * k5_v
3328 + B6_ERR * k6_v
3329 + B7_ERR * k7_v);
3330
3331 let error = cli_rk45_error_norm(position, velocity, &new_pos, &new_vel, &pos_err, &vel_err);
3333
3334 let dt_new = if error < tolerance {
3336 dt * (tolerance / error).powf(0.2).min(2.0)
3337 } else {
3338 dt * (tolerance / error).powf(0.25).max(0.1)
3339 };
3340
3341 Rk45Trial {
3342 position: new_pos,
3343 velocity: new_vel,
3344 suggested_dt: dt_new,
3345 error,
3346 }
3347 }
3348
3349 fn apply_cluster_bc_correction(&self, base_bc: f64, velocity_fps: f64) -> f64 {
3350 if let Some(ref cluster_bc) = self.cluster_bc {
3351 cluster_bc.apply_correction_for_drag_model(
3352 base_bc,
3353 self.inputs.caliber_inches,
3354 self.inputs.weight_grains,
3355 velocity_fps,
3356 self.inputs.bc_type,
3357 )
3358 } else {
3359 base_bc
3360 }
3361 }
3362
3363 fn calculate_acceleration(
3364 &self,
3365 position: &Vector3<f64>,
3366 velocity: &Vector3<f64>,
3367 wind_vector: &Vector3<f64>,
3368 resolved_atmo: (f64, f64, f64), ) -> Vector3<f64> {
3370 let actual_wind = if let Some(ref sock) = self.wind_sock {
3376 sock.vector_for_range_stateless(position.x)
3377 } else if self.inputs.enable_wind_shear {
3378 self.get_wind_at_altitude(position.y)
3379 } else {
3380 *wind_vector
3381 };
3382 let actual_wind =
3383 crate::derivatives::level_vector_to_shot_frame(actual_wind, self.inputs.shooting_angle);
3384
3385 let relative_velocity = velocity - actual_wind;
3386 let velocity_magnitude = relative_velocity.magnitude();
3387
3388 if velocity_magnitude < 0.001 {
3389 return self.gravity_acceleration();
3390 }
3391
3392 let (base_temp_c, base_press_hpa, station_ratio) = resolved_atmo;
3403
3404 let (drag_base_temp_c, drag_base_press_hpa, drag_base_ratio, drag_humidity_percent) =
3412 if let Some(ref sock) = self.atmo_sock {
3413 let (zone_temp_c, zone_press_hpa, zone_humidity) = sock.atmo_for_range(position.x);
3414 let zone_base_ratio = crate::atmosphere::calculate_air_density_cimp(
3415 zone_temp_c,
3416 zone_press_hpa,
3417 zone_humidity,
3418 ) / 1.225;
3419 (zone_temp_c, zone_press_hpa, zone_base_ratio, zone_humidity)
3420 } else {
3421 (
3422 base_temp_c,
3423 base_press_hpa,
3424 station_ratio,
3425 self.atmosphere.humidity,
3426 )
3427 };
3428 let local_alt = crate::atmosphere::shot_frame_altitude(
3429 self.atmosphere.altitude,
3430 position.x,
3431 position.y,
3432 self.inputs.shooting_angle,
3433 );
3434 let (air_density, speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
3435 local_alt,
3436 self.atmosphere.altitude,
3437 drag_base_temp_c,
3438 drag_base_press_hpa,
3439 drag_base_ratio,
3440 drag_humidity_percent,
3441 );
3442
3443 let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
3447
3448 let velocity_fps = velocity_magnitude * 3.28084;
3450
3451 let cd_to_retard = crate::constants::CD_TO_RETARD;
3456 let standard_factor = cd * cd_to_retard;
3457 let density_scale = air_density / 1.225; let a_drag_ft_s2 =
3461 (velocity_fps * velocity_fps) * standard_factor * density_scale / retard_denom;
3462 let a_drag_m_s2 = a_drag_ft_s2 * 0.3048; let drag_acceleration = -a_drag_m_s2 * (relative_velocity / velocity_magnitude);
3466
3467 let mut accel = drag_acceleration + self.gravity_acceleration();
3470
3471 if self.inputs.enable_coriolis {
3474 if let Some(lat_deg) = self.inputs.latitude {
3475 let omega_earth = 7.2921159e-5_f64; let lat = lat_deg.to_radians();
3477 let az = self.inputs.shot_azimuth; let omega = Vector3::new(
3484 omega_earth * lat.cos() * az.cos(), omega_earth * lat.sin(), -omega_earth * lat.cos() * az.sin(), );
3488 let omega = crate::derivatives::level_vector_to_shot_frame(
3489 omega,
3490 self.inputs.shooting_angle,
3491 );
3492 accel += -2.0 * omega.cross(velocity);
3497 }
3498 }
3499
3500 if self.inputs.enable_magnus
3507 && !self.inputs.use_enhanced_spin_drift
3508 && self.inputs.bullet_diameter > 0.0
3509 && self.inputs.twist_rate > 0.0
3510 {
3511 let diameter_m = self.inputs.bullet_diameter;
3512 let (spin_rad_s, spin_param) = crate::spin_drift::calculate_magnus_spin_state(
3513 self.inputs.muzzle_velocity,
3514 velocity_magnitude,
3515 self.inputs.twist_rate,
3516 diameter_m,
3517 );
3518 let mach = velocity_magnitude / speed_of_sound;
3520
3521 let d_in = self.inputs.bullet_diameter / 0.0254;
3523 let m_gr = self.inputs.bullet_mass / crate::constants::GRAINS_TO_KG;
3524 let l_in = if self.inputs.bullet_length > 0.0 {
3525 self.inputs.bullet_length / 0.0254
3526 } else {
3527 let est_m = crate::stability::estimate_bullet_length_m(
3529 self.inputs.bullet_diameter,
3530 self.inputs.bullet_mass,
3531 );
3532 if est_m > 0.0 {
3533 est_m / 0.0254
3534 } else {
3535 4.5 * d_in
3536 }
3537 };
3538 let sg = crate::spin_drift::calculate_dynamic_stability(
3542 m_gr,
3543 velocity_magnitude,
3544 spin_rad_s,
3545 d_in,
3546 l_in,
3547 air_density,
3548 );
3549
3550 let (yaw_rad, _) = crate::spin_drift::calculate_yaw_of_repose(
3552 sg,
3553 velocity_magnitude,
3554 spin_rad_s,
3555 0.0, 0.0, air_density,
3558 d_in,
3559 l_in,
3560 m_gr,
3561 mach,
3562 "match",
3563 false,
3564 );
3565
3566 let c_np = crate::derivatives::calculate_magnus_moment_coefficient(mach);
3568 let area = std::f64::consts::PI * (diameter_m / 2.0).powi(2);
3569 let magnus_force = 0.5
3570 * air_density
3571 * velocity_magnitude.powi(2)
3572 * area
3573 * c_np
3574 * spin_param
3575 * yaw_rad.sin();
3576
3577 if magnus_force.abs() > 1e-12 {
3581 if let Some(dir) = crate::derivatives::yaw_of_repose_magnus_direction(
3582 relative_velocity,
3583 self.gravity_acceleration(),
3584 self.inputs.is_twist_right,
3585 ) {
3586 accel += (magnus_force / self.inputs.bullet_mass) * dir;
3587 }
3588 }
3589 }
3590
3591 accel
3592 }
3593
3594 fn drag_terms(&self, velocity_magnitude: f64, speed_of_sound: f64) -> (f64, f64) {
3606 let cd = self.calculate_drag_coefficient(velocity_magnitude, speed_of_sound);
3607
3608 let velocity_fps = velocity_magnitude * 3.28084;
3609
3610 let (base_bc, bc_from_segments) = if let Some(segments) = self
3615 .inputs
3616 .bc_segments_data
3617 .as_ref()
3618 .filter(|segments| self.inputs.use_bc_segments && !segments.is_empty())
3619 {
3620 (
3622 crate::bc_estimation::velocity_segment_bc(
3623 velocity_fps,
3624 segments,
3625 self.inputs.bc_value,
3626 ),
3627 true,
3628 )
3629 } else if let Some(segments) = self
3630 .inputs
3631 .bc_segments
3632 .as_ref()
3633 .filter(|segments| !segments.is_empty())
3634 {
3635 (
3636 crate::derivatives::interpolated_bc(
3637 velocity_magnitude / speed_of_sound,
3638 segments,
3639 Some(&self.inputs),
3640 ),
3641 true,
3642 )
3643 } else {
3644 (self.inputs.bc_value, false)
3645 };
3646
3647 let effective_bc = if bc_from_segments {
3652 base_bc
3653 } else {
3654 self.apply_cluster_bc_correction(base_bc, velocity_fps)
3655 };
3656 let effective_bc = effective_bc.max(1e-6);
3659
3660 let retard_denom = if self.inputs.custom_drag_table.is_some() {
3665 self.inputs.custom_drag_denominator(effective_bc)
3666 } else {
3667 effective_bc
3668 };
3669
3670 (cd, retard_denom)
3671 }
3672
3673 pub fn effective_drag_coefficient(
3694 &self,
3695 velocity_magnitude: f64,
3696 speed_of_sound: f64,
3697 ) -> Option<f64> {
3698 if !velocity_magnitude.is_finite() || speed_of_sound <= 1e-9 {
3699 return None;
3700 }
3701 let sectional_density = self.inputs.sectional_density_lb_in2()?;
3702 let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
3703 if retard_denom <= 0.0 {
3704 return None;
3705 }
3706 let effective = cd * sectional_density / retard_denom;
3707 effective.is_finite().then_some(effective)
3708 }
3709
3710 fn annotate_drag_coefficients(&self, points: &mut [TrajectoryPoint], speed_of_sound: f64) {
3720 for point in points.iter_mut() {
3721 point.drag_coefficient =
3722 self.effective_drag_coefficient(point.velocity_magnitude, speed_of_sound);
3723 }
3724 }
3725
3726 fn calculate_drag_coefficient(&self, velocity: f64, speed_of_sound: f64) -> f64 {
3727 let mach = velocity / speed_of_sound;
3728
3729 if let Some(ref table) = self.inputs.custom_drag_table {
3733 return table.interpolate(mach) * self.inputs.cd_scale;
3738 }
3739
3740 crate::drag::get_drag_coefficient(mach, &self.inputs.bc_type)
3743 }
3744}
3745
3746#[derive(Debug, Clone)]
3748pub struct MonteCarloParams {
3749 pub num_simulations: usize,
3750 pub velocity_std_dev: f64,
3751 pub angle_std_dev: f64,
3752 pub bc_std_dev: f64,
3753 pub wind_speed_std_dev: f64,
3754 pub target_distance: Option<f64>,
3755 pub base_wind_speed: f64,
3756 pub base_wind_direction: f64,
3757 pub azimuth_std_dev: f64, }
3759
3760impl Default for MonteCarloParams {
3761 fn default() -> Self {
3762 Self {
3763 num_simulations: 1000,
3764 velocity_std_dev: 1.0,
3765 angle_std_dev: 0.001,
3766 bc_std_dev: 0.01,
3767 wind_speed_std_dev: 1.0,
3768 target_distance: None,
3769 base_wind_speed: 0.0,
3770 base_wind_direction: 0.0,
3771 azimuth_std_dev: 0.001, }
3773 }
3774}
3775
3776#[derive(Debug, Clone)]
3778pub struct MonteCarloResults {
3779 pub ranges: Vec<f64>,
3780 pub impact_velocities: Vec<f64>,
3781 pub impact_positions: Vec<Vector3<f64>>,
3787}
3788
3789pub const DEFAULT_HIT_RADIUS_M: f64 = 0.3;
3792
3793pub const TARGET_NOT_REACHED_SENTINEL_M: f64 = -1.0e9;
3799
3800impl MonteCarloResults {
3801 pub fn position_reached_target(position: &Vector3<f64>) -> bool {
3803 position.iter().all(|component| component.is_finite())
3804 && position.y != TARGET_NOT_REACHED_SENTINEL_M
3805 }
3806
3807 pub fn target_arrival_count(&self) -> usize {
3809 self.impact_positions
3810 .iter()
3811 .filter(|position| Self::position_reached_target(position))
3812 .count()
3813 }
3814
3815 pub fn target_shortfall_fraction(&self) -> f64 {
3818 if self.impact_positions.is_empty() {
3819 return 0.0;
3820 }
3821 (self.impact_positions.len() - self.target_arrival_count()) as f64
3822 / self.impact_positions.len() as f64
3823 }
3824
3825 pub fn target_plane_cep(&self) -> Option<f64> {
3831 let mut radial_misses: Vec<f64> = self
3832 .impact_positions
3833 .iter()
3834 .filter(|position| Self::position_reached_target(position))
3835 .map(Vector3::norm)
3836 .filter(|miss| miss.is_finite())
3837 .collect();
3838 radial_misses.sort_by(f64::total_cmp);
3839 if radial_misses.is_empty() {
3840 None
3841 } else {
3842 Some(radial_misses[radial_misses.len() / 2])
3843 }
3844 }
3845
3846 pub fn hit_probability(&self, hit_radius_m: f64) -> f64 {
3855 if self.impact_positions.is_empty() {
3856 return 0.0;
3857 }
3858 let hits = self
3859 .impact_positions
3860 .iter()
3861 .filter(|position| Self::position_is_hit(position, hit_radius_m))
3862 .count();
3863 hits as f64 / self.impact_positions.len() as f64
3864 }
3865
3866 pub fn position_is_hit(position: &Vector3<f64>, hit_radius_m: f64) -> bool {
3878 Self::position_reached_target(position) && position.norm() < hit_radius_m
3879 }
3880
3881 pub fn hit_probability_wilson(
3903 &self,
3904 hit_radius_m: f64,
3905 level: ConfidenceLevel,
3906 ) -> (f64, (f64, f64), u64) {
3907 let trials = self.impact_positions.len() as u64;
3908 let hits = self
3909 .impact_positions
3910 .iter()
3911 .filter(|position| Self::position_is_hit(position, hit_radius_m))
3912 .count() as u64;
3913 (
3914 self.hit_probability(hit_radius_m),
3915 wilson_interval(hits, trials, level),
3916 trials,
3917 )
3918 }
3919
3920 pub fn rect_hit_probability(&self, width_m: f64, height_m: f64) -> f64 {
3932 let dimensions_invalid =
3933 width_m.is_nan() || width_m <= 0.0 || height_m.is_nan() || height_m <= 0.0;
3934 if self.impact_positions.is_empty() || dimensions_invalid {
3935 return 0.0;
3936 }
3937 let half_width = width_m / 2.0;
3938 let half_height = height_m / 2.0;
3939 let hits = self
3940 .impact_positions
3941 .iter()
3942 .filter(|position| {
3943 Self::position_reached_target(position)
3944 && position.z.abs() <= half_width
3945 && position.y.abs() <= half_height
3946 })
3947 .count();
3948 hits as f64 / self.impact_positions.len() as f64
3949 }
3950}
3951
3952fn wind_from_signed_speed_sample(
3953 signed_speed: f64,
3954 sampled_direction: f64,
3955 vertical_speed: f64,
3956) -> WindConditions {
3957 if signed_speed < 0.0 {
3962 WindConditions {
3963 speed: -signed_speed,
3964 direction: sampled_direction + std::f64::consts::PI,
3965 vertical_speed,
3966 }
3967 } else {
3968 WindConditions {
3969 speed: signed_speed,
3970 direction: sampled_direction,
3971 vertical_speed,
3972 }
3973 }
3974}
3975
3976struct MonteCarloWindSampler {
3977 speed: rand_distr::Normal<f64>,
3978 direction: rand_distr::Normal<f64>,
3979 vertical_speed: f64,
3981}
3982
3983impl MonteCarloWindSampler {
3984 fn new(
3985 base_wind: &WindConditions,
3986 wind_speed_std_dev: f64,
3987 wind_direction_std_dev: f64,
3988 ) -> Result<Self, BallisticsError> {
3989 use rand_distr::Normal;
3990
3991 if !wind_direction_std_dev.is_finite() || wind_direction_std_dev < 0.0 {
3992 return Err("Wind direction standard deviation must be finite and non-negative".into());
3993 }
3994
3995 let speed = Normal::new(base_wind.speed, wind_speed_std_dev)
3996 .map_err(|e| format!("Invalid wind speed distribution: {e}"))?;
3997 let direction = Normal::new(base_wind.direction, wind_direction_std_dev)
3998 .map_err(|e| format!("Invalid wind direction distribution: {e}"))?;
3999 Ok(Self {
4000 speed,
4001 direction,
4002 vertical_speed: base_wind.vertical_speed,
4003 })
4004 }
4005
4006 fn sample<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> WindConditions {
4007 use rand_distr::Distribution;
4008
4009 wind_from_signed_speed_sample(
4010 self.speed.sample(rng),
4011 self.direction.sample(rng),
4012 self.vertical_speed,
4013 )
4014 }
4015}
4016
4017#[derive(Debug, Clone, Copy)]
4020struct TrialOutcome {
4021 range: f64,
4024 impact_velocity: f64,
4026 impact_position: Vector3<f64>,
4029}
4030
4031struct MonteCarloTrialSampler {
4051 base_inputs: BallisticInputs,
4052 atmosphere: AtmosphericConditions,
4053 solver_max_range: f64,
4054 target_distance: f64,
4057 baseline_at_target: Vector3<f64>,
4060 velocity_delta_dist: rand_distr::Normal<f64>,
4061 angle_dist: rand_distr::Normal<f64>,
4062 bc_dist: rand_distr::Normal<f64>,
4063 wind_sampler: MonteCarloWindSampler,
4064 azimuth_dist: rand_distr::Normal<f64>,
4065}
4066
4067impl MonteCarloTrialSampler {
4068 fn new(
4072 base_inputs: BallisticInputs,
4073 base_wind: &WindConditions,
4074 params: &MonteCarloParams,
4075 wind_direction_std_dev: f64,
4076 ) -> Result<Self, BallisticsError> {
4077 use rand_distr::Normal;
4078
4079 let atmosphere = AtmosphericConditions {
4080 temperature: base_inputs.temperature,
4081 pressure: base_inputs.pressure,
4082 humidity: base_inputs.humidity_percent(),
4083 altitude: base_inputs.altitude,
4084 };
4085 let target_hint = params
4086 .target_distance
4087 .unwrap_or(base_inputs.target_distance);
4088 let solver_max_range = target_hint.max(1000.0) * 2.0;
4089
4090 let mut baseline_solver =
4092 TrajectorySolver::new(base_inputs.clone(), base_wind.clone(), atmosphere.clone());
4093 baseline_solver.set_max_range(solver_max_range);
4094 let baseline_result = baseline_solver.solve()?;
4095
4096 let target_distance = params.target_distance.unwrap_or(baseline_result.max_range);
4098
4099 let baseline_at_target = baseline_result
4101 .position_at_range(target_distance)
4102 .ok_or("Could not interpolate baseline at target distance")?;
4103
4104 let velocity_delta_dist = Normal::new(0.0, params.velocity_std_dev)
4109 .map_err(|e| format!("Invalid velocity distribution: {}", e))?;
4110 let angle_dist = Normal::new(base_inputs.muzzle_angle, params.angle_std_dev)
4111 .map_err(|e| format!("Invalid angle distribution: {}", e))?;
4112 let bc_dist = Normal::new(base_inputs.bc_value, params.bc_std_dev)
4113 .map_err(|e| format!("Invalid BC distribution: {}", e))?;
4114 let wind_sampler = MonteCarloWindSampler::new(
4117 base_wind,
4118 params.wind_speed_std_dev,
4119 wind_direction_std_dev,
4120 )?;
4121 let azimuth_dist = Normal::new(base_inputs.azimuth_angle, params.azimuth_std_dev)
4122 .map_err(|e| format!("Invalid azimuth distribution: {}", e))?;
4123
4124 Ok(Self {
4125 base_inputs,
4126 atmosphere,
4127 solver_max_range,
4128 target_distance,
4129 baseline_at_target,
4130 velocity_delta_dist,
4131 angle_dist,
4132 bc_dist,
4133 wind_sampler,
4134 azimuth_dist,
4135 })
4136 }
4137
4138 fn sample_one_trial<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> Option<TrialOutcome> {
4148 use rand_distr::Distribution;
4149
4150 let mut inputs = self.base_inputs.clone();
4152 let muzzle_velocity_delta = self.velocity_delta_dist.sample(&mut *rng);
4153 inputs.muzzle_angle = self.angle_dist.sample(&mut *rng);
4154 inputs.bc_value = self.bc_dist.sample(&mut *rng).max(0.01);
4155 inputs.azimuth_angle = self.azimuth_dist.sample(&mut *rng); let wind = self.wind_sampler.sample(&mut *rng);
4159
4160 let mut solver = TrajectorySolver::new(inputs, wind, self.atmosphere.clone());
4164 solver.inputs.muzzle_velocity =
4165 (solver.inputs.muzzle_velocity + muzzle_velocity_delta).max(0.0);
4166 solver.set_max_range(self.solver_max_range);
4167 let result = solver.solve().ok()?;
4169
4170 let impact_position = if result.max_range < self.target_distance {
4176 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)
4179 } else {
4180 let pos_at_target = result.position_at_range(self.target_distance)?;
4182 Vector3::new(
4187 0.0,
4188 pos_at_target.y - self.baseline_at_target.y,
4189 pos_at_target.z - self.baseline_at_target.z,
4190 )
4191 };
4192
4193 Some(TrialOutcome {
4194 range: result.max_range,
4195 impact_velocity: result.impact_velocity,
4196 impact_position,
4197 })
4198 }
4199}
4200
4201pub fn run_monte_carlo(
4203 base_inputs: BallisticInputs,
4204 params: MonteCarloParams,
4205) -> Result<MonteCarloResults, BallisticsError> {
4206 run_monte_carlo_with_direction_std_dev(base_inputs, params, 0.0)
4207}
4208
4209pub fn run_monte_carlo_with_direction_std_dev(
4214 base_inputs: BallisticInputs,
4215 params: MonteCarloParams,
4216 wind_direction_std_dev: f64,
4217) -> Result<MonteCarloResults, BallisticsError> {
4218 let base_wind = WindConditions {
4219 speed: params.base_wind_speed,
4220 direction: params.base_wind_direction,
4221 vertical_speed: 0.0,
4222 };
4223 run_monte_carlo_with_wind_and_direction_std_dev(
4224 base_inputs,
4225 base_wind,
4226 params,
4227 wind_direction_std_dev,
4228 )
4229}
4230
4231pub fn run_monte_carlo_with_wind(
4233 base_inputs: BallisticInputs,
4234 base_wind: WindConditions,
4235 params: MonteCarloParams,
4236) -> Result<MonteCarloResults, BallisticsError> {
4237 run_monte_carlo_with_wind_and_direction_std_dev(base_inputs, base_wind, params, 0.0)
4238}
4239
4240pub fn run_monte_carlo_with_wind_and_direction_std_dev(
4245 base_inputs: BallisticInputs,
4246 base_wind: WindConditions,
4247 params: MonteCarloParams,
4248 wind_direction_std_dev: f64,
4249) -> Result<MonteCarloResults, BallisticsError> {
4250 let mut rng = rand::rng();
4251 run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
4252 base_inputs,
4253 base_wind,
4254 params,
4255 wind_direction_std_dev,
4256 &mut rng,
4257 )
4258}
4259
4260pub fn run_monte_carlo_with_wind_and_direction_std_dev_seeded(
4267 base_inputs: BallisticInputs,
4268 base_wind: WindConditions,
4269 params: MonteCarloParams,
4270 wind_direction_std_dev: f64,
4271 seed: u64,
4272) -> Result<MonteCarloResults, BallisticsError> {
4273 use rand::{rngs::StdRng, SeedableRng};
4274 let mut rng = StdRng::seed_from_u64(seed);
4275 run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
4276 base_inputs,
4277 base_wind,
4278 params,
4279 wind_direction_std_dev,
4280 &mut rng,
4281 )
4282}
4283
4284fn run_monte_carlo_with_wind_and_direction_std_dev_using_rng<R: rand::Rng + ?Sized>(
4285 base_inputs: BallisticInputs,
4286 base_wind: WindConditions,
4287 params: MonteCarloParams,
4288 wind_direction_std_dev: f64,
4289 rng: &mut R,
4290) -> Result<MonteCarloResults, BallisticsError> {
4291 let mut ranges = Vec::new();
4292 let mut impact_velocities = Vec::new();
4293 let mut impact_positions = Vec::new();
4294
4295 let sampler =
4296 MonteCarloTrialSampler::new(base_inputs, &base_wind, ¶ms, wind_direction_std_dev)?;
4297
4298 for _ in 0..params.num_simulations {
4299 if let Some(outcome) = sampler.sample_one_trial(rng) {
4303 ranges.push(outcome.range);
4304 impact_velocities.push(outcome.impact_velocity);
4305 impact_positions.push(outcome.impact_position);
4306 }
4307 }
4308
4309 if ranges.is_empty() {
4310 return Err("No successful simulations".into());
4311 }
4312
4313 Ok(MonteCarloResults {
4314 ranges,
4315 impact_velocities,
4316 impact_positions,
4317 })
4318}
4319
4320pub const MC_ADAPTIVE_SCHEMA_VERSION_V1: u32 = 1;
4322
4323pub const MC_ADAPTIVE_METHOD_V1: &str = "anytime_beta_binomial_mixture_cs_v1";
4330
4331pub const MC_ADAPTIVE_ASSUMPTIONS_V1: [&str; 4] = [
4345 "Sampling uncertainty only: intervals cover Monte Carlo sampling error, not model error in the trajectory solver or its inputs.",
4346 "Anytime-valid stopping: the beta-binomial mixture confidence sequence keeps its coverage guarantee despite stopping the moment the target half-width is met.",
4347 "Input dispersions are the independent normal distributions declared in MonteCarloParams; correlations between inputs are not modeled.",
4348 "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.",
4349];
4350
4351#[derive(Debug, Clone)]
4358pub struct McConvergence {
4359 pub level: ConfidenceLevel,
4361 pub target_half_width: f64,
4364 pub min_samples: u64,
4368 pub max_samples: u64,
4371 pub batch_size: u64,
4374}
4375
4376impl Default for McConvergence {
4377 fn default() -> Self {
4378 Self {
4379 level: ConfidenceLevel::P95,
4380 target_half_width: 0.02,
4381 min_samples: 1_000,
4382 max_samples: 100_000,
4383 batch_size: 500,
4384 }
4385 }
4386}
4387
4388impl McConvergence {
4389 pub fn validate(&self) -> Result<(), String> {
4398 if !self.target_half_width.is_finite() || self.target_half_width <= 0.0 {
4399 return Err(format!(
4400 "McConvergence.target_half_width must be a finite value greater than zero (got {})",
4401 self.target_half_width
4402 ));
4403 }
4404 if self.batch_size == 0 {
4405 return Err("McConvergence.batch_size must be greater than zero".to_string());
4406 }
4407 if self.max_samples == 0 {
4408 return Err("McConvergence.max_samples must be greater than zero".to_string());
4409 }
4410 if self.max_samples < self.min_samples {
4411 return Err(format!(
4412 "McConvergence.max_samples ({}) must be at least McConvergence.min_samples ({})",
4413 self.max_samples, self.min_samples
4414 ));
4415 }
4416 Ok(())
4417 }
4418}
4419
4420#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize)]
4422#[serde(rename_all = "snake_case")]
4423pub enum McStopReason {
4424 TargetHalfWidthMet,
4427 MaxSamplesReached,
4430}
4431
4432#[derive(Debug, Clone, serde::Serialize)]
4439pub struct AdaptiveMcReportV1 {
4440 pub schema_version: u32,
4442 pub method: String,
4444 pub assumptions: Vec<String>,
4446 pub confidence_percent: u32,
4448 pub hit_probability: f64,
4451 pub ci_low: f64,
4453 pub ci_high: f64,
4455 pub samples: u64,
4457 pub attempts: u64,
4468 pub arrivals: u64,
4483 pub stop_reason: McStopReason,
4485 pub hit_radius_m: f64,
4487 pub target_distance_m: f64,
4489 pub mean_impact_velocity_mps: f64,
4494 pub std_impact_velocity_mps: f64,
4496 pub mean_drop_at_target_m: f64,
4504 pub std_drop_at_target_m: f64,
4507 pub mean_wind_drift_at_target_m: f64,
4510 pub std_wind_drift_at_target_m: f64,
4513}
4514
4515pub fn run_monte_carlo_adaptive_seeded(
4590 base_inputs: &BallisticInputs,
4591 base_wind: &WindConditions,
4592 params: &MonteCarloParams,
4593 convergence: &McConvergence,
4594 hit_radius_m: f64,
4595 seed: u64,
4596) -> Result<AdaptiveMcReportV1, String> {
4597 use rand::{rngs::StdRng, SeedableRng};
4598
4599 convergence.validate()?;
4600
4601 let sampler = MonteCarloTrialSampler::new(base_inputs.clone(), base_wind, params, 0.0)
4604 .map_err(|e| e.to_string())?;
4605
4606 let mut rng = StdRng::seed_from_u64(seed);
4607 let mut hits_cs = BernoulliConfidenceSequence::new(convergence.level);
4608 let mut impact_velocity = Welford::new();
4609 let mut drop_at_target = Welford::new();
4610 let mut drift_at_target = Welford::new();
4611
4612 let mut attempts: u64 = 0;
4613 let mut stop_reason = McStopReason::MaxSamplesReached;
4614
4615 while attempts < convergence.max_samples {
4616 let batch = convergence
4619 .batch_size
4620 .min(convergence.max_samples - attempts);
4621 let mut batch_hits: u64 = 0;
4622 let mut batch_trials: u64 = 0;
4623
4624 for _ in 0..batch {
4625 attempts += 1;
4626 let Some(outcome) = sampler.sample_one_trial(&mut rng) else {
4627 continue; };
4629 batch_trials += 1;
4630 if MonteCarloResults::position_is_hit(&outcome.impact_position, hit_radius_m) {
4631 batch_hits += 1;
4632 }
4633 if MonteCarloResults::position_reached_target(&outcome.impact_position) {
4637 impact_velocity.push(outcome.impact_velocity);
4638 drop_at_target.push(outcome.impact_position.y);
4639 drift_at_target.push(outcome.impact_position.z);
4640 }
4641 }
4642
4643 hits_cs.update_batch(batch_hits, batch_trials);
4644
4645 if hits_cs.trials() >= convergence.min_samples
4646 && hits_cs.half_width() <= convergence.target_half_width
4647 {
4648 stop_reason = McStopReason::TargetHalfWidthMet;
4649 break;
4650 }
4651 }
4652
4653 let samples = hits_cs.trials();
4654 if samples == 0 {
4655 return Err("No successful simulations".to_string());
4656 }
4657 let (ci_low, ci_high) = hits_cs.bounds();
4658
4659 Ok(AdaptiveMcReportV1 {
4660 schema_version: MC_ADAPTIVE_SCHEMA_VERSION_V1,
4661 method: MC_ADAPTIVE_METHOD_V1.to_string(),
4662 assumptions: MC_ADAPTIVE_ASSUMPTIONS_V1
4663 .iter()
4664 .map(|s| s.to_string())
4665 .collect(),
4666 confidence_percent: convergence.level.as_percent(),
4667 hit_probability: hits_cs.successes() as f64 / samples as f64,
4668 ci_low,
4669 ci_high,
4670 samples,
4671 attempts,
4672 arrivals: drop_at_target.count(),
4676 stop_reason,
4677 hit_radius_m,
4678 target_distance_m: sampler.target_distance,
4679 mean_impact_velocity_mps: impact_velocity.mean(),
4680 std_impact_velocity_mps: impact_velocity.sample_std(),
4681 mean_drop_at_target_m: drop_at_target.mean(),
4682 std_drop_at_target_m: drop_at_target.sample_std(),
4683 mean_wind_drift_at_target_m: drift_at_target.mean(),
4684 std_wind_drift_at_target_m: drift_at_target.sample_std(),
4685 })
4686}
4687
4688pub fn calculate_zero_angle(
4690 inputs: BallisticInputs,
4691 target_distance: f64,
4692 target_height: f64,
4693) -> Result<f64, BallisticsError> {
4694 calculate_zero_angle_with_conditions(
4695 inputs,
4696 target_distance,
4697 target_height,
4698 WindConditions::default(),
4699 AtmosphericConditions::default(),
4700 )
4701}
4702
4703pub fn calculate_zero_angle_with_conditions(
4704 inputs: BallisticInputs,
4705 target_distance: f64,
4706 target_height: f64,
4707 wind: WindConditions,
4708 atmosphere: AtmosphericConditions,
4709) -> Result<f64, BallisticsError> {
4710 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
4711 solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
4712}
4713
4714pub fn calculate_zero_angle_with_resolved_conditions(
4720 inputs: BallisticInputs,
4721 target_distance: f64,
4722 target_height: f64,
4723 wind: WindConditions,
4724 atmosphere: AtmosphericConditions,
4725) -> Result<f64, BallisticsError> {
4726 let mut solver =
4727 TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
4728 solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
4729}
4730
4731pub const ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M: f64 = 2000.0;
4739
4740pub fn calculate_zero_range_from_angle_with_conditions(
4755 inputs: BallisticInputs,
4756 zero_angle_rad: f64,
4757 target_height: f64,
4758 wind: WindConditions,
4759 atmosphere: AtmosphericConditions,
4760) -> Result<ZeroCrossings, BallisticsError> {
4761 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
4762 solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
4763 solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
4764}
4765
4766pub fn calculate_zero_range_from_angle_with_resolved_conditions(
4772 inputs: BallisticInputs,
4773 zero_angle_rad: f64,
4774 target_height: f64,
4775 wind: WindConditions,
4776 atmosphere: AtmosphericConditions,
4777) -> Result<ZeroCrossings, BallisticsError> {
4778 let mut solver =
4779 TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
4780 solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
4781 solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
4782}
4783
4784#[derive(Debug, Clone, Copy, PartialEq, Eq)]
4786pub enum BcFitMode {
4787 Drop,
4789 Velocity,
4792}
4793
4794#[derive(Debug, Clone, Copy)]
4796pub struct BcEstimate {
4797 pub bc: f64,
4799 pub rms_error: f64,
4801 pub drag_model: DragModel,
4803 pub mode: BcFitMode,
4805 pub at_bound: bool,
4809}
4810
4811fn fit_value_at(
4819 points: &[TrajectoryPoint],
4820 target_dist: f64,
4821 mode: BcFitMode,
4822 drop_offset: f64,
4823) -> Option<f64> {
4824 let val = |p: &TrajectoryPoint| match mode {
4825 BcFitMode::Drop => drop_offset - p.position.y,
4826 BcFitMode::Velocity => p.velocity_magnitude,
4827 };
4828 for i in 0..points.len() {
4829 if points[i].position.x >= target_dist {
4830 if i == 0 {
4831 return Some(val(&points[0]));
4832 }
4833 let p1 = &points[i - 1];
4834 let p2 = &points[i];
4835 let dx = p2.position.x - p1.position.x;
4836 if dx.abs() < 1e-9 {
4837 return Some(val(p2));
4838 }
4839 let t = (target_dist - p1.position.x) / dx;
4840 return Some(val(p1) + t * (val(p2) - val(p1)));
4841 }
4842 }
4843 None
4844}
4845
4846fn fit_residual_sse(
4847 trajectory: &[TrajectoryPoint],
4848 observations: &[(f64, f64)],
4849 mode: BcFitMode,
4850 drop_offset: f64,
4851) -> Option<f64> {
4852 if observations.is_empty() {
4853 return None;
4854 }
4855 let mut total = 0.0;
4856 for (target_dist, target_val) in observations {
4857 let value = fit_value_at(trajectory, *target_dist, mode, drop_offset)?;
4860 let error = value - target_val;
4861 total += error * error;
4862 }
4863 Some(total)
4864}
4865
4866#[allow(clippy::too_many_arguments)] pub fn estimate_bc_fit(
4884 velocity: f64,
4885 mass: f64,
4886 diameter: f64,
4887 points: &[(f64, f64)],
4888 drag_model: DragModel,
4889 mode: BcFitMode,
4890 atmosphere: AtmosphericConditions,
4891 zero_range: Option<f64>,
4892 sight_height: f64,
4893) -> Result<BcEstimate, BallisticsError> {
4894 if points.is_empty() {
4895 return Err(BallisticsError::from(
4896 "No data points provided for BC estimation.".to_string(),
4897 ));
4898 }
4899 let max_dist = points.iter().map(|(d, _)| *d).fold(0.0_f64, f64::max);
4900 let drop_offset = if zero_range.is_some() {
4903 sight_height
4904 } else {
4905 0.0
4906 };
4907
4908 let sse = |bc_value: f64| -> Option<f64> {
4910 let mut inputs = BallisticInputs {
4911 muzzle_velocity: velocity,
4912 bc_value,
4913 bc_type: drag_model,
4914 bullet_mass: mass,
4915 bullet_diameter: diameter,
4916 sight_height,
4917 ..Default::default()
4918 };
4919 if let Some(zr) = zero_range {
4922 let za = calculate_zero_angle_with_conditions(
4928 inputs.clone(),
4929 zr,
4930 sight_height,
4931 WindConditions::default(),
4932 atmosphere.clone(),
4933 )
4934 .ok()?;
4935 inputs.muzzle_angle = za;
4936 }
4937 let mut solver =
4938 TrajectorySolver::new(inputs, WindConditions::default(), atmosphere.clone());
4939 solver.set_max_range(max_dist * 1.5);
4940 let result = solver.solve().ok()?;
4941 fit_residual_sse(&result.points, points, mode, drop_offset)
4942 };
4943
4944 let (bc_min, bc_max) = match drag_model {
4948 DragModel::G7 => (0.05, 0.70),
4949 _ => (0.10, 1.20),
4950 };
4951
4952 let mut best_bc = f64::NAN;
4954 let mut best_sse = f64::MAX;
4955 let mut bc = bc_min;
4956 while bc <= bc_max + 1e-9 {
4957 if let Some(s) = sse(bc) {
4958 if s < best_sse {
4959 best_sse = s;
4960 best_bc = bc;
4961 }
4962 }
4963 bc += 0.01;
4964 }
4965 if !best_bc.is_finite() {
4966 return Err(BallisticsError::from(
4967 "Unable to estimate BC from provided data. Check that the values and units are correct."
4968 .to_string(),
4969 ));
4970 }
4971
4972 let lo = (best_bc - 0.01).max(bc_min);
4974 let hi = (best_bc + 0.01).min(bc_max);
4975 let mut bc = lo;
4976 while bc <= hi + 1e-9 {
4977 if let Some(s) = sse(bc) {
4978 if s < best_sse {
4979 best_sse = s;
4980 best_bc = bc;
4981 }
4982 }
4983 bc += 0.001;
4984 }
4985
4986 let at_bound = best_bc <= bc_min + 0.011 || best_bc >= bc_max - 0.011;
4989 let rms_error = (best_sse / points.len() as f64).sqrt();
4992 Ok(BcEstimate {
4993 bc: best_bc,
4994 rms_error,
4995 drag_model,
4996 mode,
4997 at_bound,
4998 })
4999}
5000
5001pub fn estimate_bc_from_trajectory(
5004 velocity: f64,
5005 mass: f64,
5006 diameter: f64,
5007 points: &[(f64, f64)], ) -> Result<f64, BallisticsError> {
5009 estimate_bc_fit(
5010 velocity,
5011 mass,
5012 diameter,
5013 points,
5014 DragModel::G1,
5015 BcFitMode::Drop,
5016 AtmosphericConditions::default(),
5017 None,
5018 0.05,
5019 )
5020 .map(|e| e.bc)
5021}
5022
5023use rand;
5025use rand_distr;
5026
5027#[cfg(test)]
5028mod mba737_powder_resolution_tests {
5029 use super::*;
5030
5031 #[test]
5032 fn linear_model_cold_powder_subtracts() {
5033 let v = resolve_powder_adjusted_velocity(823.0, 11.1, true, 0.5486, 21.1, None, None);
5035 assert!((v - (823.0 + 0.5486 * (11.1 - 21.1))).abs() < 1e-12);
5036 assert!(v < 823.0);
5037 }
5038
5039 #[test]
5040 fn linear_model_hot_powder_adds() {
5041 let v = resolve_powder_adjusted_velocity(823.0, 31.1, true, 0.5486, 21.1, None, None);
5042 assert!((v - (823.0 + 0.5486 * 10.0)).abs() < 1e-12);
5043 }
5044
5045 #[test]
5046 fn disabled_flag_is_passthrough() {
5047 let v = resolve_powder_adjusted_velocity(823.0, 40.0, false, 0.5486, 21.1, None, None);
5048 assert_eq!(v, 823.0);
5049 }
5050
5051 #[test]
5052 fn curve_overrides_linear_and_interpolates_at_powder_temp() {
5053 let curve = [(4.4, 798.6), (21.1, 823.0), (37.8, 841.2)];
5054 let v = resolve_powder_adjusted_velocity(
5056 823.0,
5057 30.0,
5058 true,
5059 99.0,
5060 21.1,
5061 Some(&curve),
5062 Some(4.4),
5063 );
5064 assert!((v - 798.6).abs() < 1e-9);
5065 }
5066
5067 #[test]
5068 fn curve_falls_back_to_ambient_and_clamps() {
5069 let curve = [(4.4, 798.6), (37.8, 841.2)];
5070 let v = resolve_powder_adjusted_velocity(823.0, -40.0, true, 1.0, 21.1, Some(&curve), None);
5072 assert!((v - 798.6).abs() < 1e-9);
5073 let v_hot =
5074 resolve_powder_adjusted_velocity(823.0, 60.0, true, 1.0, 21.1, Some(&curve), None);
5075 assert!((v_hot - 841.2).abs() < 1e-9);
5076 }
5077
5078 #[test]
5079 fn empty_curve_suppresses_linear_fallback() {
5080 let v = resolve_powder_adjusted_velocity(823.0, 40.0, true, 0.5486, 21.1, Some(&[]), None);
5082 assert_eq!(v, 823.0);
5083 }
5084
5085 #[test]
5086 fn sweep_huge_range_errors_instead_of_overflowing() {
5087 assert!(parse_powder_sweep("0:1e20:1").is_err());
5090 assert!(parse_powder_sweep("0:1e308:1e-3").is_err());
5091 }
5092
5093 #[test]
5094 fn sweep_fractional_step_keeps_end_row() {
5095 let rows = parse_powder_sweep("0:0.3:0.1").unwrap();
5097 assert_eq!(rows.len(), 4);
5098 assert!((rows[3] - 0.3).abs() < 1e-9);
5099 }
5100
5101 #[test]
5102 fn solver_and_helper_agree_on_linear_model() {
5103 let inputs = BallisticInputs {
5105 use_powder_sensitivity: true,
5106 powder_temp_sensitivity: 0.5486,
5107 powder_temp: 21.1,
5108 temperature: 4.4,
5109 ..Default::default()
5110 };
5111 let expected = resolve_powder_adjusted_velocity(
5112 inputs.muzzle_velocity,
5113 inputs.temperature,
5114 true,
5115 0.5486,
5116 21.1,
5117 None,
5118 None,
5119 );
5120 let solver = TrajectorySolver::new(
5121 inputs,
5122 WindConditions::default(),
5123 AtmosphericConditions::default(),
5124 );
5125 assert!((solver.inputs.muzzle_velocity - expected).abs() < 1e-12);
5126 }
5127}
5128
5129#[cfg(test)]
5130mod mba1302_solver_seam_tests {
5131 use super::*;
5132 use crate::wind::WindSegment;
5133
5134 #[test]
5135 fn authoritative_station_atmosphere_preserves_explicit_standard_values_at_altitude() {
5136 let atmosphere = AtmosphericConditions {
5137 temperature: 15.0,
5138 pressure: 1013.25,
5139 humidity: 50.0,
5140 altitude: 2_000.0,
5141 };
5142 let legacy = TrajectorySolver::new(
5143 BallisticInputs::default(),
5144 WindConditions::default(),
5145 atmosphere.clone(),
5146 );
5147 let authoritative = TrajectorySolver::new_with_resolved_station_atmosphere(
5148 BallisticInputs::default(),
5149 WindConditions::default(),
5150 atmosphere,
5151 );
5152
5153 let (legacy_density, _, legacy_temp_c, legacy_pressure_hpa) = legacy.resolved_atmosphere();
5154 let (authoritative_density, _, authoritative_temp_c, authoritative_pressure_hpa) =
5155 authoritative.resolved_atmosphere();
5156 let (icao_temp_k, icao_pressure_pa) =
5157 crate::atmosphere::calculate_icao_standard_atmosphere(2_000.0);
5158 let (expected_authoritative_density, _) =
5159 crate::atmosphere::calculate_atmosphere(2_000.0, Some(15.0), Some(1013.25), 50.0);
5160
5161 assert!((legacy_temp_c - (icao_temp_k - 273.15)).abs() < 1e-12);
5162 assert!((legacy_pressure_hpa - icao_pressure_pa / 100.0).abs() < 1e-12);
5163 assert_eq!(authoritative_temp_c.to_bits(), 15.0_f64.to_bits());
5164 assert_eq!(authoritative_pressure_hpa.to_bits(), 1013.25_f64.to_bits());
5165 assert_eq!(
5166 authoritative_density.to_bits(),
5167 expected_authoritative_density.to_bits()
5168 );
5169 assert!(
5170 (authoritative_density - legacy_density).abs() > 0.1,
5171 "explicit standard values at altitude must differ from ICAO-at-altitude: explicit={authoritative_density}, ICAO={legacy_density}"
5172 );
5173 }
5174
5175 #[test]
5184 fn precomputed_absolute_resolution_via_authoritative_matches_legacy_new() {
5185 for (temperature, pressure, altitude) in [
5186 (15.0, 1013.25, 0.0), (15.0, 1013.25, 2000.0), (-5.0, 850.0, 2000.0), (22.0, 950.0, 500.0),
5190 ] {
5191 let atmosphere = AtmosphericConditions {
5192 temperature,
5193 pressure,
5194 humidity: 50.0,
5195 altitude,
5196 };
5197 let legacy = TrajectorySolver::new(
5198 BallisticInputs::default(),
5199 WindConditions::default(),
5200 atmosphere.clone(),
5201 );
5202
5203 let (resolved_temp_c, resolved_pressure_hpa) =
5204 crate::atmosphere::resolve_station_conditions_with_pressure_mode(
5205 temperature,
5206 pressure,
5207 altitude,
5208 crate::atmosphere::PressureReferenceMode::Absolute,
5209 );
5210 let precomputed_atmosphere = AtmosphericConditions {
5211 temperature: resolved_temp_c,
5212 pressure: resolved_pressure_hpa,
5213 humidity: 50.0,
5214 altitude,
5215 };
5216 let precomputed = TrajectorySolver::new_with_resolved_station_atmosphere(
5217 BallisticInputs::default(),
5218 WindConditions::default(),
5219 precomputed_atmosphere,
5220 );
5221
5222 let (legacy_density, legacy_sos, legacy_temp_c, legacy_pressure_hpa) =
5223 legacy.resolved_atmosphere();
5224 let (pre_density, pre_sos, pre_temp_c, pre_pressure_hpa) =
5225 precomputed.resolved_atmosphere();
5226
5227 assert_eq!(
5228 legacy_temp_c.to_bits(),
5229 pre_temp_c.to_bits(),
5230 "temperature=({temperature}, {pressure}, {altitude})"
5231 );
5232 assert_eq!(
5233 legacy_pressure_hpa.to_bits(),
5234 pre_pressure_hpa.to_bits(),
5235 "pressure=({temperature}, {pressure}, {altitude})"
5236 );
5237 assert_eq!(legacy_density.to_bits(), pre_density.to_bits());
5238 assert_eq!(legacy_sos.to_bits(), pre_sos.to_bits());
5239 }
5240 }
5241
5242 fn configured_euler_zero(vertical_wind_mps: f64, time_step_s: f64) -> TrajectorySolver {
5243 let inputs = BallisticInputs {
5244 muzzle_velocity: 800.0,
5245 bc_value: 0.5,
5246 bc_type: DragModel::G7,
5247 bullet_mass: 0.0109,
5248 bullet_diameter: 0.00782,
5249 bullet_length: 0.0309,
5250 sight_height: 0.05,
5251 ground_threshold: -100.0,
5252 use_rk4: false,
5253 use_adaptive_rk45: false,
5254 ..BallisticInputs::default()
5255 };
5256 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(
5257 inputs,
5258 WindConditions::default(),
5259 AtmosphericConditions::default(),
5260 );
5261 solver.set_max_range(300.0);
5262 solver.set_time_step(time_step_s);
5263 if vertical_wind_mps != 0.0 {
5264 solver.set_wind_segments(vec![WindSegment {
5265 speed_kmh: 0.0,
5266 angle_deg: 0.0,
5267 until_m: 400.0,
5268 vertical_mps: vertical_wind_mps,
5269 }]);
5270 }
5271 solver
5272 }
5273
5274 #[test]
5275 fn inclined_shot_zeroes_like_a_level_rifle() {
5276 const ZERO_DISTANCE_M: f64 = 91.44; const SIGHT_HEIGHT_M: f64 = 0.0381; let inputs = BallisticInputs {
5285 bc_value: 0.5,
5286 bullet_mass: 150.0 * 0.06479891 / 1000.0,
5287 muzzle_velocity: 2700.0 * 0.3048,
5288 sight_height: SIGHT_HEIGHT_M,
5289 ..Default::default()
5290 };
5291
5292 let mut level = inputs.clone();
5293 level.shooting_angle = 0.0;
5294 let level_angle = TrajectorySolver::new(level, Default::default(), Default::default())
5295 .find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
5296 .expect("level zero must solve");
5297
5298 let mut inclined = inputs;
5299 inclined.shooting_angle = 5.71_f64.to_radians();
5300 let inclined_angle =
5301 TrajectorySolver::new(inclined, Default::default(), Default::default())
5302 .find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
5303 .expect("MBA-1412: a 5.71 deg incline at a 100 yd zero must be solvable");
5304
5305 assert!(
5306 (inclined_angle - level_angle).abs() < 1e-9,
5307 "zeroing is level-rifle sight geometry; incline must not move the solved zero: \
5308 level={level_angle}, inclined={inclined_angle}"
5309 );
5310 }
5311
5312 #[test]
5313 fn configured_zero_keeps_segments_method_and_time_step_then_sets_base_angle() {
5314 const TARGET_DISTANCE_M: f64 = 150.0;
5315 const TARGET_HEIGHT_M: f64 = 0.05;
5316
5317 let mut segmented = configured_euler_zero(-10.0, 0.02);
5320 let coarse_height = segmented
5321 .zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
5322 .expect("coarse configured trial")
5323 .expect("coarse trial reaches target");
5324 let mut fine = segmented.clone();
5325 fine.set_time_step(0.001);
5326 let fine_height = fine
5327 .zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
5328 .expect("fine configured trial")
5329 .expect("fine trial reaches target");
5330 assert!(
5331 (coarse_height - fine_height).abs() > 1e-5,
5332 "configured Euler step must affect zero trials: coarse={coarse_height}, fine={fine_height}"
5333 );
5334
5335 let segmented_angle = segmented
5336 .calculate_and_set_zero_angle(
5337 TARGET_DISTANCE_M,
5338 TARGET_HEIGHT_M,
5339 ZeroTargetFrame::SightLine,
5340 )
5341 .expect("segmented zero");
5342 assert_eq!(
5343 segmented.inputs.muzzle_angle.to_bits(),
5344 segmented_angle.to_bits(),
5345 "successful zero must install its angle on the configured solver"
5346 );
5347 assert_eq!(segmented.time_step.to_bits(), 0.02_f64.to_bits());
5348 assert_eq!(segmented.max_range.to_bits(), 300.0_f64.to_bits());
5349 assert!(segmented.wind_sock.is_some());
5350 assert_eq!(
5351 segmented.station_atmosphere_resolution,
5352 StationAtmosphereResolution::Authoritative
5353 );
5354 let zero_height = segmented
5355 .zero_trial_height_at(
5356 segmented_angle,
5357 TARGET_DISTANCE_M,
5358 ZeroTargetFrame::SightLine,
5359 )
5360 .expect("verify segmented zero")
5361 .expect("zeroed trial reaches target");
5362 assert!(
5363 (zero_height - TARGET_HEIGHT_M).abs() < 0.0001,
5364 "configured zero missed target: height={zero_height}"
5365 );
5366
5367 let mut calm = configured_euler_zero(0.0, 0.02);
5368 let calm_angle = calm
5369 .calculate_and_set_zero_angle(
5370 TARGET_DISTANCE_M,
5371 TARGET_HEIGHT_M,
5372 ZeroTargetFrame::SightLine,
5373 )
5374 .expect("calm zero");
5375 assert!(
5376 (segmented_angle - calm_angle).abs() > 1e-5,
5377 "segmented vertical wind must participate in zero trials: segmented={segmented_angle}, calm={calm_angle}"
5378 );
5379 }
5380}
5381
5382#[cfg(test)]
5383mod result_sanity_tests {
5384 use super::*;
5385
5386 fn default_solver() -> TrajectorySolver {
5387 TrajectorySolver::new(
5388 BallisticInputs::default(),
5389 WindConditions::default(),
5390 AtmosphericConditions::default(),
5391 )
5392 }
5393
5394 fn minimal_result() -> TrajectoryResult {
5395 TrajectoryResult {
5396 max_range: 100.0,
5397 max_height: 1.0,
5398 time_of_flight: 0.5,
5399 impact_velocity: 700.0,
5400 impact_energy: 2450.0,
5401 projectile_mass_kg: 0.01,
5402 line_of_sight_height_m: 1.5,
5403 station_speed_of_sound_mps: 340.0,
5404 termination: TrajectoryTermination::MaxRange,
5405 points: vec![],
5406 sampled_points: None,
5407 min_pitch_damping: None,
5408 transonic_mach: None,
5409 angular_state: None,
5410 max_yaw_angle: None,
5411 max_precession_angle: None,
5412 aerodynamic_jump: None,
5413 mach_1_2_distance_m: None,
5414 mach_1_0_distance_m: None,
5415 mach_0_9_distance_m: None,
5416 }
5417 }
5418
5419 #[test]
5420 fn mba1293_negative_scalars_fail_the_result_postcondition() {
5421 let solver = default_solver();
5422 solver
5423 .validate_result_sanity(&minimal_result())
5424 .expect("a sane result must pass");
5425
5426 for (name, mutate) in [
5427 (
5428 "max_range",
5429 (|r| r.max_range = -50.588) as fn(&mut TrajectoryResult),
5430 ),
5431 ("time_of_flight", |r| r.time_of_flight = -1.0),
5432 ("impact_velocity", |r| r.impact_velocity = -700.0),
5433 ("impact_energy", |r| r.impact_energy = -1.0),
5434 ] {
5435 let mut result = minimal_result();
5436 mutate(&mut result);
5437 let error = solver
5438 .validate_result_sanity(&result)
5439 .expect_err("negative scalar must fail");
5440 assert!(
5441 error.to_string().contains(name),
5442 "error for {name} did not name the field: {error}"
5443 );
5444 }
5445 }
5446
5447 #[test]
5448 fn mba1293_speed_budget_bounds_legitimate_states_and_rejects_divergence() {
5449 let solver = default_solver();
5450 let mv = solver.inputs.muzzle_velocity;
5451
5452 let position = Vector3::new(10.0, 0.0, 0.0);
5454 solver
5455 .validate_integration_state(&position, &Vector3::new(mv, 0.0, 0.0), 0.01)
5456 .expect("muzzle-speed state must pass");
5457
5458 let error = solver
5460 .validate_integration_state(&position, &Vector3::new(-13.0 * mv, 0.0, 0.0), 0.01)
5461 .expect_err("13x muzzle speed must fail the budget");
5462 assert!(error.to_string().contains("diverged"), "{error}");
5463
5464 let after_fall = mv + crate::constants::G_ACCEL_MPS2 * 60.0;
5466 solver
5467 .validate_integration_state(&position, &Vector3::new(0.0, -after_fall, 0.0), 60.0)
5468 .expect("gravity-accelerated speed within g*t must pass");
5469 }
5470}
5471
5472#[cfg(test)]
5473mod trajectory_point_budget_tests {
5474 use super::*;
5475 use crate::MAX_TRAJECTORY_SAMPLES;
5476
5477 fn solver_with_budget(
5478 use_rk4: bool,
5479 use_adaptive_rk45: bool,
5480 point_budget: usize,
5481 max_range: f64,
5482 ) -> TrajectorySolver {
5483 let inputs = BallisticInputs {
5484 use_rk4,
5485 use_adaptive_rk45,
5486 ground_threshold: f64::NEG_INFINITY,
5487 ..BallisticInputs::default()
5488 };
5489 let mut solver = TrajectorySolver::new(
5490 inputs,
5491 WindConditions::default(),
5492 AtmosphericConditions::default(),
5493 );
5494 solver.max_trajectory_points = point_budget;
5495 solver.set_max_range(max_range);
5496 solver.set_time_step(0.001);
5497 solver
5498 }
5499
5500 #[test]
5501 fn mba1283_every_solver_errors_instead_of_exceeding_point_budget() {
5502 for (mode, use_rk4, use_adaptive_rk45) in [
5503 ("Euler", false, false),
5504 ("RK4", true, false),
5505 ("RK45", true, true),
5506 ] {
5507 let error = solver_with_budget(use_rk4, use_adaptive_rk45, 3, 10.0)
5508 .solve()
5509 .expect_err("a solve requiring more than three points must fail");
5510 assert!(
5511 error.to_string().contains("point limit of 3"),
5512 "unexpected {mode} point-budget error: {error}"
5513 );
5514 }
5515 }
5516
5517 #[test]
5518 fn mba1283_interpolated_endpoint_counts_toward_point_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 result = solver_with_budget(use_rk4, use_adaptive_rk45, 2, 0.1)
5525 .solve()
5526 .expect("the initial point plus exact endpoint fit a two-point budget");
5527 assert_eq!(result.points.len(), 2, "unexpected {mode} point count");
5528
5529 let error = solver_with_budget(use_rk4, use_adaptive_rk45, 1, 0.1)
5530 .solve()
5531 .expect_err("the exact endpoint must not exceed a one-point budget");
5532 assert!(
5533 error.to_string().contains("point limit of 1"),
5534 "unexpected {mode} endpoint-budget error: {error}"
5535 );
5536 }
5537 }
5538
5539 #[test]
5540 fn mba1299_every_solver_preflights_the_sample_budget() {
5541 for (mode, use_rk4, use_adaptive_rk45) in [
5542 ("Euler", false, false),
5543 ("RK4", true, false),
5544 ("RK45", true, true),
5545 ] {
5546 let inputs = BallisticInputs {
5547 use_rk4,
5548 use_adaptive_rk45,
5549 enable_trajectory_sampling: true,
5550 sample_interval: 1.0,
5551 ground_threshold: f64::NEG_INFINITY,
5552 ..BallisticInputs::default()
5553 };
5554 let mut solver = TrajectorySolver::new(
5555 inputs,
5556 WindConditions::default(),
5557 AtmosphericConditions::default(),
5558 );
5559 solver.set_max_range(MAX_TRAJECTORY_SAMPLES as f64);
5560 solver.max_trajectory_points = 0;
5563
5564 let error = solver
5565 .solve()
5566 .expect_err("an over-limit sample grid must fail before integration");
5567 assert!(
5568 error
5569 .to_string()
5570 .contains("trajectory sample limit of 250000 exceeded"),
5571 "unexpected {mode} sample-budget error: {error}"
5572 );
5573 }
5574 }
5575
5576 #[test]
5577 fn mba1299_normal_sampling_does_not_change_solver_results() {
5578 for (mode, use_rk4, use_adaptive_rk45) in [
5579 ("Euler", false, false),
5580 ("RK4", true, false),
5581 ("RK45", true, true),
5582 ] {
5583 let solve = |enable_trajectory_sampling| {
5584 let inputs = BallisticInputs {
5585 use_rk4,
5586 use_adaptive_rk45,
5587 enable_trajectory_sampling,
5588 sample_interval: 0.5,
5589 ground_threshold: f64::NEG_INFINITY,
5590 ..BallisticInputs::default()
5591 };
5592 let mut solver = TrajectorySolver::new(
5593 inputs,
5594 WindConditions::default(),
5595 AtmosphericConditions::default(),
5596 );
5597 solver.set_max_range(2.0);
5598 solver.solve().expect("normal short-range solve")
5599 };
5600
5601 let baseline = solve(false);
5602 let sampled = solve(true);
5603 for (field, left, right) in [
5604 ("max_range", baseline.max_range, sampled.max_range),
5605 ("max_height", baseline.max_height, sampled.max_height),
5606 (
5607 "time_of_flight",
5608 baseline.time_of_flight,
5609 sampled.time_of_flight,
5610 ),
5611 (
5612 "impact_velocity",
5613 baseline.impact_velocity,
5614 sampled.impact_velocity,
5615 ),
5616 (
5617 "impact_energy",
5618 baseline.impact_energy,
5619 sampled.impact_energy,
5620 ),
5621 ] {
5622 assert_eq!(
5623 left.to_bits(),
5624 right.to_bits(),
5625 "{mode} sampling changed {field}"
5626 );
5627 }
5628 assert_eq!(baseline.points.len(), sampled.points.len());
5629 for (index, (left, right)) in baseline.points.iter().zip(&sampled.points).enumerate() {
5630 assert_eq!(
5631 left.time.to_bits(),
5632 right.time.to_bits(),
5633 "{mode} point {index}"
5634 );
5635 assert_eq!(
5636 left.position.map(f64::to_bits),
5637 right.position.map(f64::to_bits),
5638 "{mode} point {index} position"
5639 );
5640 assert_eq!(
5641 left.velocity_magnitude.to_bits(),
5642 right.velocity_magnitude.to_bits(),
5643 "{mode} point {index} velocity"
5644 );
5645 assert_eq!(
5646 left.kinetic_energy.to_bits(),
5647 right.kinetic_energy.to_bits(),
5648 "{mode} point {index} energy"
5649 );
5650 }
5651 assert!(baseline.sampled_points.is_none());
5652 let samples = sampled
5653 .sampled_points
5654 .expect("sampling-enabled solve should return observations");
5655 assert_eq!(
5656 samples
5657 .iter()
5658 .map(|sample| sample.distance_m)
5659 .collect::<Vec<_>>(),
5660 vec![0.0, 0.5, 1.0, 1.5, 2.0],
5661 "{mode} normal sampling grid changed"
5662 );
5663 }
5664 }
5665}
5666
5667#[cfg(test)]
5668mod monte_carlo_result_tests {
5669 use super::*;
5670
5671 fn make_results(impact_positions: Vec<Vector3<f64>>) -> MonteCarloResults {
5672 let count = impact_positions.len();
5673 MonteCarloResults {
5674 ranges: vec![500.0; count],
5675 impact_velocities: vec![300.0; count],
5676 impact_positions,
5677 }
5678 }
5679
5680 #[test]
5681 fn target_plane_cep_excludes_shortfall_markers() {
5682 let mut positions: Vec<Vector3<f64>> = (1..=5)
5683 .map(|radius| Vector3::new(0.0, radius as f64, 0.0))
5684 .collect();
5685 positions.extend((0..5).map(|_| Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)));
5686 let results = make_results(positions);
5687
5688 assert_eq!(results.target_arrival_count(), 5);
5689 assert_eq!(results.target_shortfall_fraction(), 0.5);
5690 assert_eq!(results.target_plane_cep(), Some(3.0));
5691
5692 let one_shortfall = make_results(vec![
5693 Vector3::new(0.0, 1.0, 0.0),
5694 Vector3::new(0.0, 2.0, 0.0),
5695 Vector3::new(0.0, 3.0, 0.0),
5696 Vector3::new(0.0, 4.0, 0.0),
5697 Vector3::new(0.0, 5.0, 0.0),
5698 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5699 ]);
5700 assert_eq!(one_shortfall.target_plane_cep(), Some(3.0));
5701 }
5702
5703 #[test]
5704 fn all_shortfalls_have_no_cep_but_still_count_as_misses() {
5705 let all_shortfalls = make_results(vec![
5706 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5707 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5708 ]);
5709 assert_eq!(all_shortfalls.target_arrival_count(), 0);
5710 assert_eq!(all_shortfalls.target_shortfall_fraction(), 1.0);
5711 assert_eq!(all_shortfalls.target_plane_cep(), None);
5712 assert_eq!(all_shortfalls.hit_probability(0.3), 0.0);
5713
5714 let one_hit_one_shortfall = make_results(vec![
5715 Vector3::new(0.0, 0.1, 0.0),
5716 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5717 ]);
5718 assert_eq!(one_hit_one_shortfall.hit_probability(0.3), 0.5);
5719 }
5720
5721 #[test]
5723 fn rect_hit_probability_checks_independent_axis_halves() {
5724 let results = make_results(vec![
5725 Vector3::new(0.0, 0.1, 0.1),
5727 Vector3::new(0.0, 0.0, 0.2),
5729 Vector3::new(0.0, 0.0, 0.201),
5731 Vector3::new(0.0, 0.301, 0.0),
5733 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5735 ]);
5736 assert!((results.rect_hit_probability(0.4, 0.6) - 0.4).abs() < 1e-12);
5738 }
5739
5740 #[test]
5741 fn rect_hit_probability_matches_circular_hit_probability_for_a_centered_hit() {
5742 let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
5743 assert_eq!(results.rect_hit_probability(0.5, 0.5), 1.0);
5744 assert_eq!(results.hit_probability(0.3), 1.0);
5745 }
5746
5747 #[test]
5748 fn rect_hit_probability_is_zero_for_empty_or_nonpositive_dimensions() {
5749 let empty = make_results(vec![]);
5750 assert_eq!(empty.rect_hit_probability(1.0, 1.0), 0.0);
5751
5752 let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
5753 assert_eq!(results.rect_hit_probability(0.0, 1.0), 0.0);
5754 assert_eq!(results.rect_hit_probability(1.0, 0.0), 0.0);
5755 assert_eq!(results.rect_hit_probability(-1.0, 1.0), 0.0);
5756 }
5757}
5758
5759#[cfg(test)]
5760mod monte_carlo_seeded_tests {
5761 use super::*;
5762
5763 fn seeded_test_fixture() -> (BallisticInputs, WindConditions) {
5772 (
5773 BallisticInputs {
5774 muzzle_velocity: 800.0,
5775 ..BallisticInputs::default()
5776 },
5777 WindConditions::default(),
5778 )
5779 }
5780
5781 fn loose_params() -> MonteCarloParams {
5791 MonteCarloParams {
5792 num_simulations: 1, velocity_std_dev: 3.0,
5794 angle_std_dev: 3.5e-4,
5795 bc_std_dev: 0.01,
5796 wind_speed_std_dev: 1.0,
5797 target_distance: Some(300.0),
5798 base_wind_speed: 0.0,
5799 base_wind_direction: 0.0,
5800 azimuth_std_dev: 3.5e-4,
5801 }
5802 }
5803
5804 fn mixed_arrival_params() -> MonteCarloParams {
5812 MonteCarloParams {
5813 num_simulations: 1,
5814 target_distance: Some(1920.0),
5815 ..MonteCarloParams::default()
5816 }
5817 }
5818
5819 #[test]
5835 fn legacy_seeded_estimates_are_pinned_bit_for_bit() {
5836 let (inputs, wind) = seeded_test_fixture();
5837 let params = MonteCarloParams {
5838 num_simulations: 200,
5839 target_distance: Some(500.0),
5840 ..MonteCarloParams::default()
5841 };
5842
5843 let results = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5844 inputs,
5845 wind,
5846 params,
5847 0.01,
5848 0x1352_5EED,
5849 )
5850 .expect("seeded legacy run");
5851
5852 assert_eq!(results.ranges.len(), 200, "ranges length");
5857 assert_eq!(
5858 results.impact_velocities.len(),
5859 200,
5860 "impact_velocities length"
5861 );
5862 assert_eq!(
5863 results.impact_positions.len(),
5864 200,
5865 "impact_positions length"
5866 );
5867
5868 assert_eq!(
5870 results.hit_probability(DEFAULT_HIT_RADIUS_M).to_bits(),
5871 0.14_f64.to_bits(),
5872 "hit_probability = {:?}",
5873 results.hit_probability(DEFAULT_HIT_RADIUS_M)
5874 );
5875
5876 let expected_ranges: [f64; 3] = [1907.972891143359, 1936.408435469319, 1912.8150447617645];
5877 let expected_velocities: [f64; 3] =
5878 [238.6187151542299, 239.91923651600106, 243.14112455427164];
5879 let expected_positions: [(f64, f64, f64); 3] = [
5880 (0.0, -0.0643556039548101, 0.7344970252014579),
5881 (0.0, 0.5769422971539162, 0.27227201756386726),
5882 (0.0, -0.7440425792472842, 0.1541804446282822),
5883 ];
5884
5885 for (i, expected) in expected_ranges.iter().enumerate() {
5886 assert_eq!(
5887 results.ranges[i].to_bits(),
5888 expected.to_bits(),
5889 "ranges[{i}] = {:?}, pinned {expected:?}",
5890 results.ranges[i]
5891 );
5892 }
5893 for (i, expected) in expected_velocities.iter().enumerate() {
5894 assert_eq!(
5895 results.impact_velocities[i].to_bits(),
5896 expected.to_bits(),
5897 "impact_velocities[{i}] = {:?}, pinned {expected:?}",
5898 results.impact_velocities[i]
5899 );
5900 }
5901 for (i, (x, y, z)) in expected_positions.iter().enumerate() {
5902 let actual = results.impact_positions[i];
5903 assert_eq!(
5904 actual.x.to_bits(),
5905 x.to_bits(),
5906 "impact_positions[{i}].x = {:?}",
5907 actual.x
5908 );
5909 assert_eq!(
5910 actual.y.to_bits(),
5911 y.to_bits(),
5912 "impact_positions[{i}].y = {:?}",
5913 actual.y
5914 );
5915 assert_eq!(
5916 actual.z.to_bits(),
5917 z.to_bits(),
5918 "impact_positions[{i}].z = {:?}",
5919 actual.z
5920 );
5921 }
5922 }
5923
5924 #[test]
5925 fn seeded_runs_are_deterministic_and_match_the_using_rng_path() {
5926 let inputs = BallisticInputs {
5927 muzzle_velocity: 800.0,
5928 ..BallisticInputs::default()
5929 };
5930 let params = MonteCarloParams {
5931 num_simulations: 64,
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.01,
5941 42,
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.01,
5949 42,
5950 )
5951 .expect("seeded run b");
5952
5953 assert_eq!(a.ranges.len(), b.ranges.len());
5954 for (ra, rb) in a.ranges.iter().zip(b.ranges.iter()) {
5955 assert_eq!(ra.to_bits(), rb.to_bits());
5956 }
5957 for (pa, pb) in a.impact_positions.iter().zip(b.impact_positions.iter()) {
5958 assert_eq!(pa.x.to_bits(), pb.x.to_bits());
5959 assert_eq!(pa.y.to_bits(), pb.y.to_bits());
5960 assert_eq!(pa.z.to_bits(), pb.z.to_bits());
5961 }
5962 }
5963
5964 #[test]
5965 fn different_seeds_generally_produce_different_draws() {
5966 let inputs = BallisticInputs {
5967 muzzle_velocity: 800.0,
5968 ..BallisticInputs::default()
5969 };
5970 let params = MonteCarloParams {
5971 num_simulations: 32,
5972 velocity_std_dev: 5.0,
5973 target_distance: Some(200.0),
5974 ..MonteCarloParams::default()
5975 };
5976
5977 let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5978 inputs.clone(),
5979 WindConditions::default(),
5980 params.clone(),
5981 0.0,
5982 1,
5983 )
5984 .expect("seeded run a");
5985 let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5986 inputs,
5987 WindConditions::default(),
5988 params,
5989 0.0,
5990 2,
5991 )
5992 .expect("seeded run b");
5993
5994 assert_ne!(a.impact_velocities, b.impact_velocities);
5995 }
5996
5997 #[test]
5998 fn adaptive_stops_at_target_half_width_on_an_easy_case() {
5999 let (inputs, wind) = seeded_test_fixture();
6000 let conv = McConvergence {
6001 target_half_width: 0.05,
6002 ..Default::default()
6003 };
6004 let r = run_monte_carlo_adaptive_seeded(
6005 &inputs,
6006 &wind,
6007 &loose_params(),
6008 &conv,
6009 DEFAULT_HIT_RADIUS_M,
6010 0x1352_ADA9,
6011 )
6012 .unwrap();
6013
6014 assert_eq!(r.stop_reason, McStopReason::TargetHalfWidthMet);
6015 assert!(
6016 (r.ci_high - r.ci_low) / 2.0 <= 0.05 + 1e-12,
6017 "half-width {} exceeds the requested 0.05",
6018 (r.ci_high - r.ci_low) / 2.0
6019 );
6020 assert!(r.samples >= conv.min_samples, "stopped below min_samples");
6021 assert!(r.samples < conv.max_samples, "did not actually stop early");
6022 assert!(r.samples.is_multiple_of(conv.batch_size) || r.samples == conv.min_samples);
6023 assert!(r.ci_low <= r.hit_probability && r.hit_probability <= r.ci_high);
6024
6025 assert_eq!(r.hit_radius_m, DEFAULT_HIT_RADIUS_M);
6027 assert_eq!(r.target_distance_m, 300.0);
6028 assert_eq!(r.confidence_percent, 95);
6029 assert!(
6031 r.samples > 1,
6032 "params.num_simulations must be ignored by the adaptive driver"
6033 );
6034 assert!(
6037 r.mean_impact_velocity_mps > 0.0,
6038 "no impact velocity accumulated"
6039 );
6040 assert!(
6041 r.std_drop_at_target_m > 0.0 && r.std_wind_drift_at_target_m > 0.0,
6042 "dispersion collapsed: drop sd {} drift sd {}",
6043 r.std_drop_at_target_m,
6044 r.std_wind_drift_at_target_m
6045 );
6046 }
6047
6048 #[test]
6049 fn adaptive_caps_at_max_samples_on_an_impossible_target() {
6050 let (inputs, wind) = seeded_test_fixture();
6051 let conv = McConvergence {
6052 target_half_width: 1e-6,
6053 max_samples: 3_000,
6054 batch_size: 500,
6055 min_samples: 1_000,
6056 level: ConfidenceLevel::P95,
6057 };
6058 let r = run_monte_carlo_adaptive_seeded(
6059 &inputs,
6060 &wind,
6061 &loose_params(),
6062 &conv,
6063 DEFAULT_HIT_RADIUS_M,
6064 7,
6065 )
6066 .unwrap();
6067
6068 assert_eq!(r.stop_reason, McStopReason::MaxSamplesReached);
6069 assert_eq!(r.samples, 3_000);
6070 }
6071
6072 #[test]
6081 fn adaptive_stops_between_the_floor_and_the_ceiling() {
6082 let (inputs, wind) = seeded_test_fixture();
6083 let conv = McConvergence {
6084 level: ConfidenceLevel::P95,
6085 target_half_width: 0.03,
6086 min_samples: 0,
6087 max_samples: 10_000,
6088 batch_size: 100,
6089 };
6090 let r = run_monte_carlo_adaptive_seeded(
6091 &inputs,
6092 &wind,
6093 &loose_params(),
6094 &conv,
6095 DEFAULT_HIT_RADIUS_M,
6096 0x1352_5A1D,
6097 )
6098 .unwrap();
6099
6100 assert_eq!(r.stop_reason, McStopReason::TargetHalfWidthMet);
6101 assert!(r.samples > 0);
6102 assert!(
6103 r.samples.is_multiple_of(conv.batch_size),
6104 "samples {} is not a whole number of batches",
6105 r.samples
6106 );
6107 assert!(
6108 r.samples > conv.min_samples,
6109 "stopped on the floor, not on the data"
6110 );
6111 assert!(
6112 r.samples < conv.max_samples,
6113 "ran to the ceiling, so nothing adaptive was exercised"
6114 );
6115 assert!(
6116 r.samples > conv.batch_size,
6117 "stopped on the very first look ({} samples); the multi-batch path is untested",
6118 r.samples
6119 );
6120 assert!((r.ci_high - r.ci_low) / 2.0 <= 0.03 + 1e-12);
6121 assert!(r.ci_low <= r.hit_probability && r.hit_probability <= r.ci_high);
6122 assert_eq!(r.attempts, r.samples, "no trial should have been dropped");
6123 }
6124
6125 #[test]
6133 fn adaptive_runs_a_truncated_final_batch_up_to_max_samples() {
6134 let (inputs, wind) = seeded_test_fixture();
6135 let conv = McConvergence {
6136 level: ConfidenceLevel::P95,
6137 target_half_width: 1e-6,
6138 min_samples: 0,
6139 max_samples: 750,
6140 batch_size: 500,
6141 };
6142 let r = run_monte_carlo_adaptive_seeded(
6143 &inputs,
6144 &wind,
6145 &loose_params(),
6146 &conv,
6147 DEFAULT_HIT_RADIUS_M,
6148 0x1352_7B10,
6149 )
6150 .unwrap();
6151
6152 assert_eq!(r.stop_reason, McStopReason::MaxSamplesReached);
6153 assert_eq!(
6154 r.samples, 750,
6155 "the 250-trial final batch did not run, or was not truncated"
6156 );
6157 assert_eq!(r.attempts, 750);
6158 assert!(!r.samples.is_multiple_of(conv.batch_size));
6159 }
6160
6161 #[test]
6162 fn adaptive_is_deterministic_for_a_seed() {
6163 let (inputs, wind) = seeded_test_fixture();
6164 let conv = McConvergence::default();
6165 let a = run_monte_carlo_adaptive_seeded(
6166 &inputs,
6167 &wind,
6168 &loose_params(),
6169 &conv,
6170 DEFAULT_HIT_RADIUS_M,
6171 99,
6172 )
6173 .unwrap();
6174 let b = run_monte_carlo_adaptive_seeded(
6175 &inputs,
6176 &wind,
6177 &loose_params(),
6178 &conv,
6179 DEFAULT_HIT_RADIUS_M,
6180 99,
6181 )
6182 .unwrap();
6183
6184 assert_eq!(a.hit_probability.to_bits(), b.hit_probability.to_bits());
6185 assert_eq!(a.samples, b.samples);
6186 assert_eq!(a.ci_low.to_bits(), b.ci_low.to_bits());
6187 assert_eq!(a.ci_high.to_bits(), b.ci_high.to_bits());
6188 assert_eq!(
6191 a.mean_impact_velocity_mps.to_bits(),
6192 b.mean_impact_velocity_mps.to_bits()
6193 );
6194 assert_eq!(
6195 a.std_drop_at_target_m.to_bits(),
6196 b.std_drop_at_target_m.to_bits()
6197 );
6198 }
6199
6200 #[test]
6201 fn adaptive_report_carries_schema_method_and_all_four_assumptions() {
6202 let (inputs, wind) = seeded_test_fixture();
6203 let conv = McConvergence {
6207 min_samples: 0,
6208 max_samples: 50,
6209 batch_size: 50,
6210 target_half_width: 1.0,
6211 level: ConfidenceLevel::P90,
6212 };
6213 let r = run_monte_carlo_adaptive_seeded(
6214 &inputs,
6215 &wind,
6216 &mixed_arrival_params(),
6217 &conv,
6218 DEFAULT_HIT_RADIUS_M,
6219 0x1352_D0C5,
6220 )
6221 .unwrap();
6222
6223 assert_eq!(r.schema_version, MC_ADAPTIVE_SCHEMA_VERSION_V1);
6224 assert_eq!(r.schema_version, 1);
6225 assert_eq!(r.method, "anytime_beta_binomial_mixture_cs_v1");
6226 assert_eq!(r.confidence_percent, 90);
6227
6228 assert_eq!(r.attempts, 50, "one full batch was drawn");
6234 assert_eq!(r.samples, 50, "no trial was dropped by the solver");
6235 assert!(
6236 r.arrivals > 0 && r.arrivals < r.samples,
6237 "fixture must split the run: arrivals {} of samples {}",
6238 r.arrivals,
6239 r.samples
6240 );
6241 assert!(
6244 r.arrivals >= 2,
6245 "arrivals {} too few for a sample sd",
6246 r.arrivals
6247 );
6248 assert!(r.std_drop_at_target_m > 0.0 && r.std_impact_velocity_mps > 0.0);
6249 assert!(r.attempts >= r.samples && r.samples >= r.arrivals);
6251
6252 assert_eq!(r.assumptions.len(), 4, "exactly four assumptions expected");
6258 assert_eq!(
6259 r.assumptions[0],
6260 "Sampling uncertainty only: intervals cover Monte Carlo sampling error, not model error in the trajectory solver or its inputs."
6261 );
6262 assert_eq!(
6263 r.assumptions[1],
6264 "Anytime-valid stopping: the beta-binomial mixture confidence sequence keeps its coverage guarantee despite stopping the moment the target half-width is met."
6265 );
6266 assert_eq!(
6267 r.assumptions[2],
6268 "Input dispersions are the independent normal distributions declared in MonteCarloParams; correlations between inputs are not modeled."
6269 );
6270 assert_eq!(
6271 r.assumptions[3],
6272 "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."
6273 );
6274
6275 assert_eq!(
6277 serde_json::to_string(&McStopReason::TargetHalfWidthMet).unwrap(),
6278 "\"target_half_width_met\""
6279 );
6280 assert_eq!(
6281 serde_json::to_string(&McStopReason::MaxSamplesReached).unwrap(),
6282 "\"max_samples_reached\""
6283 );
6284 }
6285
6286 #[test]
6287 fn adaptive_rejects_nonsense_convergence() {
6288 let (inputs, wind) = seeded_test_fixture();
6289 let run = |conv: McConvergence| {
6290 run_monte_carlo_adaptive_seeded(
6291 &inputs,
6292 &wind,
6293 &loose_params(),
6294 &conv,
6295 DEFAULT_HIT_RADIUS_M,
6296 1,
6297 )
6298 .unwrap_err()
6299 };
6300
6301 for bad_width in [0.0, -0.01, f64::NAN] {
6302 let err = run(McConvergence {
6303 target_half_width: bad_width,
6304 ..Default::default()
6305 });
6306 assert!(
6307 err.contains("target_half_width"),
6308 "error must name the field, got: {err}"
6309 );
6310 }
6311
6312 let err = run(McConvergence {
6313 batch_size: 0,
6314 ..Default::default()
6315 });
6316 assert!(err.contains("batch_size"), "got: {err}");
6317
6318 let err = run(McConvergence {
6319 min_samples: 5_000,
6320 max_samples: 1_000,
6321 ..Default::default()
6322 });
6323 assert!(err.contains("max_samples"), "got: {err}");
6324 assert!(err.contains("min_samples"), "got: {err}");
6325
6326 let err = run(McConvergence {
6327 max_samples: 0,
6328 min_samples: 0,
6329 ..Default::default()
6330 });
6331 assert!(err.contains("max_samples"), "got: {err}");
6332
6333 let err = McConvergence {
6335 batch_size: 0,
6336 ..Default::default()
6337 }
6338 .validate()
6339 .unwrap_err();
6340 assert!(err.contains("batch_size"));
6341 }
6342
6343 #[test]
6344 fn wilson_companion_matches_hit_probability_and_wilson_interval() {
6345 let (inputs, wind) = seeded_test_fixture();
6346 let params = MonteCarloParams {
6347 num_simulations: 128,
6348 target_distance: Some(500.0),
6349 ..MonteCarloParams::default()
6350 };
6351 let results = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
6352 inputs,
6353 wind,
6354 params,
6355 0.01,
6356 0x1352_C0DE,
6357 )
6358 .expect("seeded legacy run");
6359
6360 for level in [
6361 ConfidenceLevel::P90,
6362 ConfidenceLevel::P95,
6363 ConfidenceLevel::P99,
6364 ] {
6365 let (p_hat, (lo, hi), n) = results.hit_probability_wilson(DEFAULT_HIT_RADIUS_M, level);
6366
6367 assert_eq!(
6370 p_hat.to_bits(),
6371 results.hit_probability(DEFAULT_HIT_RADIUS_M).to_bits()
6372 );
6373 assert_eq!(n, results.impact_positions.len() as u64);
6374
6375 let hits = results
6379 .impact_positions
6380 .iter()
6381 .filter(|p| MonteCarloResults::position_is_hit(p, DEFAULT_HIT_RADIUS_M))
6382 .count() as u64;
6383 let (want_lo, want_hi) = wilson_interval(hits, n, level);
6384 assert_eq!(lo.to_bits(), want_lo.to_bits());
6385 assert_eq!(hi.to_bits(), want_hi.to_bits());
6386 assert!(lo <= p_hat && p_hat <= hi, "interval excludes p_hat");
6387 }
6388
6389 let empty = MonteCarloResults {
6391 ranges: Vec::new(),
6392 impact_velocities: Vec::new(),
6393 impact_positions: Vec::new(),
6394 };
6395 assert_eq!(
6396 empty.hit_probability_wilson(DEFAULT_HIT_RADIUS_M, ConfidenceLevel::P95),
6397 (0.0, (0.0, 1.0), 0)
6398 );
6399 }
6400}
6401
6402#[cfg(test)]
6403mod monte_carlo_powder_curve_tests {
6404 use super::*;
6405 use rand::{rngs::StdRng, SeedableRng};
6406
6407 #[test]
6408 fn powder_curve_preserves_sampled_muzzle_velocity_dispersion() {
6409 let inputs = BallisticInputs {
6410 muzzle_velocity: 700.0,
6411 powder_temp_curve: Some(vec![(15.0, 800.0)]),
6412 powder_curve_temp_c: Some(15.0),
6413 ..BallisticInputs::default()
6414 };
6415 let params = MonteCarloParams {
6416 num_simulations: 16,
6417 velocity_std_dev: 20.0,
6418 angle_std_dev: 1e-12,
6419 bc_std_dev: 1e-12,
6420 wind_speed_std_dev: 1e-12,
6421 target_distance: Some(100.0),
6422 azimuth_std_dev: 1e-12,
6423 ..MonteCarloParams::default()
6424 };
6425
6426 let mut rng = StdRng::seed_from_u64(0x5EED_1176);
6427 let results = run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
6428 inputs,
6429 WindConditions::default(),
6430 params,
6431 0.0,
6432 &mut rng,
6433 )
6434 .expect("Monte Carlo solve");
6435 let min_velocity = results
6436 .impact_velocities
6437 .iter()
6438 .copied()
6439 .fold(f64::INFINITY, f64::min);
6440 let max_velocity = results
6441 .impact_velocities
6442 .iter()
6443 .copied()
6444 .fold(f64::NEG_INFINITY, f64::max);
6445
6446 assert!(
6447 max_velocity - min_velocity > 1.0,
6448 "20 m/s muzzle spread collapsed after curve resolution: impact-velocity span={} m/s",
6449 max_velocity - min_velocity
6450 );
6451 }
6452}
6453
6454#[cfg(test)]
6455mod monte_carlo_wind_sampling_tests {
6456 use super::*;
6457 use rand::{rngs::StdRng, SeedableRng};
6458
6459 #[test]
6460 fn wind_speed_sigma_does_not_change_seeded_direction_draws() {
6461 let base_wind = WindConditions {
6462 speed: 100.0,
6463 direction: 0.37,
6464 vertical_speed: 0.0,
6465 };
6466 let narrow_speed = MonteCarloWindSampler::new(&base_wind, 0.5, 0.2).unwrap();
6467 let wide_speed = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
6468 let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
6469 let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
6470 let mut speed_changed = false;
6471
6472 for _ in 0..32 {
6473 let narrow = narrow_speed.sample(&mut narrow_rng);
6474 let wide = wide_speed.sample(&mut wide_rng);
6475 assert!(narrow.speed > 0.0 && wide.speed > 0.0);
6476 assert_eq!(narrow.direction.to_bits(), wide.direction.to_bits());
6477 speed_changed |= narrow.speed.to_bits() != wide.speed.to_bits();
6478 }
6479 assert!(
6480 speed_changed,
6481 "different speed sigmas must still vary speed draws"
6482 );
6483 }
6484
6485 #[test]
6486 fn zero_direction_sigma_has_no_angular_jitter() {
6487 let base_wind = WindConditions {
6488 speed: 100.0,
6489 direction: 0.37,
6490 vertical_speed: 0.0,
6491 };
6492 let sampler = MonteCarloWindSampler::new(&base_wind, 4.0, 0.0).unwrap();
6493 let mut rng = StdRng::seed_from_u64(0x5EED_1223);
6494 let mut speed_changed = false;
6495
6496 for _ in 0..32 {
6497 let wind = sampler.sample(&mut rng);
6498 speed_changed |= wind.speed.to_bits() != base_wind.speed.to_bits();
6499 assert_eq!(wind.direction.to_bits(), base_wind.direction.to_bits());
6500 }
6501 assert!(speed_changed, "speed uncertainty should remain active");
6502 }
6503
6504 #[test]
6505 fn direction_sigma_controls_seeded_angular_spread_in_radians() {
6506 let base_wind = WindConditions {
6507 speed: 100.0,
6508 direction: 0.37,
6509 vertical_speed: 0.0,
6510 };
6511 let narrow = MonteCarloWindSampler::new(&base_wind, 4.0, 0.1).unwrap();
6512 let wide = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
6513 let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
6514 let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
6515 let mut nonzero_direction_draw = false;
6516
6517 for _ in 0..32 {
6518 let narrow_wind = narrow.sample(&mut narrow_rng);
6519 let wide_wind = wide.sample(&mut wide_rng);
6520 assert_eq!(narrow_wind.speed.to_bits(), wide_wind.speed.to_bits());
6521
6522 let narrow_delta = narrow_wind.direction - base_wind.direction;
6523 let wide_delta = wide_wind.direction - base_wind.direction;
6524 assert!((wide_delta - 2.0 * narrow_delta).abs() < 1e-12);
6525 nonzero_direction_draw |= narrow_delta.abs() > 1e-6;
6526 }
6527 assert!(
6528 nonzero_direction_draw,
6529 "positive radians sigma must vary direction"
6530 );
6531 }
6532
6533 #[test]
6534 fn direction_sigma_rejects_negative_or_nonfinite_values() {
6535 let base_wind = WindConditions::default();
6536 for sigma in [-0.1, f64::NAN, f64::INFINITY] {
6537 assert!(MonteCarloWindSampler::new(&base_wind, 1.0, sigma).is_err());
6538 }
6539 }
6540
6541 #[test]
6542 fn base_vertical_wind_rides_into_every_mc_sample() {
6543 use rand::SeedableRng;
6547 let base_wind = WindConditions {
6548 vertical_speed: 4.2,
6549 ..Default::default()
6550 };
6551 let sampler = MonteCarloWindSampler::new(&base_wind, 1.0, 0.2).unwrap();
6552 let mut rng = rand::rngs::StdRng::seed_from_u64(7);
6553 for _ in 0..32 {
6554 let w = sampler.sample(&mut rng);
6555 assert_eq!(w.vertical_speed, 4.2);
6556 }
6557 }
6558
6559 #[test]
6560 fn negative_speed_sample_reverses_wind_direction() {
6561 let direction = 0.25;
6562 let signed_speed = -2.5;
6563 let wind = wind_from_signed_speed_sample(signed_speed, direction, 0.0);
6564 let positive_wind = wind_from_signed_speed_sample(2.5, direction, 0.0);
6565
6566 assert_eq!(wind.speed, 2.5);
6567 assert!(
6568 (wind.direction - (direction + std::f64::consts::PI)).abs() < f64::EPSILON,
6569 "negative speed must reverse direction by pi: got {}",
6570 wind.direction
6571 );
6572 assert_eq!(positive_wind.speed, 2.5);
6573 assert_eq!(positive_wind.direction, direction);
6574
6575 let normalized_x = -wind.speed * wind.direction.cos();
6576 let normalized_z = -wind.speed * wind.direction.sin();
6577 let signed_x = -signed_speed * direction.cos();
6578 let signed_z = -signed_speed * direction.sin();
6579 assert!((normalized_x - signed_x).abs() < 1e-12);
6580 assert!((normalized_z - signed_z).abs() < 1e-12);
6581 }
6582}
6583
6584#[cfg(test)]
6585mod bc_fit_objective_tests {
6586 use super::*;
6587
6588 fn velocity_point(range_m: f64, velocity_mps: f64) -> TrajectoryPoint {
6589 TrajectoryPoint {
6590 time: 0.0,
6591 position: Vector3::new(range_m, 0.0, 0.0),
6592 velocity_magnitude: velocity_mps,
6593 kinetic_energy: 0.0,
6594 drag_coefficient: None,
6595 }
6596 }
6597
6598 #[test]
6599 fn candidate_that_misses_an_observation_has_no_score() {
6600 let trajectory = vec![velocity_point(0.0, 800.0), velocity_point(100.0, 700.0)];
6601 let observations = vec![(50.0, 750.0), (150.0, 600.0)];
6602
6603 assert!(
6604 fit_residual_sse(&trajectory, &observations, BcFitMode::Velocity, 0.0).is_none(),
6605 "a candidate that reaches only one of two observations must not compete on partial SSE"
6606 );
6607
6608 let complete_observations = vec![(50.0, 740.0), (100.0, 680.0)];
6609 assert_eq!(
6610 fit_residual_sse(
6611 &trajectory,
6612 &complete_observations,
6613 BcFitMode::Velocity,
6614 0.0,
6615 ),
6616 Some(500.0)
6617 );
6618 }
6619}
6620
6621#[cfg(test)]
6622mod cluster_bc_reference_space_tests {
6623 use super::*;
6624
6625 fn acceleration_at_1100_fps(inputs: BallisticInputs) -> Vector3<f64> {
6626 let solver = TrajectorySolver::new(
6627 inputs,
6628 WindConditions::default(),
6629 AtmosphericConditions::default(),
6630 );
6631 let position = Vector3::zeros();
6632 let velocity = Vector3::new(1100.0 / 3.28084, 0.0, 0.0);
6633 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6634 solver.calculate_acceleration(
6635 &position,
6636 &velocity,
6637 &Vector3::zeros(),
6638 (temp_c, pressure_hpa, density / 1.225),
6639 )
6640 }
6641
6642 #[test]
6643 fn solver_passes_g7_reference_model_to_cluster_classifier() {
6644 let inputs = BallisticInputs {
6645 bc_value: 0.190,
6646 bc_type: DragModel::G7,
6647 bullet_mass: 77.0 * crate::constants::GRAINS_TO_KG,
6648 bullet_diameter: 0.224 * 0.0254,
6649 use_cluster_bc: true,
6650 ..BallisticInputs::default()
6651 };
6652
6653 let solver = TrajectorySolver::new(
6654 inputs,
6655 WindConditions::default(),
6656 AtmosphericConditions::default(),
6657 );
6658 let corrected = solver.apply_cluster_bc_correction(0.190, 2800.0);
6659
6660 assert!(
6661 (corrected / 0.190 - 1.004).abs() < 1e-12,
6662 "solver selected the wrong G7 cluster multiplier: {}",
6663 corrected / 0.190
6664 );
6665 }
6666
6667 #[test]
6668 fn velocity_bc_segments_are_not_cluster_corrected_twice() {
6669 let segmented_clustered = BallisticInputs {
6670 bc_value: 0.5,
6671 bc_type: DragModel::G7,
6672 use_bc_segments: true,
6673 bc_segments_data: Some(vec![
6674 crate::BCSegmentData {
6675 velocity_min: 0.0,
6676 velocity_max: 1_600.0,
6677 bc_value: 0.4,
6678 },
6679 crate::BCSegmentData {
6680 velocity_min: 1_600.0,
6681 velocity_max: 5_000.0,
6682 bc_value: 0.45,
6683 },
6684 ]),
6685 use_cluster_bc: true,
6686 ..BallisticInputs::default()
6687 };
6688 let mut segmented_only = segmented_clustered.clone();
6689 segmented_only.use_cluster_bc = false;
6690 let mut constant_clustered = segmented_clustered.clone();
6691 constant_clustered.bc_value = 0.4;
6692 constant_clustered.bc_segments_data = None;
6693
6694 let stacked = acceleration_at_1100_fps(segmented_clustered);
6695 let segment_only = acceleration_at_1100_fps(segmented_only);
6696 let cluster_only = acceleration_at_1100_fps(constant_clustered);
6697
6698 assert!(
6699 (stacked.x - segment_only.x).abs() < 1e-12,
6700 "segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
6701 stacked.x,
6702 segment_only.x
6703 );
6704 assert!(
6705 (cluster_only.x - segment_only.x).abs() > 1e-6,
6706 "cluster correction must remain active for a constant BC"
6707 );
6708 }
6709
6710 #[test]
6711 fn mach_bc_segments_are_not_cluster_corrected_twice() {
6712 let mach_segmented_clustered = BallisticInputs {
6713 bc_value: 0.5,
6714 bc_type: DragModel::G7,
6715 use_bc_segments: false,
6716 bc_segments: Some(vec![(0.5, 0.3), (1.5, 0.5)]),
6717 use_cluster_bc: true,
6718 ..BallisticInputs::default()
6719 };
6720 let mut mach_segmented_only = mach_segmented_clustered.clone();
6721 mach_segmented_only.use_cluster_bc = false;
6722
6723 let stacked = acceleration_at_1100_fps(mach_segmented_clustered);
6724 let segment_only = acceleration_at_1100_fps(mach_segmented_only);
6725
6726 assert!(
6727 (stacked.x - segment_only.x).abs() < 1e-12,
6728 "Mach segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
6729 stacked.x,
6730 segment_only.x
6731 );
6732 }
6733}
6734
6735#[cfg(test)]
6736mod velocity_bc_flag_tests {
6737 use super::*;
6738
6739 fn acceleration_at_600_mps(inputs: BallisticInputs) -> Vector3<f64> {
6740 let solver = TrajectorySolver::new(
6741 inputs,
6742 WindConditions::default(),
6743 AtmosphericConditions::default(),
6744 );
6745 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6746 solver.calculate_acceleration(
6747 &Vector3::zeros(),
6748 &Vector3::new(600.0, 0.0, 0.0),
6749 &Vector3::zeros(),
6750 (temp_c, pressure_hpa, density / 1.225),
6751 )
6752 }
6753
6754 #[test]
6755 fn velocity_bc_data_requires_opt_in_in_trajectory_solver() {
6756 let scalar_inputs = BallisticInputs {
6757 bc_value: 0.5,
6758 bc_type: DragModel::G7,
6759 ..BallisticInputs::default()
6760 };
6761 let mut disabled_inputs = scalar_inputs.clone();
6762 disabled_inputs.bc_segments_data = Some(vec![crate::BCSegmentData {
6763 velocity_min: 0.0,
6764 velocity_max: 4_000.0,
6765 bc_value: 0.46,
6766 }]);
6767 disabled_inputs.use_bc_segments = false;
6768 let mut enabled_inputs = disabled_inputs.clone();
6769 enabled_inputs.use_bc_segments = true;
6770 let mut mach_only_inputs = scalar_inputs.clone();
6771 mach_only_inputs.bc_segments = Some(vec![(0.0, 0.4), (3.0, 0.4)]);
6772 let mut disabled_with_both = mach_only_inputs.clone();
6773 disabled_with_both.bc_segments_data = disabled_inputs.bc_segments_data.clone();
6774
6775 let scalar = acceleration_at_600_mps(scalar_inputs);
6776 let disabled = acceleration_at_600_mps(disabled_inputs);
6777 let enabled = acceleration_at_600_mps(enabled_inputs);
6778 let mach_only = acceleration_at_600_mps(mach_only_inputs);
6779 let disabled_with_both = acceleration_at_600_mps(disabled_with_both);
6780
6781 assert_eq!(
6782 disabled.x.to_bits(),
6783 scalar.x.to_bits(),
6784 "a populated velocity table must not change drag while use_bc_segments is false"
6785 );
6786 assert!(
6787 enabled.x < disabled.x - 1.0,
6788 "enabling the lower BC table must increase drag: disabled ax={} enabled ax={}",
6789 disabled.x,
6790 enabled.x
6791 );
6792 assert_eq!(
6793 disabled_with_both.x.to_bits(),
6794 mach_only.x.to_bits(),
6795 "disabling velocity data must fall through to an explicit Mach table"
6796 );
6797 }
6798}
6799
6800#[cfg(test)]
6801mod mach_bc_segment_tests {
6802 use super::*;
6803
6804 #[test]
6805 fn trajectory_solver_interpolates_explicit_mach_bc_segments() {
6806 let segmented_inputs = BallisticInputs {
6807 bc_value: 0.8,
6808 use_bc_segments: false,
6809 bc_segments: Some(vec![(1.0, 0.2), (2.0, 0.4)]),
6810 bc_segments_data: None,
6811 ..BallisticInputs::default()
6812 };
6813
6814 let mut expected_inputs = segmented_inputs.clone();
6815 expected_inputs.bc_value = 0.3;
6816 expected_inputs.bc_segments = None;
6817
6818 let atmosphere = AtmosphericConditions::default();
6819 let segmented_solver = TrajectorySolver::new(
6820 segmented_inputs,
6821 WindConditions::default(),
6822 atmosphere.clone(),
6823 );
6824 let expected_solver =
6825 TrajectorySolver::new(expected_inputs, WindConditions::default(), atmosphere);
6826 let position = Vector3::zeros();
6827 let (density, _, temp_c, pressure_hpa) = segmented_solver.resolved_atmosphere();
6828 let (_, local_speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
6829 segmented_solver.atmosphere.altitude,
6830 segmented_solver.atmosphere.altitude,
6831 temp_c,
6832 pressure_hpa,
6833 density / 1.225,
6834 segmented_solver.atmosphere.humidity,
6835 );
6836 let velocity = Vector3::new(1.5 * local_speed_of_sound, 0.0, 0.0);
6837 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
6838
6839 let segmented_acceleration = segmented_solver.calculate_acceleration(
6840 &position,
6841 &velocity,
6842 &Vector3::zeros(),
6843 resolved_atmo,
6844 );
6845 let expected_acceleration = expected_solver.calculate_acceleration(
6846 &position,
6847 &velocity,
6848 &Vector3::zeros(),
6849 resolved_atmo,
6850 );
6851
6852 assert!(
6853 (segmented_acceleration.x - expected_acceleration.x).abs() < 1e-12,
6854 "Mach 1.5 must interpolate BC 0.3: segmented ax={} expected ax={}",
6855 segmented_acceleration.x,
6856 expected_acceleration.x
6857 );
6858 }
6859}
6860
6861#[cfg(test)]
6862mod custom_drag_table_validation_tests {
6863 use super::*;
6864
6865 #[test]
6866 fn solve_accepts_zero_bc_when_custom_table_present() {
6867 let inputs = BallisticInputs {
6868 bc_value: 0.0, bullet_mass: 0.0106,
6870 bullet_diameter: 0.00782,
6871 muzzle_velocity: 850.0,
6872 custom_drag_table: Some(crate::drag::DragTable::new(
6873 vec![0.5, 1.0, 2.0, 3.0],
6874 vec![0.23, 0.40, 0.30, 0.26],
6875 )),
6876 ..BallisticInputs::default()
6877 };
6878 let solver = TrajectorySolver::new(
6879 inputs,
6880 WindConditions::default(),
6881 AtmosphericConditions::default(),
6882 );
6883 assert!(solver.solve().is_ok());
6885 }
6886
6887 #[test]
6888 fn solve_still_requires_bc_without_table() {
6889 let inputs = BallisticInputs {
6890 bc_value: 0.0,
6891 bullet_mass: 0.0106,
6892 bullet_diameter: 0.00782,
6893 muzzle_velocity: 850.0,
6894 ..BallisticInputs::default()
6895 };
6896 let solver = TrajectorySolver::new(
6897 inputs,
6898 WindConditions::default(),
6899 AtmosphericConditions::default(),
6900 );
6901 assert!(solver.solve().is_err());
6902 }
6903}
6904
6905#[cfg(test)]
6907mod cd_scale_tests {
6908 use super::*;
6909
6910 fn deck() -> crate::drag::DragTable {
6911 crate::drag::DragTable::new(vec![0.5, 1.0, 2.0, 3.0], vec![0.23, 0.40, 0.30, 0.26])
6912 }
6913
6914 fn deck_inputs(cd_scale: f64) -> BallisticInputs {
6915 BallisticInputs {
6916 bullet_mass: 0.0106,
6917 bullet_diameter: 0.00782,
6918 muzzle_velocity: 850.0,
6919 custom_drag_table: Some(deck()),
6920 cd_scale,
6921 ..BallisticInputs::default()
6922 }
6923 }
6924
6925 #[test]
6926 fn default_cd_scale_is_one() {
6927 assert_eq!(BallisticInputs::default().cd_scale, 1.0);
6928 }
6929
6930 #[test]
6933 fn cd_scale_absent_is_byte_identical_to_explicit_one() {
6934 let omitted = BallisticInputs {
6935 bullet_mass: 0.0106,
6936 bullet_diameter: 0.00782,
6937 muzzle_velocity: 850.0,
6938 custom_drag_table: Some(deck()),
6939 ..BallisticInputs::default()
6940 };
6941 let explicit = BallisticInputs {
6942 cd_scale: 1.0,
6943 ..omitted.clone()
6944 };
6945
6946 let solver_omitted = TrajectorySolver::new(
6947 omitted,
6948 WindConditions::default(),
6949 AtmosphericConditions::default(),
6950 );
6951 let solver_explicit = TrajectorySolver::new(
6952 explicit,
6953 WindConditions::default(),
6954 AtmosphericConditions::default(),
6955 );
6956
6957 let cd_omitted = solver_omitted.calculate_drag_coefficient(700.0, 340.0);
6958 let cd_explicit = solver_explicit.calculate_drag_coefficient(700.0, 340.0);
6959 assert_eq!(
6960 cd_omitted.to_bits(),
6961 cd_explicit.to_bits(),
6962 "default cd_scale must be bit-identical to an explicit 1.0"
6963 );
6964
6965 let result = solver_omitted.solve();
6968 assert!(
6969 result.is_ok(),
6970 "existing custom-deck solves must pass unchanged"
6971 );
6972 }
6973
6974 #[test]
6976 fn cd_scale_multiplies_the_interpolated_cd_exactly() {
6977 let velocity = 700.0;
6978 let speed_of_sound = 340.0;
6979 let mach = velocity / speed_of_sound;
6980 let expected_unscaled = deck().interpolate(mach);
6981
6982 for &scale in &[0.90, 1.0, 1.10, 1.5] {
6983 let solver = TrajectorySolver::new(
6984 deck_inputs(scale),
6985 WindConditions::default(),
6986 AtmosphericConditions::default(),
6987 );
6988 let cd = solver.calculate_drag_coefficient(velocity, speed_of_sound);
6989 assert!(
6990 (cd - expected_unscaled * scale).abs() < 1e-12,
6991 "scale={scale}: cd={cd} expected={}",
6992 expected_unscaled * scale
6993 );
6994 }
6995 }
6996
6997 #[test]
7001 fn cd_scale_direction_on_cli_api_solver() {
7002 let solve = |scale: f64| {
7003 TrajectorySolver::new(
7004 deck_inputs(scale),
7005 WindConditions::default(),
7006 AtmosphericConditions::default(),
7007 )
7008 .solve()
7009 .expect("custom-deck solve should succeed")
7010 };
7011
7012 let baseline = solve(1.0);
7013 let scaled_up = solve(1.10);
7014 let scaled_down = solve(0.90);
7015
7016 assert!(
7017 scaled_up.impact_velocity < baseline.impact_velocity,
7018 "cd_scale=1.10 must increase drag -> lower impact velocity: base={} up={}",
7019 baseline.impact_velocity,
7020 scaled_up.impact_velocity
7021 );
7022 assert!(
7023 scaled_down.impact_velocity > baseline.impact_velocity,
7024 "cd_scale=0.90 must decrease drag -> higher impact velocity: base={} down={}",
7025 baseline.impact_velocity,
7026 scaled_down.impact_velocity
7027 );
7028 }
7029
7030 #[test]
7032 fn validate_for_solve_rejects_invalid_cd_scale() {
7033 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
7034 let solver = TrajectorySolver::new(
7035 deck_inputs(bad),
7036 WindConditions::default(),
7037 AtmosphericConditions::default(),
7038 );
7039 assert!(
7040 solver.solve().is_err(),
7041 "cd_scale={bad} must be rejected by validate_for_solve"
7042 );
7043 }
7044 }
7045
7046 #[test]
7053 fn validate_for_solve_rejects_invalid_cd_scale_without_a_custom_drag_table() {
7054 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
7055 let inputs = BallisticInputs {
7056 bc_value: 0.5,
7057 bc_type: crate::DragModel::G1,
7058 bullet_mass: 0.0106,
7059 bullet_diameter: 0.00782,
7060 muzzle_velocity: 850.0,
7061 cd_scale: bad,
7062 ..BallisticInputs::default()
7063 };
7064 assert!(
7065 inputs.custom_drag_table.is_none(),
7066 "precondition: no custom deck"
7067 );
7068 let solver = TrajectorySolver::new(
7069 inputs,
7070 WindConditions::default(),
7071 AtmosphericConditions::default(),
7072 );
7073 assert!(
7074 solver.solve().is_err(),
7075 "cd_scale={bad} must be rejected by validate_for_solve even without a custom \
7076 drag table"
7077 );
7078 }
7079 }
7080
7081 #[test]
7085 fn cd_scale_is_inert_without_a_custom_drag_table() {
7086 let make = |cd_scale: f64| BallisticInputs {
7087 bc_value: 0.5,
7088 bc_type: crate::DragModel::G1,
7089 bullet_mass: 0.0106,
7090 bullet_diameter: 0.00782,
7091 muzzle_velocity: 850.0,
7092 cd_scale,
7093 ..BallisticInputs::default()
7094 };
7095 let solver_neutral = TrajectorySolver::new(
7096 make(1.0),
7097 WindConditions::default(),
7098 AtmosphericConditions::default(),
7099 );
7100 let solver_far = TrajectorySolver::new(
7101 make(1.5),
7102 WindConditions::default(),
7103 AtmosphericConditions::default(),
7104 );
7105 let cd_neutral = solver_neutral.calculate_drag_coefficient(700.0, 340.0);
7106 let cd_far = solver_far.calculate_drag_coefficient(700.0, 340.0);
7107 assert_eq!(
7108 cd_neutral.to_bits(),
7109 cd_far.to_bits(),
7110 "cd_scale must not affect the G-model/BC drag path"
7111 );
7112 }
7113
7114 #[test]
7117 fn cd_scale_shifts_all_three_solver_paths_in_the_same_direction() {
7118 let cli_solve = |scale: f64| {
7120 TrajectorySolver::new(
7121 deck_inputs(scale),
7122 WindConditions::default(),
7123 AtmosphericConditions::default(),
7124 )
7125 .solve()
7126 .expect("cli_api custom-deck solve should succeed")
7127 };
7128 let cli_baseline = cli_solve(1.0);
7129 let cli_scaled = cli_solve(1.10);
7130 assert!(
7131 cli_scaled.impact_velocity < cli_baseline.impact_velocity,
7132 "cli_api: cd_scale=1.10 must lower impact velocity"
7133 );
7134
7135 let derivatives_accel_x = |scale: f64| {
7137 let inputs = deck_inputs(scale);
7138 crate::derivatives::compute_derivatives(
7139 nalgebra::Vector3::zeros(),
7140 nalgebra::Vector3::new(700.0, 0.0, 0.0),
7141 &inputs,
7142 nalgebra::Vector3::zeros(),
7143 (1.225, 340.0, 0.0, 0.0),
7144 inputs.bc_value,
7145 None,
7146 0.0,
7147 None,
7148 )[3]
7149 };
7150 let deriv_baseline = derivatives_accel_x(1.0);
7151 let deriv_scaled = derivatives_accel_x(1.10);
7152 assert!(
7153 deriv_scaled < deriv_baseline,
7154 "derivatives: cd_scale=1.10 must make x-acceleration more negative (more drag): \
7155 base={deriv_baseline} scaled={deriv_scaled}"
7156 );
7157
7158 let fast_final_speed = |scale: f64| {
7160 let inputs = deck_inputs(scale);
7161 let wind_sock = crate::wind::WindSock::new(vec![]);
7162 let params = crate::fast_trajectory::FastIntegrationParams {
7163 horiz: 500.0,
7164 vert: 0.0,
7165 initial_state: [0.0, 0.0, 0.0, 850.0, 0.0, 0.0],
7166 t_span: (0.0, 5.0),
7167 atmo_params: (0.0, 15.0, 1013.25, 1.0),
7168 atmo_sock: None,
7169 };
7170 let solution = crate::fast_trajectory::fast_integrate(&inputs, &wind_sock, params);
7171 assert!(
7172 solution.success,
7173 "fast_integrate must succeed for scale={scale}"
7174 );
7175 let last = solution.t.len() - 1;
7176 let (vx, vy, vz) = (
7177 solution.y[3][last],
7178 solution.y[4][last],
7179 solution.y[5][last],
7180 );
7181 (vx * vx + vy * vy + vz * vz).sqrt()
7182 };
7183 let fast_baseline = fast_final_speed(1.0);
7184 let fast_scaled = fast_final_speed(1.10);
7185 assert!(
7186 fast_scaled < fast_baseline,
7187 "fast_trajectory: cd_scale=1.10 must lower final speed: base={fast_baseline} scaled={fast_scaled}"
7188 );
7189 }
7190}
7191
7192#[cfg(test)]
7193mod humid_local_mach_tests {
7194 use super::*;
7195
7196 fn solver_with_station_humidity(humidity_percent: f64) -> TrajectorySolver {
7197 let inputs = BallisticInputs {
7198 custom_drag_table: Some(crate::drag::DragTable::new(vec![0.5, 1.5], vec![0.1, 1.1])),
7199 ..BallisticInputs::default()
7200 };
7201 TrajectorySolver::new(
7202 inputs,
7203 WindConditions::default(),
7204 AtmosphericConditions {
7205 temperature: 30.0,
7206 pressure: 1013.25,
7207 humidity: humidity_percent,
7208 altitude: 0.0,
7209 },
7210 )
7211 }
7212
7213 fn acceleration(solver: &TrajectorySolver, base_ratio: f64) -> Vector3<f64> {
7214 solver.calculate_acceleration(
7215 &Vector3::zeros(),
7216 &Vector3::new(350.0, 0.0, 0.0),
7217 &Vector3::zeros(),
7218 (30.0, 1013.25, base_ratio),
7219 )
7220 }
7221
7222 #[test]
7223 fn local_mach_uses_station_humidity_when_density_is_held_constant() {
7224 let dry = acceleration(&solver_with_station_humidity(0.0), 1.0);
7225 let humid = acceleration(&solver_with_station_humidity(100.0), 1.0);
7226
7227 assert!(
7228 humid.x > dry.x,
7229 "humid sound speed should lower Mach and drag on the rising test curve: dry ax={} humid ax={}",
7230 dry.x,
7231 humid.x
7232 );
7233 }
7234
7235 #[test]
7236 fn active_atmosphere_zone_uses_zone_humidity_instead_of_station_humidity() {
7237 let zone_humidity = 80.0;
7238 let zone_ratio =
7239 crate::atmosphere::calculate_air_density_cimp(30.0, 1013.25, zone_humidity) / 1.225;
7240 let station_solver = solver_with_station_humidity(zone_humidity);
7241 let mut zoned_solver = solver_with_station_humidity(0.0);
7242 zoned_solver.set_atmo_segments(vec![(30.0, 1013.25, zone_humidity, 1_000.0)]);
7243
7244 let station = acceleration(&station_solver, zone_ratio);
7245 let zoned = acceleration(&zoned_solver, zone_ratio);
7246
7247 assert!(
7248 (zoned - station).norm() < 1e-12,
7249 "active zone T/P/RH should override the station atmosphere: station={station:?} zoned={zoned:?}"
7250 );
7251 }
7252}
7253
7254#[cfg(test)]
7255mod inclined_atmosphere_frame_tests {
7256 use super::*;
7257
7258 fn expected_shot_frame_vector(level: Vector3<f64>, angle: f64) -> Vector3<f64> {
7259 let (sin_angle, cos_angle) = angle.sin_cos();
7260 Vector3::new(
7261 level.x * cos_angle + level.y * sin_angle,
7262 -level.x * sin_angle + level.y * cos_angle,
7263 level.z,
7264 )
7265 }
7266
7267 #[test]
7268 fn inclined_positions_at_same_world_altitude_have_same_solver_acceleration() {
7269 let angle = std::f64::consts::FRAC_PI_6;
7270 let inputs = BallisticInputs {
7271 shooting_angle: angle,
7272 ..BallisticInputs::default()
7273 };
7274 let atmosphere = AtmosphericConditions {
7275 altitude: 100.0,
7276 ..AtmosphericConditions::default()
7277 };
7278 let solver = TrajectorySolver::new(inputs, WindConditions::default(), atmosphere);
7279 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7280 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
7281 let velocity = Vector3::new(600.0, 0.0, 0.0);
7282 let along_slant = Vector3::new(1_000.0, 0.0, 0.0);
7283 let across_slant = Vector3::new(0.0, 500.0 / angle.cos(), 0.0);
7284
7285 let a = solver.calculate_acceleration(
7286 &along_slant,
7287 &velocity,
7288 &Vector3::zeros(),
7289 resolved_atmo,
7290 );
7291 let b = solver.calculate_acceleration(
7292 &across_slant,
7293 &velocity,
7294 &Vector3::zeros(),
7295 resolved_atmo,
7296 );
7297
7298 assert!(
7299 (a - b).norm() < 1e-10,
7300 "solver acceleration differs at equal world altitude: {a:?} vs {b:?}"
7301 );
7302 }
7303
7304 #[test]
7305 fn inclined_headwind_is_rotated_into_solver_frame() {
7306 let angle = std::f64::consts::FRAC_PI_6;
7307 let inputs = BallisticInputs {
7308 shooting_angle: angle,
7309 ..BallisticInputs::default()
7310 };
7311 let solver = TrajectorySolver::new(
7312 inputs,
7313 WindConditions::default(),
7314 AtmosphericConditions::default(),
7315 );
7316 let level_headwind = Vector3::new(-100.0, 0.0, 0.0);
7317 let velocity = expected_shot_frame_vector(level_headwind, angle);
7318 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7319 let actual = solver.calculate_acceleration(
7320 &Vector3::zeros(),
7321 &velocity,
7322 &level_headwind,
7323 (temp_c, pressure_hpa, density / 1.225),
7324 );
7325
7326 assert!(
7327 (actual - solver.gravity_acceleration()).norm() < 1e-12,
7328 "co-moving horizontal wind must leave only shot-frame gravity: {actual:?}"
7329 );
7330 }
7331
7332 #[test]
7333 fn inclined_coriolis_is_rotated_into_solver_frame() {
7334 let angle = std::f64::consts::FRAC_PI_6;
7335 let latitude_deg = 45.0_f64;
7336 let shot_azimuth = 0.4_f64;
7337 let velocity = Vector3::new(600.0, 20.0, 5.0);
7338 let base_inputs = BallisticInputs {
7339 shooting_angle: angle,
7340 latitude: Some(latitude_deg),
7341 shot_azimuth,
7342 ..BallisticInputs::default()
7343 };
7344 let acceleration = |enable_coriolis| {
7345 let mut inputs = base_inputs.clone();
7346 inputs.enable_coriolis = enable_coriolis;
7347 let solver = TrajectorySolver::new(
7348 inputs,
7349 WindConditions::default(),
7350 AtmosphericConditions::default(),
7351 );
7352 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7353 solver.calculate_acceleration(
7354 &Vector3::zeros(),
7355 &velocity,
7356 &Vector3::zeros(),
7357 (temp_c, pressure_hpa, density / 1.225),
7358 )
7359 };
7360
7361 let omega_earth = 7.2921159e-5_f64;
7362 let latitude = latitude_deg.to_radians();
7363 let level_omega = Vector3::new(
7364 omega_earth * latitude.cos() * shot_azimuth.cos(),
7365 omega_earth * latitude.sin(),
7366 -omega_earth * latitude.cos() * shot_azimuth.sin(),
7367 );
7368 let expected = -2.0 * expected_shot_frame_vector(level_omega, angle).cross(&velocity);
7369 let actual = acceleration(true) - acceleration(false);
7370
7371 assert!(
7372 (actual - expected).norm() < 1e-12,
7373 "inclined Coriolis mismatch: actual={actual:?}, expected={expected:?}"
7374 );
7375 }
7376}
7377
7378#[cfg(test)]
7379mod terminal_range_interpolation_tests {
7380 use super::*;
7381
7382 #[test]
7383 fn terminal_finalizer_selects_the_earliest_crossed_boundary() {
7384 let inputs = BallisticInputs {
7385 ground_threshold: 0.0,
7386 ..BallisticInputs::default()
7387 };
7388 let mut solver = TrajectorySolver::new(
7389 inputs,
7390 WindConditions::default(),
7391 AtmosphericConditions::default(),
7392 );
7393 solver.set_max_range(120.0);
7394
7395 let previous_speed = 700.0;
7396 let mut points = vec![TrajectoryPoint {
7397 time: 99.0,
7398 position: Vector3::new(90.0, 1.0, -1.0),
7399 velocity_magnitude: previous_speed,
7400 kinetic_energy: 0.5 * solver.inputs.bullet_mass * previous_speed.powi(2),
7401 drag_coefficient: None,
7402 }];
7403 let mut max_height = 1.0;
7404 let termination = solver
7405 .append_terminal_endpoint(
7406 &mut points,
7407 Vector3::new(130.0, -3.0, 3.0),
7408 Vector3::new(600.0, 0.0, 0.0),
7409 101.0,
7410 &mut max_height,
7411 )
7412 .expect("the final step brackets supported boundaries");
7413
7414 assert_eq!(termination, TrajectoryTermination::GroundThreshold);
7415 assert_eq!(points.len(), 2);
7416 let terminal = points.last().expect("terminal point");
7417 assert_eq!(terminal.time, 99.5);
7418 assert_eq!(terminal.position, Vector3::new(100.0, 0.0, 0.0));
7419 assert_eq!(terminal.velocity_magnitude, 675.0);
7420 assert_eq!(
7421 terminal.kinetic_energy,
7422 0.5 * solver.inputs.bullet_mass * 675.0_f64.powi(2)
7423 );
7424
7425 solver.set_max_range(100.0);
7427 let mut tied_points = vec![points[0].clone()];
7428 assert_eq!(
7429 solver
7430 .append_terminal_endpoint(
7431 &mut tied_points,
7432 Vector3::new(130.0, -3.0, 3.0),
7433 Vector3::new(600.0, 0.0, 0.0),
7434 101.0,
7435 &mut max_height,
7436 )
7437 .expect("tied boundaries remain a valid terminal"),
7438 TrajectoryTermination::GroundThreshold
7439 );
7440 }
7441
7442 #[test]
7443 fn sub_ulp_terminal_crossing_replaces_instead_of_duplicating_range() {
7444 let ground_threshold = f64::from_bits(1.0_f64.to_bits() - 1);
7445 let inputs = BallisticInputs {
7446 ground_threshold,
7447 ..BallisticInputs::default()
7448 };
7449 let mut solver = TrajectorySolver::new(
7450 inputs,
7451 WindConditions::default(),
7452 AtmosphericConditions::default(),
7453 );
7454 solver.set_max_range(1_000.0);
7455
7456 let speed = 700.0;
7457 let mut points = vec![TrajectoryPoint {
7458 time: 0.0,
7459 position: Vector3::new(100.0, 1.0, 0.0),
7460 velocity_magnitude: speed,
7461 kinetic_energy: 0.5 * solver.inputs.bullet_mass * speed.powi(2),
7462 drag_coefficient: None,
7463 }];
7464 let mut max_height = 1.0;
7465 let termination = solver
7466 .append_terminal_endpoint(
7467 &mut points,
7468 Vector3::new(101.0, 0.0, 0.0),
7469 Vector3::new(699.0, 0.0, 0.0),
7470 1.0,
7471 &mut max_height,
7472 )
7473 .expect("sub-ULP ground crossing remains representable as one terminal state");
7474
7475 assert_eq!(termination, TrajectoryTermination::GroundThreshold);
7476 assert_eq!(points.len(), 1);
7477 assert_eq!(points[0].position.x, 100.0);
7478 assert_eq!(points[0].position.y.to_bits(), ground_threshold.to_bits());
7479 assert!(points[0].time > 0.0);
7480 }
7481
7482 #[test]
7483 fn every_solver_appends_an_exact_max_range_endpoint() {
7484 let target_range = 0.1;
7485 let modes = [
7486 ("Euler", false, false),
7487 ("RK4", true, false),
7488 ("RK45", true, true),
7489 ];
7490
7491 for (name, use_rk4, use_adaptive_rk45) in modes {
7492 let inputs = BallisticInputs {
7493 use_rk4,
7494 use_adaptive_rk45,
7495 ground_threshold: f64::NEG_INFINITY,
7496 enable_trajectory_sampling: true,
7497 sample_interval: target_range,
7498 ..BallisticInputs::default()
7499 };
7500 let mut solver = TrajectorySolver::new(
7501 inputs,
7502 WindConditions::default(),
7503 AtmosphericConditions::default(),
7504 );
7505 solver.set_max_range(target_range);
7506
7507 let result = solver.solve().expect("short-range solve should succeed");
7508 let terminal = result.points.last().expect("terminal point is missing");
7509 let muzzle = result.points.first().expect("muzzle point is missing");
7510
7511 assert_eq!(result.termination, TrajectoryTermination::MaxRange);
7512 assert_eq!(
7513 terminal.position.x.to_bits(),
7514 target_range.to_bits(),
7515 "{name} did not terminate exactly at max_range"
7516 );
7517 assert_eq!(result.max_range.to_bits(), target_range.to_bits());
7518 assert!(
7519 result.time_of_flight > 0.0 && result.time_of_flight < solver.time_step,
7520 "{name} terminal time was not interpolated within the crossing step: {}",
7521 result.time_of_flight
7522 );
7523 assert_eq!(result.time_of_flight.to_bits(), terminal.time.to_bits());
7524 assert_eq!(
7525 result.impact_velocity.to_bits(),
7526 terminal.velocity_magnitude.to_bits()
7527 );
7528 assert_eq!(
7529 result.impact_energy.to_bits(),
7530 terminal.kinetic_energy.to_bits()
7531 );
7532 let expected_energy = 0.5 * solver.inputs.bullet_mass * result.impact_velocity.powi(2);
7533 assert!((result.impact_energy - expected_energy).abs() < 1e-12);
7534 assert!(terminal.velocity_magnitude < muzzle.velocity_magnitude);
7535 assert!(terminal.kinetic_energy < muzzle.kinetic_energy);
7536
7537 let terminal_sample = result
7538 .sampled_points
7539 .as_ref()
7540 .and_then(|samples| samples.last())
7541 .expect("terminal trajectory sample is missing");
7542 assert_eq!(
7543 terminal_sample.distance_m.to_bits(),
7544 target_range.to_bits(),
7545 "{name} sampling did not include max_range"
7546 );
7547 assert_eq!(
7548 terminal_sample.time_s.to_bits(),
7549 result.time_of_flight.to_bits()
7550 );
7551 assert_eq!(
7552 terminal_sample.velocity_mps.to_bits(),
7553 result.impact_velocity.to_bits()
7554 );
7555 assert!((terminal_sample.energy_j - result.impact_energy).abs() < 1e-12);
7556 }
7557 }
7558}
7559
7560#[cfg(test)]
7561mod precession_inertia_wiring_tests {
7562 use super::*;
7563
7564 #[test]
7565 fn solver_uses_projectile_specific_moments_of_inertia() {
7566 let mass_kg = 55.0 * crate::constants::GRAINS_TO_KG;
7567 let caliber_m = 0.224 * 0.0254;
7568 let length_m = 0.75 * 0.0254;
7569 let inputs = BallisticInputs {
7570 bullet_mass: mass_kg,
7571 bullet_diameter: caliber_m,
7572 bullet_length: length_m,
7573 muzzle_velocity: 800.0,
7574 twist_rate: 7.0,
7575 enable_precession_nutation: true,
7576 use_rk4: false,
7577 use_adaptive_rk45: false,
7578 ..BallisticInputs::default()
7579 };
7580 let mut solver = TrajectorySolver::new(
7581 inputs,
7582 WindConditions::default(),
7583 AtmosphericConditions::default(),
7584 );
7585 solver.set_max_range(0.1);
7586
7587 let (air_density, speed_of_sound, _, _) = solver.resolved_atmosphere();
7588 let velocity_mps = solver.inputs.muzzle_velocity;
7589 let velocity_fps = velocity_mps * 3.28084;
7590 let twist_rate_ft = solver.inputs.twist_rate / 12.0;
7591 let spin_rate_rad_s = (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI;
7592 let initial_state = AngularState {
7593 pitch_angle: 0.001,
7594 yaw_angle: 0.001,
7595 pitch_rate: 0.0,
7596 yaw_rate: 0.0,
7597 precession_angle: 0.0,
7598 nutation_phase: 0.0,
7599 };
7600 let params = PrecessionNutationParams {
7601 mass_kg,
7602 caliber_m,
7603 length_m,
7604 spin_rate_rad_s,
7605 spin_inertia: crate::spin_decay::calculate_moment_of_inertia(
7606 mass_kg, caliber_m, length_m, "ogive",
7607 ),
7608 transverse_inertia: crate::pitch_damping::calculate_transverse_moment_of_inertia(
7609 mass_kg, caliber_m, length_m, "ogive",
7610 ),
7611 velocity_mps,
7612 air_density_kg_m3: air_density,
7613 mach: velocity_mps / speed_of_sound,
7614 pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
7615 nutation_damping_factor: 0.05,
7616 };
7617 let expected = calculate_combined_angular_motion(
7618 ¶ms,
7619 &initial_state,
7620 0.0,
7621 solver.time_step,
7622 0.001,
7623 );
7624 let actual = solver
7625 .solve()
7626 .expect("one-step solve should succeed")
7627 .angular_state
7628 .expect("precession/nutation was enabled");
7629
7630 assert!(
7631 (actual.precession_angle - expected.precession_angle).abs() < 1e-15,
7632 "precession phase used the wrong inertia: actual={}, expected={}",
7633 actual.precession_angle,
7634 expected.precession_angle
7635 );
7636 assert!(
7637 (actual.nutation_phase - expected.nutation_phase).abs() < 1e-15,
7638 "nutation phase used the wrong inertia: actual={}, expected={}",
7639 actual.nutation_phase,
7640 expected.nutation_phase
7641 );
7642 }
7643}
7644
7645#[cfg(test)]
7646mod form_factor_drag_tests {
7647 use super::*;
7648
7649 fn acceleration_with_form_factor_flag(enabled: bool) -> Vector3<f64> {
7650 let inputs = BallisticInputs {
7651 bc_value: 0.462,
7652 bc_type: DragModel::G1,
7653 bullet_model: Some("168gr SMK Match".to_string()),
7654 use_form_factor: enabled,
7655 ..BallisticInputs::default()
7656 };
7657 let solver = TrajectorySolver::new(
7658 inputs,
7659 WindConditions::default(),
7660 AtmosphericConditions::default(),
7661 );
7662 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7663 solver.calculate_acceleration(
7664 &Vector3::zeros(),
7665 &Vector3::new(600.0, 0.0, 0.0),
7666 &Vector3::zeros(),
7667 (temp_c, pressure_hpa, density / 1.225),
7668 )
7669 }
7670
7671 #[test]
7672 fn measured_bc_drag_does_not_apply_name_based_form_factor_again() {
7673 let baseline = acceleration_with_form_factor_flag(false);
7674 let flagged = acceleration_with_form_factor_flag(true);
7675
7676 assert!(
7677 (flagged - baseline).norm() < 1e-12,
7678 "published BC already encodes form factor: baseline={baseline:?} flagged={flagged:?}"
7679 );
7680 }
7681}
7682
7683#[cfg(test)]
7684mod rk45_adaptivity_tests {
7685 use super::*;
7686
7687 #[test]
7688 fn cli_rk45_error_norm_scales_components_independently() {
7689 let position = Vector3::new(1.0e9, 0.0, 0.0);
7690 let velocity = Vector3::new(800.0, 0.0, 0.0);
7691 let fifth_position = position;
7692 let fifth_velocity = velocity;
7693 let fourth_position = position;
7694 let fourth_velocity = Vector3::new(800.0, 1.0e-3, 0.0);
7695
7696 let error = cli_rk45_error_norm(
7697 &position,
7698 &velocity,
7699 &fifth_position,
7700 &fifth_velocity,
7701 &fourth_position,
7702 &fourth_velocity,
7703 );
7704 let expected = 1.0e-3 / 6.0_f64.sqrt();
7705
7706 assert!(
7707 (error - expected).abs() <= 1e-15,
7708 "large downrange position masked a velocity-component error: {error}"
7709 );
7710 }
7711
7712 fn discontinuous_wind_solver() -> TrajectorySolver {
7713 let inputs = BallisticInputs::default();
7714 let mut solver = TrajectorySolver::new(
7715 inputs,
7716 WindConditions::default(),
7717 AtmosphericConditions::default(),
7718 );
7719 solver.set_wind_segments(vec![
7720 crate::wind::WindSegment::new(0.0, 90.0, 4.0),
7721 crate::wind::WindSegment::new(1_000.0, 90.0, 10_000.0),
7722 ]);
7723 solver
7724 }
7725
7726 #[test]
7727 fn rk45_retries_discontinuous_trial_before_advancing() {
7728 let solver = discontinuous_wind_solver();
7729 let position = Vector3::new(0.0, solver.inputs.muzzle_height, 0.0);
7730 let velocity = Vector3::new(solver.inputs.muzzle_velocity, 0.0, 0.0);
7731 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7732 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
7733 let dt = 0.01;
7734
7735 let rejected_trial = solver.rk45_step(
7736 &position,
7737 &velocity,
7738 dt,
7739 &Vector3::zeros(),
7740 RK45_TOLERANCE,
7741 resolved_atmo,
7742 );
7743 assert!(
7744 rejected_trial.error > RK45_TOLERANCE,
7745 "discontinuous full step must exceed tolerance, got {}",
7746 rejected_trial.error
7747 );
7748
7749 let accepted =
7750 solver.adaptive_rk45_step(&position, &velocity, dt, &Vector3::zeros(), resolved_atmo);
7751 assert!(accepted.used_dt < dt, "oversized trial was not retried");
7752 assert!(
7753 accepted.error <= RK45_TOLERANCE || accepted.used_dt <= RK45_MIN_DT,
7754 "accepted error {} exceeds tolerance at dt {}",
7755 accepted.error,
7756 accepted.used_dt
7757 );
7758
7759 let accepted_trial = solver.rk45_step(
7760 &position,
7761 &velocity,
7762 accepted.used_dt,
7763 &Vector3::zeros(),
7764 RK45_TOLERANCE,
7765 resolved_atmo,
7766 );
7767 assert_eq!(accepted.position, accepted_trial.position);
7768 assert_eq!(accepted.velocity, accepted_trial.velocity);
7769 assert!((RK45_MIN_DT..=RK45_MAX_DT).contains(&accepted.next_dt));
7770 }
7771}
7772
7773#[cfg(test)]
7774mod ground_termination_tests {
7775 use super::*;
7776 use crate::trajectory_observation::TrajectoryObservationFlag;
7777
7778 #[test]
7779 fn every_solver_reports_one_exact_early_ground_endpoint() {
7780 for (name, use_rk4, use_adaptive_rk45) in [
7781 ("Euler", false, false),
7782 ("RK4", true, false),
7783 ("RK45", true, true),
7784 ] {
7785 let inputs = BallisticInputs {
7786 muzzle_height: 1.0,
7787 muzzle_angle: -0.2,
7788 ground_threshold: 0.0,
7789 use_rk4,
7790 use_adaptive_rk45,
7791 ..BallisticInputs::default()
7792 };
7793 let mut solver = TrajectorySolver::new(
7794 inputs,
7795 WindConditions::default(),
7796 AtmosphericConditions::default(),
7797 );
7798 solver.set_max_range(1_000.0);
7799
7800 let result = solver.solve().expect("early-ground solve should succeed");
7801 let terminal = result.points.last().expect("terminal point is missing");
7802
7803 assert_eq!(result.termination, TrajectoryTermination::GroundThreshold);
7804 assert_eq!(terminal.position.y.to_bits(), 0.0_f64.to_bits());
7805 assert!(
7806 terminal.position.x < 1_000.0,
7807 "{name} incorrectly reached max range"
7808 );
7809 assert_eq!(result.max_range.to_bits(), terminal.position.x.to_bits());
7810 assert_eq!(
7811 result
7812 .points
7813 .iter()
7814 .filter(|point| point.position.y == 0.0)
7815 .count(),
7816 1,
7817 "{name} did not retain exactly one ground endpoint"
7818 );
7819
7820 let observations = result
7821 .sample_observations(1.0, 100)
7822 .expect("checked early-ground sampling should succeed");
7823 assert!(observations[..observations.len() - 1]
7824 .iter()
7825 .all(|observation| observation.distance_m < terminal.position.x));
7826 let terminal_observation = observations.last().expect("terminal observation");
7827 assert_eq!(
7828 terminal_observation.distance_m.to_bits(),
7829 terminal.position.x.to_bits()
7830 );
7831 assert!(terminal_observation
7832 .flags
7833 .contains(&TrajectoryObservationFlag::Terminal));
7834 assert!(terminal_observation
7835 .flags
7836 .contains(&TrajectoryObservationFlag::GroundThreshold));
7837 assert_eq!(
7838 observations
7839 .iter()
7840 .filter(|observation| observation
7841 .flags
7842 .contains(&TrajectoryObservationFlag::Terminal))
7843 .count(),
7844 1,
7845 "{name} repeated the terminal observation"
7846 );
7847 }
7848 }
7849
7850 #[test]
7855 fn rk4_and_rk45_descend_to_ground_threshold() {
7856 for adaptive in [false, true] {
7857 let inputs = BallisticInputs {
7858 muzzle_angle: 0.1, use_rk4: true,
7860 use_adaptive_rk45: adaptive,
7861 ..BallisticInputs::default()
7862 };
7863 assert_eq!(
7864 inputs.ground_threshold, -100.0,
7865 "default ground_threshold is -100 m"
7866 );
7867
7868 let mut solver = TrajectorySolver::new(
7869 inputs,
7870 WindConditions::default(),
7871 AtmosphericConditions::default(),
7872 );
7873 solver.set_max_range(1.0e7);
7875
7876 let result = solver.solve().expect("solve should succeed");
7877 let final_y = result
7878 .points
7879 .last()
7880 .expect("trajectory has points")
7881 .position
7882 .y;
7883 assert!(
7884 final_y < -1.0,
7885 "adaptive_rk45={adaptive}: final y = {final_y} m; a lofted shot should descend \
7886 past launch level toward the ground_threshold floor, not stop at y = 0"
7887 );
7888 }
7889 }
7890}
7891
7892#[cfg(test)]
7893mod magnus_stability_tests {
7894 use super::*;
7895
7896 #[test]
7897 fn yaw_of_repose_magnus_force_is_vertical_and_twist_signed() {
7898 let acceleration = |enable_magnus, is_twist_right| {
7899 let inputs = BallisticInputs {
7900 muzzle_velocity: 822.96,
7901 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
7902 bullet_diameter: 0.308 * 0.0254,
7903 bullet_length: 1.215 * 0.0254,
7904 twist_rate: 10.0,
7905 is_twist_right,
7906 enable_magnus,
7907 ..BallisticInputs::default()
7908 };
7909 let solver = TrajectorySolver::new(
7910 inputs,
7911 WindConditions::default(),
7912 AtmosphericConditions::default(),
7913 );
7914 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7915 solver.calculate_acceleration(
7916 &Vector3::zeros(),
7917 &Vector3::new(822.96, 0.0, 0.0),
7918 &Vector3::zeros(),
7919 (temp_c, pressure_hpa, density / 1.225),
7920 )
7921 };
7922
7923 let baseline = acceleration(false, true);
7924 let right_twist = acceleration(true, true) - baseline;
7925 let left_twist = acceleration(true, false) - baseline;
7926
7927 assert!(
7928 right_twist.y < 0.0,
7929 "right-hand Magnus must point down, got {right_twist:?}"
7930 );
7931 assert!(
7932 left_twist.y > 0.0,
7933 "left-hand Magnus must point up, got {left_twist:?}"
7934 );
7935 assert!((right_twist.y + left_twist.y).abs() < 1e-12);
7936 assert!(right_twist.x.abs() < 1e-12 && right_twist.z.abs() < 1e-12);
7937 assert!(left_twist.x.abs() < 1e-12 && left_twist.z.abs() < 1e-12);
7938 }
7939
7940 #[test]
7941 fn magnus_uses_velocity_corrected_muzzle_stability_gate() {
7942 let muzzle_velocity = 1_400.0 / 3.28084;
7943 let inputs = BallisticInputs {
7944 muzzle_velocity,
7945 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
7946 bullet_diameter: 0.308 * 0.0254,
7947 bullet_length: 1.215 * 0.0254,
7948 twist_rate: 15.0,
7949 enable_magnus: true,
7950 ..BallisticInputs::default()
7951 };
7952 let solver = TrajectorySolver::new(
7953 inputs.clone(),
7954 WindConditions::default(),
7955 AtmosphericConditions::default(),
7956 );
7957
7958 let bare_sg = crate::spin_drift::miller_stability(0.308, 168.0, 15.0, 1.215);
7959 let canonical_sg = solver.effective_spin_drift_sg();
7960 assert!(bare_sg > 1.0, "test requires bare Sg above the Magnus gate");
7961 assert!(
7962 canonical_sg < 1.0,
7963 "velocity-corrected Sg must be below the gate, got {canonical_sg}"
7964 );
7965
7966 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
7967 let acceleration = solver.calculate_acceleration(
7968 &Vector3::zeros(),
7969 &Vector3::new(muzzle_velocity, 0.0, 0.0),
7970 &Vector3::zeros(),
7971 (temp_c, pressure_hpa, density / 1.225),
7972 );
7973 let mut baseline_inputs = inputs;
7974 baseline_inputs.enable_magnus = false;
7975 let baseline_solver = TrajectorySolver::new(
7976 baseline_inputs,
7977 WindConditions::default(),
7978 AtmosphericConditions::default(),
7979 );
7980 let baseline = baseline_solver.calculate_acceleration(
7981 &Vector3::zeros(),
7982 &Vector3::new(muzzle_velocity, 0.0, 0.0),
7983 &Vector3::zeros(),
7984 (temp_c, pressure_hpa, density / 1.225),
7985 );
7986
7987 assert_eq!(
7988 acceleration, baseline,
7989 "canonical Sg below 1 must suppress every Magnus acceleration component"
7990 );
7991 }
7992
7993 #[test]
7994 fn magnus_force_grows_as_fixed_spin_projectile_slows() {
7995 let inputs = BallisticInputs {
7996 muzzle_velocity: 800.0,
7997 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
7998 bullet_diameter: 0.308 * 0.0254,
7999 bullet_length: 1.215 * 0.0254,
8000 twist_rate: 12.0,
8001 enable_magnus: true,
8002 ..BallisticInputs::default()
8003 };
8004
8005 let magnus_acceleration = |speed_mps| {
8006 let evaluate = |enable_magnus| {
8007 let mut run_inputs = inputs.clone();
8008 run_inputs.enable_magnus = enable_magnus;
8009 let solver = TrajectorySolver::new(
8010 run_inputs,
8011 WindConditions::default(),
8012 AtmosphericConditions::default(),
8013 );
8014 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
8015 solver
8016 .calculate_acceleration(
8017 &Vector3::zeros(),
8018 &Vector3::new(speed_mps, 0.0, 0.0),
8019 &Vector3::zeros(),
8020 (temp_c, pressure_hpa, density / 1.225),
8021 )
8022 .y
8023 };
8024 (evaluate(true) - evaluate(false)).abs()
8025 };
8026
8027 let fast = magnus_acceleration(200.0);
8028 let slow = magnus_acceleration(100.0);
8029 let ratio = slow / fast;
8030 let expected_ratio = 2.0_f64.powf(5.0 / 3.0);
8031
8032 assert!(fast > 0.0 && slow > 0.0, "fast={fast}, slow={slow}");
8033 assert!(
8034 (ratio - expected_ratio).abs() < 1e-3,
8035 "fixed-spin Magnus acceleration must grow downrange; slow/fast={ratio}, \
8036 expected={expected_ratio}"
8037 );
8038 }
8039}
8040
8041#[cfg(test)]
8042mod coriolis_direction_tests {
8043 use super::*;
8044 use std::f64::consts::FRAC_PI_2;
8045
8046 #[test]
8047 fn supersonic_crossing_flags_a_positive_range_sample() {
8048 use crate::trajectory_sampling::TrajectoryFlag;
8052
8053 for (solver_name, use_rk4, use_adaptive_rk45) in [
8054 ("Euler", false, false),
8055 ("RK4", true, false),
8056 ("RK45", true, true),
8057 ] {
8058 let inputs = BallisticInputs {
8059 muzzle_velocity: 850.0,
8060 bc_value: 0.2,
8061 bc_type: DragModel::G7,
8062 muzzle_angle: 0.03,
8063 enable_trajectory_sampling: true,
8064 sample_interval: 50.0,
8065 use_rk4,
8066 use_adaptive_rk45,
8067 ..BallisticInputs::default()
8068 };
8069 let mut solver = TrajectorySolver::new(
8070 inputs,
8071 WindConditions::default(),
8072 AtmosphericConditions::default(),
8073 );
8074 solver.set_max_range(2000.0);
8075 let samples = solver
8076 .solve()
8077 .expect("supersonic solve should succeed")
8078 .sampled_points
8079 .expect("sampling was enabled");
8080 let flagged_distances: Vec<_> = samples
8081 .iter()
8082 .filter(|sample| sample.flags.contains(&TrajectoryFlag::MachTransition))
8083 .map(|sample| sample.distance_m)
8084 .collect();
8085
8086 assert!(
8087 !flagged_distances.is_empty()
8088 && flagged_distances.iter().all(|distance| *distance > 0.0),
8089 "{solver_name} must flag genuine crossings only at positive range: {flagged_distances:?}"
8090 );
8091 }
8092 }
8093
8094 #[test]
8095 fn subsonic_launch_does_not_flag_a_muzzle_transition() {
8096 use crate::trajectory_sampling::TrajectoryFlag;
8097
8098 for (solver_name, use_rk4, use_adaptive_rk45) in [
8099 ("Euler", false, false),
8100 ("RK4", true, false),
8101 ("RK45", true, true),
8102 ] {
8103 let inputs = BallisticInputs {
8104 muzzle_velocity: 250.0,
8105 muzzle_angle: 0.02,
8106 enable_trajectory_sampling: true,
8107 sample_interval: 25.0,
8108 use_rk4,
8109 use_adaptive_rk45,
8110 ..BallisticInputs::default()
8111 };
8112 let mut solver = TrajectorySolver::new(
8113 inputs,
8114 WindConditions::default(),
8115 AtmosphericConditions::default(),
8116 );
8117 solver.set_max_range(300.0);
8118 let samples = solver
8119 .solve()
8120 .expect("subsonic solve should succeed")
8121 .sampled_points
8122 .expect("sampling was enabled");
8123
8124 assert!(
8125 samples
8126 .iter()
8127 .all(|sample| !sample.flags.contains(&TrajectoryFlag::MachTransition)),
8128 "{solver_name} marked a Mach transition for a launch already below Mach 1"
8129 );
8130 }
8131 }
8132
8133 #[test]
8134 fn mach_transition_tracker_requires_a_downward_crossing() {
8135 fn record(mach_values: &[f64]) -> Vec<f64> {
8136 let mut tracker = MachTransitionTracker::default();
8137 let mut distances = Vec::new();
8138 for (index, mach) in mach_values.iter().copied().enumerate() {
8139 tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
8140 }
8141 distances
8142 }
8143
8144 assert!(record(&[0.9, 0.8, 0.7]).is_empty());
8145 assert_eq!(record(&[1.1, 1.05, 0.99]), vec![20.0]);
8146 assert_eq!(record(&[1.2, 1.19, 1.0, 0.99]), vec![10.0, 30.0]);
8147 assert_eq!(record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]), vec![20.0, 30.0]);
8148 assert!(record(&[1.3, f64::NAN, 1.1]).is_empty());
8149 }
8150
8151 #[test]
8152 fn mach_transition_tracker_labels_0_9_without_touching_the_flat_vec() {
8153 fn record(mach_values: &[f64]) -> (Vec<f64>, MachTransitionTracker) {
8159 let mut tracker = MachTransitionTracker::default();
8160 let mut distances = Vec::new();
8161 for (index, mach) in mach_values.iter().copied().enumerate() {
8162 tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
8163 }
8164 (distances, tracker)
8165 }
8166
8167 let (distances, tracker) = record(&[0.9, 0.8, 0.7]);
8170 assert!(distances.is_empty()); assert_eq!(tracker.mach_1_2_distance_m, None);
8172 assert_eq!(tracker.mach_1_0_distance_m, None);
8173 assert_eq!(tracker.mach_0_9_distance_m, Some(10.0));
8174
8175 let (distances, tracker) = record(&[1.1, 1.05, 0.99]);
8177 assert_eq!(distances, vec![20.0]);
8178 assert_eq!(tracker.mach_1_2_distance_m, None);
8179 assert_eq!(tracker.mach_1_0_distance_m, Some(20.0));
8180 assert_eq!(tracker.mach_0_9_distance_m, None);
8181
8182 let (distances, tracker) = record(&[1.2, 1.19, 1.0, 0.99]);
8184 assert_eq!(distances, vec![10.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(10.0));
8186 assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
8187 assert_eq!(tracker.mach_0_9_distance_m, None);
8188
8189 let (distances, tracker) = record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]);
8192 assert_eq!(distances, vec![20.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(20.0));
8194 assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
8195 assert_eq!(tracker.mach_0_9_distance_m, Some(50.0));
8196 assert!(
8197 tracker.mach_1_2_distance_m < tracker.mach_1_0_distance_m
8198 && tracker.mach_1_0_distance_m < tracker.mach_0_9_distance_m,
8199 "labeled crossings must be strictly increasing downrange"
8200 );
8201
8202 let (distances, tracker) = record(&[1.3, f64::NAN, 1.1]);
8204 assert!(distances.is_empty());
8205 assert_eq!(tracker.mach_1_2_distance_m, None);
8206 assert_eq!(tracker.mach_1_0_distance_m, None);
8207 assert_eq!(tracker.mach_0_9_distance_m, None);
8208 }
8209
8210 #[test]
8211 fn humidity_percent_converts_and_clamps() {
8212 let mut i = BallisticInputs {
8214 humidity: 0.5,
8215 ..BallisticInputs::default()
8216 };
8217 assert!((i.humidity_percent() - 50.0).abs() < 1e-9, "0.5 -> 50%");
8218 i.humidity = 0.0;
8219 assert_eq!(i.humidity_percent(), 0.0);
8220 i.humidity = 1.0;
8221 assert_eq!(i.humidity_percent(), 100.0);
8222 i.humidity = 1.5; assert_eq!(i.humidity_percent(), 100.0);
8224 }
8225
8226 fn vertical_at(shot_azimuth: f64, range_m: f64) -> f64 {
8229 let inputs = BallisticInputs {
8230 muzzle_velocity: 800.0,
8231 bc_value: 0.5,
8232 bc_type: DragModel::G7,
8233 muzzle_angle: 0.02, enable_coriolis: true,
8235 latitude: Some(45.0),
8236 shot_azimuth,
8237 ground_threshold: f64::NEG_INFINITY, ..BallisticInputs::default()
8239 };
8240 let mut solver = TrajectorySolver::new(
8241 inputs,
8242 WindConditions::default(),
8243 AtmosphericConditions::default(),
8244 );
8245 solver.set_max_range(range_m + 50.0);
8246 let r = solver.solve().expect("solve");
8247 let pts = &r.points;
8248 for i in 1..pts.len() {
8249 if pts[i].position.x >= range_m {
8250 let p1 = &pts[i - 1];
8251 let p2 = &pts[i];
8252 let t = (range_m - p1.position.x) / (p2.position.x - p1.position.x);
8253 return p1.position.y + t * (p2.position.y - p1.position.y);
8254 }
8255 }
8256 panic!("range {range_m} not reached");
8257 }
8258
8259 #[test]
8264 fn eotvos_east_higher_than_west() {
8265 let range = 600.0;
8266 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!(
8270 east > west,
8271 "east ({east:.5}) must be higher than west ({west:.5}) at {range} m (Eötvös)"
8272 );
8273 assert!(
8274 east > north && north > west,
8275 "north ({north:.5}) must lie between east ({east:.5}) and west ({west:.5})"
8276 );
8277 assert!(
8278 (east - west) > 1e-3,
8279 "E-W vertical separation ({:.6} m) should be physically meaningful, not FP noise",
8280 east - west
8281 );
8282 }
8283
8284 #[test]
8292 fn labeled_mach_crossings_match_pinned_pre_change_flat_vec_across_solvers() {
8293 let cases = [
8295 (
8296 "Euler",
8297 false,
8298 false,
8299 670.9878683238721_f64,
8300 805.5274119916264_f64,
8301 ),
8302 (
8303 "RK4",
8304 true,
8305 false,
8306 671.7257336844475_f64,
8307 805.933409072171_f64,
8308 ),
8309 (
8310 "RK45",
8311 true,
8312 true,
8313 672.4905711917901_f64,
8314 806.5709746782849_f64,
8315 ),
8316 ];
8317
8318 for (solver_name, use_rk4, use_adaptive_rk45, expected_1_2, expected_1_0) in cases {
8319 let inputs = BallisticInputs {
8320 muzzle_velocity: 850.0,
8321 bc_value: 0.2,
8322 bc_type: DragModel::G7,
8323 muzzle_angle: 0.03,
8324 use_rk4,
8325 use_adaptive_rk45,
8326 ..BallisticInputs::default()
8327 };
8328 let mut solver = TrajectorySolver::new(
8329 inputs,
8330 WindConditions::default(),
8331 AtmosphericConditions::default(),
8332 );
8333 solver.set_max_range(2000.0);
8334 let result = solver.solve().expect("solve should succeed");
8335
8336 assert_eq!(
8337 result.mach_1_2_distance_m,
8338 Some(expected_1_2),
8339 "{solver_name}: mach_1_2_distance_m must match the pinned pre-change flat-Vec value"
8340 );
8341 assert_eq!(
8342 result.mach_1_0_distance_m,
8343 Some(expected_1_0),
8344 "{solver_name}: mach_1_0_distance_m must match the pinned pre-change flat-Vec value"
8345 );
8346
8347 let mach_1_2 = result.mach_1_2_distance_m.expect("crosses 1.2");
8348 let mach_1_0 = result.mach_1_0_distance_m.expect("crosses 1.0");
8349 let mach_0_9 = result
8350 .mach_0_9_distance_m
8351 .expect("this trajectory also goes past 0.9 within 2000 m");
8352 assert!(
8353 mach_1_2 < mach_1_0 && mach_1_0 < mach_0_9,
8354 "{solver_name}: labeled crossings must be strictly increasing downrange \
8355 (1.2={mach_1_2}, 1.0={mach_1_0}, 0.9={mach_0_9})"
8356 );
8357 }
8358 }
8359
8360 #[test]
8363 fn labeled_mach_crossings_are_none_for_a_fully_supersonic_trajectory() {
8364 for (solver_name, use_rk4, use_adaptive_rk45) in [
8365 ("Euler", false, false),
8366 ("RK4", true, false),
8367 ("RK45", true, true),
8368 ] {
8369 let inputs = BallisticInputs {
8370 muzzle_velocity: 850.0,
8371 bc_value: 0.2,
8372 bc_type: DragModel::G7,
8373 muzzle_angle: 0.03,
8374 use_rk4,
8375 use_adaptive_rk45,
8376 ..BallisticInputs::default()
8377 };
8378 let mut solver = TrajectorySolver::new(
8379 inputs,
8380 WindConditions::default(),
8381 AtmosphericConditions::default(),
8382 );
8383 solver.set_max_range(200.0);
8385 let result = solver.solve().expect("solve should succeed");
8386
8387 assert_eq!(
8388 result.mach_1_2_distance_m, None,
8389 "{solver_name}: must not report a 1.2 crossing that never happens"
8390 );
8391 assert_eq!(
8392 result.mach_1_0_distance_m, None,
8393 "{solver_name}: must not report a 1.0 crossing that never happens"
8394 );
8395 assert_eq!(
8396 result.mach_0_9_distance_m, None,
8397 "{solver_name}: must not report a 0.9 crossing that never happens"
8398 );
8399 }
8400 }
8401}
8402
8403#[cfg(test)]
8404mod cant_tests {
8405 use super::*;
8406
8407 fn base_inputs() -> BallisticInputs {
8408 BallisticInputs {
8409 muzzle_velocity: 800.0,
8410 bc_value: 0.5,
8411 bc_type: DragModel::G7,
8412 bullet_mass: 0.0109,
8413 bullet_diameter: 0.00782,
8414 bullet_length: 0.0309,
8415 sight_height: 0.05,
8416 twist_rate: 10.0,
8417 use_rk4: true,
8418 ..BallisticInputs::default()
8419 }
8420 }
8421
8422 fn solve_with(inputs: BallisticInputs, max_range: f64) -> TrajectoryResult {
8423 let mut s = TrajectorySolver::new(
8424 inputs,
8425 WindConditions::default(),
8426 AtmosphericConditions::default(),
8427 );
8428 s.set_max_range(max_range);
8429 s.solve().expect("solve")
8430 }
8431
8432 fn yz_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
8434 let pts = &result.points;
8435 for i in 1..pts.len() {
8436 if pts[i].position.x >= x {
8437 let (p1, p2) = (&pts[i - 1], &pts[i]);
8438 let dx = p2.position.x - p1.position.x;
8439 let t = if dx.abs() < 1e-12 {
8440 0.0
8441 } else {
8442 (x - p1.position.x) / dx
8443 };
8444 return (
8445 p1.position.y + t * (p2.position.y - p1.position.y),
8446 p1.position.z + t * (p2.position.z - p1.position.z),
8447 );
8448 }
8449 }
8450 panic!("trajectory never reached {x} m");
8451 }
8452
8453 #[test]
8454 fn cant_sign_clockwise_up_offset_goes_right_and_low() {
8455 let mut level = base_inputs();
8457 level.muzzle_angle = 0.003; let mut canted = level.clone();
8459 canted.cant_angle = 10f64.to_radians();
8460
8461 let (y0, z0) = yz_at(&solve_with(level, 400.0), 300.0);
8462 let (y1, z1) = yz_at(&solve_with(canted, 400.0), 300.0);
8463 assert!(
8464 z1 > z0 + 0.01,
8465 "clockwise cant must move POI right: z0={z0} z1={z1}"
8466 );
8467 assert!(
8468 y1 < y0 - 0.001,
8469 "clockwise cant must move POI low: y0={y0} y1={y1}"
8470 );
8471 }
8472
8473 #[test]
8474 fn pure_cant_shows_bore_offset_near_range() {
8475 let mut i = base_inputs();
8478 i.muzzle_angle = 0.0;
8479 i.cant_angle = 10f64.to_radians();
8480 let sh = i.sight_height;
8481 let r = solve_with(i, 60.0);
8482 let first = &r.points[1]; let expected = -sh * 10f64.to_radians().sin();
8484 assert!(
8485 (first.position.z - expected).abs() < 0.005,
8486 "near-muzzle lateral {} should be ~bore offset {expected}",
8487 first.position.z
8488 );
8489 }
8490
8491 #[test]
8492 fn zero_angle_is_independent_of_cant() {
8493 let a = base_inputs();
8494 let mut b = base_inputs();
8495 b.cant_angle = 15f64.to_radians();
8496 let za = calculate_zero_angle(a.clone(), 100.0, 0.0).expect("zero a");
8497 let zb = calculate_zero_angle(b.clone(), 100.0, 0.0).expect("zero b");
8498 assert_eq!(
8499 za.to_bits(),
8500 zb.to_bits(),
8501 "zeroing must ignore cant: {za} vs {zb}"
8502 );
8503 let _ = (a.cant_angle, b.cant_angle);
8505 }
8506
8507 #[test]
8508 fn nonfinite_cant_is_rejected() {
8509 let mut i = base_inputs();
8510 i.cant_angle = f64::NAN;
8511 let s = TrajectorySolver::new(
8512 i,
8513 WindConditions::default(),
8514 AtmosphericConditions::default(),
8515 );
8516 assert!(s.solve().is_err());
8517 }
8518
8519 #[test]
8520 fn incline_and_cant_compose_without_breaking() {
8521 let mut flat = base_inputs();
8523 flat.muzzle_angle = 0.003;
8524 flat.shooting_angle = 15f64.to_radians();
8525 let mut canted = flat.clone();
8526 canted.cant_angle = 10f64.to_radians();
8527 let (_, z_flat) = yz_at(&solve_with(flat, 400.0), 300.0);
8528 let (_, z_cant) = yz_at(&solve_with(canted, 400.0), 300.0);
8529 assert!(
8530 z_cant > z_flat,
8531 "cant must still deflect right on an incline"
8532 );
8533 }
8534}
8535
8536#[cfg(test)]
8537mod vertical_wind_tests {
8538 use super::*;
8539
8540 fn base_inputs() -> BallisticInputs {
8541 BallisticInputs {
8542 muzzle_velocity: 800.0,
8543 bc_value: 0.5,
8544 bc_type: DragModel::G7,
8545 bullet_mass: 0.0109,
8546 bullet_diameter: 0.00782,
8547 bullet_length: 0.0309,
8548 sight_height: 0.05,
8549 twist_rate: 10.0,
8550 use_rk4: true,
8551 ..BallisticInputs::default()
8552 }
8553 }
8554
8555 fn y_at(result: &TrajectoryResult, x: f64) -> f64 {
8557 let pts = &result.points;
8558 for i in 1..pts.len() {
8559 if pts[i].position.x >= x {
8560 let (p1, p2) = (&pts[i - 1], &pts[i]);
8561 let dx = p2.position.x - p1.position.x;
8562 let t = if dx.abs() < 1e-12 {
8563 0.0
8564 } else {
8565 (x - p1.position.x) / dx
8566 };
8567 return p1.position.y + t * (p2.position.y - p1.position.y);
8568 }
8569 }
8570 panic!("trajectory never reached {x} m");
8571 }
8572
8573 fn solve_with(
8574 inputs: BallisticInputs,
8575 wind: WindConditions,
8576 max_range: f64,
8577 ) -> TrajectoryResult {
8578 let mut s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
8579 s.set_max_range(max_range);
8580 s.solve().expect("solve")
8581 }
8582
8583 #[test]
8584 fn updraft_raises_poi_downrange() {
8585 let calm_inputs = base_inputs();
8588 let calm_wind = WindConditions::default();
8589 let updraft = WindConditions {
8590 vertical_speed: 5.0,
8591 ..Default::default()
8592 };
8593
8594 let calm = solve_with(calm_inputs.clone(), calm_wind, 500.0);
8595 let updraft_result = solve_with(calm_inputs, updraft, 500.0);
8596
8597 let y_calm = y_at(&calm, 400.0);
8598 let y_updraft = y_at(&updraft_result, 400.0);
8599 assert!(
8600 y_updraft > y_calm,
8601 "5 m/s updraft must raise POI at 400m: calm={y_calm}, updraft={y_updraft}"
8602 );
8603 }
8604
8605 #[test]
8606 fn zero_vertical_is_default_and_finite_required() {
8607 assert_eq!(WindConditions::default().vertical_speed, 0.0);
8608
8609 let inputs = base_inputs();
8610 let wind = WindConditions {
8611 vertical_speed: f64::NAN,
8612 ..Default::default()
8613 };
8614 let s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
8615 assert!(
8616 s.solve().is_err(),
8617 "NaN wind.vertical_speed must be rejected by validate_for_solve"
8618 );
8619 }
8620}
8621
8622#[cfg(test)]
8624mod bc_reference_standard_tests {
8625 use super::*;
8626
8627 fn base_inputs() -> BallisticInputs {
8628 BallisticInputs {
8629 muzzle_velocity: 800.0,
8630 bc_value: 0.5,
8631 bc_type: DragModel::G7,
8632 bullet_mass: 0.0109,
8633 bullet_diameter: 0.00782,
8634 bullet_length: 0.0309,
8635 sight_height: 0.05,
8636 twist_rate: 10.0,
8637 use_rk4: true,
8638 ..BallisticInputs::default()
8639 }
8640 }
8641
8642 fn y_and_speed_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
8644 let pts = &result.points;
8645 for i in 1..pts.len() {
8646 if pts[i].position.x >= x {
8647 let (p1, p2) = (&pts[i - 1], &pts[i]);
8648 let dx = p2.position.x - p1.position.x;
8649 let t = if dx.abs() < 1e-12 {
8650 0.0
8651 } else {
8652 (x - p1.position.x) / dx
8653 };
8654 return (
8655 p1.position.y + t * (p2.position.y - p1.position.y),
8656 p1.velocity_magnitude + t * (p2.velocity_magnitude - p1.velocity_magnitude),
8657 );
8658 }
8659 }
8660 panic!("trajectory never reached {x} m");
8661 }
8662
8663 #[test]
8666 fn asm_to_icao_ratio_matches_documented_value() {
8667 assert!(
8668 (crate::constants::ASM_TO_ICAO_BC - 0.98237).abs() < 1e-5,
8669 "ASM_TO_ICAO_BC = {} must equal 0.98237 to 5 decimal places",
8670 crate::constants::ASM_TO_ICAO_BC
8671 );
8672 assert_eq!(
8677 crate::constants::ASM_TO_ICAO_BC,
8678 crate::constants::ASM_DENSITY_LB_FT3 / crate::constants::ICAO_DENSITY_LB_FT3
8679 );
8680 }
8681
8682 #[test]
8685 fn default_bc_reference_standard_is_icao() {
8686 assert_eq!(
8687 BallisticInputs::default().bc_reference_standard,
8688 BcReferenceStandard::Icao
8689 );
8690 }
8691
8692 #[test]
8697 fn icao_reference_leaves_bc_value_bit_identical() {
8698 let raw_bc: f64 = 0.4372911; let inputs = BallisticInputs {
8700 bc_value: raw_bc,
8701 bc_reference_standard: BcReferenceStandard::Icao,
8702 ..base_inputs()
8703 };
8704 let solver = TrajectorySolver::new(
8705 inputs,
8706 WindConditions::default(),
8707 AtmosphericConditions::default(),
8708 );
8709 assert_eq!(solver.inputs.bc_value.to_bits(), raw_bc.to_bits());
8710 }
8711
8712 #[test]
8713 fn default_inputs_solve_is_unaffected_by_the_new_field_existing() {
8714 let a = TrajectorySolver::new(
8718 base_inputs(),
8719 WindConditions::default(),
8720 AtmosphericConditions::default(),
8721 )
8722 .solve()
8723 .expect("solve a");
8724 let b = TrajectorySolver::new(
8725 BallisticInputs { ..base_inputs() },
8726 WindConditions::default(),
8727 AtmosphericConditions::default(),
8728 )
8729 .solve()
8730 .expect("solve b");
8731 assert_eq!(a.impact_velocity.to_bits(), b.impact_velocity.to_bits());
8732 assert_eq!(a.max_range.to_bits(), b.max_range.to_bits());
8733 }
8734
8735 #[test]
8738 fn army_standard_metro_scales_bc_value_by_exactly_the_derived_ratio() {
8739 let raw_bc = 0.5;
8740 let inputs = BallisticInputs {
8741 bc_value: raw_bc,
8742 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8743 ..base_inputs()
8744 };
8745 let solver = TrajectorySolver::new(
8746 inputs,
8747 WindConditions::default(),
8748 AtmosphericConditions::default(),
8749 );
8750 assert_eq!(
8751 solver.inputs.bc_value,
8752 raw_bc * crate::constants::ASM_TO_ICAO_BC
8753 );
8754 }
8755
8756 #[test]
8757 fn army_standard_metro_scales_mach_keyed_bc_segments() {
8758 let inputs = BallisticInputs {
8759 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8760 bc_segments: Some(vec![(0.5, 0.40), (1.5, 0.30)]),
8761 ..base_inputs()
8762 };
8763 let solver = TrajectorySolver::new(
8764 inputs,
8765 WindConditions::default(),
8766 AtmosphericConditions::default(),
8767 );
8768 let segments = solver.inputs.bc_segments.as_ref().expect("segments");
8769 assert_eq!(segments[0], (0.5, 0.40 * crate::constants::ASM_TO_ICAO_BC));
8770 assert_eq!(segments[1], (1.5, 0.30 * crate::constants::ASM_TO_ICAO_BC));
8771 }
8772
8773 #[test]
8774 fn army_standard_metro_scales_velocity_keyed_bc_segments_data() {
8775 let inputs = BallisticInputs {
8776 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8777 bc_segments_data: Some(vec![
8778 crate::BCSegmentData {
8779 velocity_min: 0.0,
8780 velocity_max: 500.0,
8781 bc_value: 0.40,
8782 },
8783 crate::BCSegmentData {
8784 velocity_min: 500.0,
8785 velocity_max: 900.0,
8786 bc_value: 0.45,
8787 },
8788 ]),
8789 ..base_inputs()
8790 };
8791 let solver = TrajectorySolver::new(
8792 inputs,
8793 WindConditions::default(),
8794 AtmosphericConditions::default(),
8795 );
8796 let segments = solver.inputs.bc_segments_data.as_ref().expect("segments");
8797 assert_eq!(
8798 segments[0].bc_value,
8799 0.40 * crate::constants::ASM_TO_ICAO_BC
8800 );
8801 assert_eq!(
8802 segments[1].bc_value,
8803 0.45 * crate::constants::ASM_TO_ICAO_BC
8804 );
8805 assert_eq!(segments[0].velocity_min, 0.0);
8807 assert_eq!(segments[1].velocity_max, 900.0);
8808 }
8809
8810 #[test]
8815 fn army_standard_metro_moves_impact_in_the_more_drag_direction() {
8816 let solve_at = |standard: BcReferenceStandard| {
8817 let inputs = BallisticInputs {
8818 bc_value: 0.475,
8819 bc_reference_standard: standard,
8820 ..base_inputs()
8821 };
8822 let mut solver = TrajectorySolver::new(
8823 inputs,
8824 WindConditions::default(),
8825 AtmosphericConditions::default(),
8826 );
8827 solver.set_max_range(500.0);
8828 solver.solve().expect("solve")
8829 };
8830
8831 let icao = solve_at(BcReferenceStandard::Icao);
8832 let asm = solve_at(BcReferenceStandard::ArmyStandardMetro);
8833
8834 let (y_icao, v_icao) = y_and_speed_at(&icao, 400.0);
8835 let (y_asm, v_asm) = y_and_speed_at(&asm, 400.0);
8836
8837 assert!(
8838 y_asm < y_icao,
8839 "ArmyStandardMetro must drop MORE (lower y) at 400m than Icao for the same raw \
8840 bc_value: icao_y={y_icao}, asm_y={y_asm}"
8841 );
8842 assert!(
8843 v_asm < v_icao,
8844 "ArmyStandardMetro must retain LESS velocity at 400m than Icao for the same raw \
8845 bc_value: icao_v={v_icao}, asm_v={v_asm}"
8846 );
8847 }
8848
8849 #[test]
8856 fn monte_carlo_inherits_the_normalized_bc_reference() {
8857 let base_inputs_asm = BallisticInputs {
8858 bc_value: 0.475,
8859 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
8860 ..base_inputs()
8861 };
8862 let wind = WindConditions::default();
8863
8864 let mut direct_solver = TrajectorySolver::new(
8870 base_inputs_asm.clone(),
8871 wind.clone(),
8872 AtmosphericConditions::default(),
8873 );
8874 direct_solver.set_max_range(base_inputs_asm.target_distance.max(1000.0) * 2.0);
8875 let direct = direct_solver.solve().expect("direct solve");
8876
8877 let mc_params = MonteCarloParams {
8878 num_simulations: 1,
8879 velocity_std_dev: 0.0,
8880 angle_std_dev: 0.0,
8881 bc_std_dev: 0.0,
8882 wind_speed_std_dev: 0.0,
8883 target_distance: None,
8884 base_wind_speed: 0.0,
8885 base_wind_direction: 0.0,
8886 azimuth_std_dev: 0.0,
8887 };
8888 let mc = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
8889 base_inputs_asm,
8890 wind,
8891 mc_params,
8892 0.0,
8893 42,
8894 )
8895 .expect("monte carlo");
8896
8897 assert_eq!(mc.ranges.len(), 1);
8898 assert_eq!(
8899 mc.ranges[0].to_bits(),
8900 direct.max_range.to_bits(),
8901 "a zero-dispersion single MC sample must match a plain solve of the same \
8902 ASM-referenced inputs bit-for-bit"
8903 );
8904 assert_eq!(
8905 mc.impact_velocities[0].to_bits(),
8906 direct.impact_velocity.to_bits()
8907 );
8908 }
8909
8910 #[test]
8918 fn estimate_bc_fit_recovers_an_icao_referenced_bc() {
8919 let known_bc = 0.475;
8920 let velocity = 800.0;
8921 let mass = 0.0109;
8922 let diameter = 0.00782;
8923 let atmosphere = AtmosphericConditions::default();
8924
8925 let synth_inputs = BallisticInputs {
8926 muzzle_velocity: velocity,
8927 bc_value: known_bc,
8928 bc_type: DragModel::G7,
8929 bullet_mass: mass,
8930 bullet_diameter: diameter,
8931 bullet_length: 0.0309,
8932 sight_height: 0.05,
8933 twist_rate: 10.0,
8934 use_rk4: true,
8935 bc_reference_standard: BcReferenceStandard::Icao,
8936 ..BallisticInputs::default()
8937 };
8938 let mut solver =
8939 TrajectorySolver::new(synth_inputs, WindConditions::default(), atmosphere.clone());
8940 solver.set_max_range(500.0);
8941 let trajectory = solver.solve().expect("synthetic solve");
8942
8943 let points: Vec<(f64, f64)> = [100.0, 200.0, 300.0, 400.0]
8944 .iter()
8945 .map(|&d| {
8946 let (y, _) = {
8947 let pts = &trajectory.points;
8948 let mut found = None;
8949 for i in 1..pts.len() {
8950 if pts[i].position.x >= d {
8951 let (p1, p2) = (&pts[i - 1], &pts[i]);
8952 let dx = p2.position.x - p1.position.x;
8953 let t = if dx.abs() < 1e-12 {
8954 0.0
8955 } else {
8956 (d - p1.position.x) / dx
8957 };
8958 found =
8959 Some((p1.position.y + t * (p2.position.y - p1.position.y), 0.0));
8960 break;
8961 }
8962 }
8963 found.expect("trajectory reached observation distance")
8964 };
8965 (d, -y) })
8967 .collect();
8968
8969 let estimate = estimate_bc_fit(
8970 velocity,
8971 mass,
8972 diameter,
8973 &points,
8974 DragModel::G7,
8975 BcFitMode::Drop,
8976 atmosphere,
8977 None,
8978 0.05,
8979 )
8980 .expect("fit should converge");
8981
8982 assert!(
8983 (estimate.bc - known_bc).abs() < 0.02,
8984 "fit should recover the known ICAO-referenced bc={known_bc}, got {}",
8985 estimate.bc
8986 );
8987 }
8988
8989 #[test]
8992 fn custom_drag_table_makes_bc_reference_standard_numerically_inert() {
8993 let table =
8994 crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0, 3.0], vec![0.3, 0.4, 0.3, 0.2])
8995 .expect("valid table");
8996
8997 let solve_with = |standard: BcReferenceStandard| {
8998 let inputs = BallisticInputs {
8999 bc_value: 0.5, bc_reference_standard: standard,
9001 custom_drag_table: Some(table.clone()),
9002 ..base_inputs()
9003 };
9004 let mut solver = TrajectorySolver::new(
9005 inputs,
9006 WindConditions::default(),
9007 AtmosphericConditions::default(),
9008 );
9009 solver.set_max_range(500.0);
9010 solver.solve().expect("solve")
9011 };
9012
9013 let icao = solve_with(BcReferenceStandard::Icao);
9014 let asm = solve_with(BcReferenceStandard::ArmyStandardMetro);
9015
9016 assert_eq!(
9017 icao.impact_velocity.to_bits(),
9018 asm.impact_velocity.to_bits(),
9019 "a custom drag table must make bc_reference_standard fully inert"
9020 );
9021 assert_eq!(icao.max_range.to_bits(), asm.max_range.to_bits());
9022 }
9023
9024 #[test]
9025 fn custom_drag_table_inert_warning_fires_only_for_army_standard_metro_with_a_table() {
9026 let table = crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0], vec![0.3, 0.4, 0.3])
9027 .expect("valid table");
9028
9029 let no_table_icao = base_inputs();
9031 assert!(no_table_icao
9032 .bc_reference_standard_inert_warning()
9033 .is_none());
9034 let no_table_asm = BallisticInputs {
9035 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
9036 ..base_inputs()
9037 };
9038 assert!(no_table_asm.bc_reference_standard_inert_warning().is_none());
9039
9040 let table_icao = BallisticInputs {
9042 custom_drag_table: Some(table.clone()),
9043 ..base_inputs()
9044 };
9045 assert!(table_icao.bc_reference_standard_inert_warning().is_none());
9046
9047 let table_asm = BallisticInputs {
9049 custom_drag_table: Some(table),
9050 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
9051 ..base_inputs()
9052 };
9053 let warning = table_asm
9054 .bc_reference_standard_inert_warning()
9055 .expect("must warn");
9056 assert!(warning.contains("--bc-reference"));
9057 assert!(warning.contains("--drag-table"));
9058 }
9059}
9060
9061#[cfg(test)]
9067mod effective_drag_coefficient_tests {
9068 use super::*;
9069
9070 fn inputs_175gr_g7() -> BallisticInputs {
9071 let mut inputs = BallisticInputs {
9072 bc_value: 0.243,
9073 bc_type: DragModel::G7,
9074 muzzle_velocity: 823.0,
9075 ..Default::default()
9076 };
9077 inputs.bullet_mass = 175.0 * crate::constants::GRAINS_TO_KG;
9081 inputs.bullet_diameter = 0.308 * 0.0254;
9082 inputs.weight_grains = 175.0;
9083 inputs.caliber_inches = 0.308;
9084 inputs
9085 }
9086
9087 fn solver(inputs: BallisticInputs) -> TrajectorySolver {
9088 TrajectorySolver::new(
9089 inputs,
9090 WindConditions::default(),
9091 AtmosphericConditions::default(),
9092 )
9093 }
9094
9095 #[test]
9099 fn reports_the_projectiles_own_cd_not_the_reference_tables() {
9100 let inputs = inputs_175gr_g7();
9101 let sd = inputs.sectional_density_lb_in2().expect("SD");
9102 let solver = solver(inputs);
9103
9104 let sos = 340.0;
9105 let velocity = 800.0;
9106 let mach = velocity / sos;
9107
9108 let reference = crate::drag::get_drag_coefficient(mach, &DragModel::G7);
9109 let reported = solver
9110 .effective_drag_coefficient(velocity, sos)
9111 .expect("mass and diameter are set");
9112
9113 let expected = reference * sd / 0.243;
9114 assert!(
9115 (reported - expected).abs() < 1e-12,
9116 "reported {reported} != Cd_ref * SD / BC {expected}"
9117 );
9118 assert!(
9121 (reported - reference).abs() > 1e-6,
9122 "form factor collapsed to 1; this fixture no longer distinguishes the two values"
9123 );
9124 }
9125
9126 #[test]
9129 fn a_custom_drag_table_passes_through_unscaled() {
9130 let mut inputs = inputs_175gr_g7();
9131 inputs.custom_drag_table = Some(crate::drag::DragTable::new(
9132 vec![0.5, 3.0],
9133 vec![0.15, 0.40],
9134 ));
9135 let solver = solver(inputs);
9136
9137 let sos = 340.0;
9138 let velocity = 0.9 * sos;
9139 let table_value = solver
9140 .inputs
9141 .custom_drag_table
9142 .as_ref()
9143 .expect("table")
9144 .interpolate(0.9);
9145
9146 let reported = solver
9147 .effective_drag_coefficient(velocity, sos)
9148 .expect("mass and diameter are set");
9149 assert!(
9150 (reported - table_value).abs() < 1e-12,
9151 "custom table Cd {table_value} was rescaled to {reported}"
9152 );
9153 }
9154
9155 #[test]
9158 fn a_velocity_segmented_bc_steps_the_reported_cd() {
9159 let mut inputs = inputs_175gr_g7();
9160 inputs.use_bc_segments = true;
9161 inputs.bc_segments_data = Some(vec![
9162 crate::BCSegmentData {
9163 velocity_min: 2400.0,
9164 velocity_max: 4000.0,
9165 bc_value: 0.243,
9166 },
9167 crate::BCSegmentData {
9168 velocity_min: 0.0,
9169 velocity_max: 2400.0,
9170 bc_value: 0.200,
9171 },
9172 ]);
9173 let solver = solver(inputs);
9174
9175 let sos = 340.0;
9176 let above = solver
9178 .effective_drag_coefficient(2500.0 / 3.28084, sos)
9179 .expect("cd");
9180 let below = solver
9181 .effective_drag_coefficient(2300.0 / 3.28084, sos)
9182 .expect("cd");
9183
9184 assert!(
9186 below > above,
9187 "expected the 0.200 band to report a higher Cd than the 0.243 band; got {below} vs {above}"
9188 );
9189 }
9190
9191 #[test]
9194 fn is_absent_when_sectional_density_is_unknown() {
9195 let mut inputs = inputs_175gr_g7();
9196 inputs.weight_grains = 0.0;
9197 inputs.bullet_mass = 0.0;
9198 let solver = solver(inputs);
9199 assert!(solver.effective_drag_coefficient(800.0, 340.0).is_none());
9200 }
9201
9202 #[test]
9207 fn the_json_emit_rule_is_flag_gated_and_absent_when_cd_is_unknown() {
9208 let mut point = TrajectoryPoint {
9209 time: 0.0,
9210 position: nalgebra::Vector3::new(0.0, 0.0, 0.0),
9211 velocity_magnitude: 800.0,
9212 kinetic_energy: 3000.0,
9213 drag_coefficient: Some(0.31),
9214 };
9215 assert_eq!(point.drag_coefficient_json_value(true), Some(0.31));
9216 assert_eq!(
9217 point.drag_coefficient_json_value(false),
9218 None,
9219 "without the flag the key must not exist, so default JSON stays byte-identical"
9220 );
9221 point.drag_coefficient = None;
9222 assert_eq!(
9223 point.drag_coefficient_json_value(true),
9224 None,
9225 "unknown sectional density must yield an ABSENT key, not null"
9226 );
9227 }
9228
9229 #[test]
9231 fn every_point_of_a_solved_trajectory_carries_the_value() {
9232 let mut solver = solver(inputs_175gr_g7());
9233 solver.set_max_range(300.0);
9234 let result = solver.solve().expect("solve");
9235 assert!(!result.points.is_empty());
9236 assert!(
9237 result.points.iter().all(|p| p.drag_coefficient.is_some()),
9238 "the post-integration pass missed at least one point"
9239 );
9240 }
9241}