1use crate::cluster_bc::ClusterBCDegradation;
3use crate::pitch_damping::{calculate_pitch_damping_coefficient, PitchDampingCoefficients};
4use crate::precession_nutation::{
5 calculate_combined_angular_motion, projectile_moments_of_inertia, AngularState,
6 PrecessionNutationParams,
7};
8use crate::trajectory_sampling::{
9 projected_sample_count, sample_trajectory, TrajectoryData, TrajectoryOutputs,
10 TrajectorySample,
11};
12use crate::trajectory_observation::TrajectoryTermination;
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 use_bc_segments: bool,
248 pub bc_segments: Option<Vec<(f64, f64)>>, pub bc_segments_data: Option<Vec<crate::BCSegmentData>>, pub use_enhanced_spin_drift: bool,
251 pub use_form_factor: bool,
254 pub enable_wind_shear: bool,
255 pub wind_shear_model: String,
256 pub enable_trajectory_sampling: bool,
257 pub sample_interval: f64, pub drops_reference: DropsReference,
267 pub enable_pitch_damping: bool,
268 pub enable_precession_nutation: bool,
269 pub enable_aerodynamic_jump: bool,
272 pub use_cluster_bc: bool, pub custom_drag_table: Option<crate::drag::DragTable>,
276 pub cd_scale: f64,
284
285 pub bc_type_str: Option<String>,
287}
288
289impl BallisticInputs {
290 pub fn humidity_percent(&self) -> f64 {
295 (self.humidity * 100.0).clamp(0.0, 100.0)
296 }
297
298 pub fn windage_zero_bias_rad(&self, zero_distance_m: f64) -> f64 {
313 if zero_distance_m > 0.0 {
314 (self.zero_poi_horizontal_m + self.sight_offset_lateral_m) / zero_distance_m
315 } else {
316 0.0
317 }
318 }
319
320 pub fn sectional_density_lb_in2(&self) -> Option<f64> {
326 let weight_gr = if self.weight_grains > 0.0 {
327 self.weight_grains
328 } else {
329 self.bullet_mass / crate::constants::GRAINS_TO_KG };
331 let diameter_in = if self.caliber_inches > 0.0 {
332 self.caliber_inches
333 } else {
334 self.bullet_diameter / 0.0254 };
336 if weight_gr > 0.0 && diameter_in > 0.0 {
337 Some(weight_gr / 7000.0 / (diameter_in * diameter_in))
338 } else {
339 None
340 }
341 }
342
343 pub fn custom_drag_denominator(&self, fallback_bc: f64) -> f64 {
355 match self.sectional_density_lb_in2() {
356 Some(sd) => sd,
357 None => {
358 static WARN_ONCE: std::sync::Once = std::sync::Once::new();
359 WARN_ONCE.call_once(|| {
360 eprintln!(
361 "Warning: custom drag table active but bullet mass/diameter are \
362 unavailable; falling back to bc_value for the retardation denominator"
363 );
364 });
365 fallback_bc
366 }
367 }
368 }
369
370 pub fn bc_reference_standard_inert_warning(&self) -> Option<&'static str> {
381 if self.custom_drag_table.is_some()
382 && matches!(self.bc_reference_standard, BcReferenceStandard::ArmyStandardMetro)
383 {
384 Some(BC_REFERENCE_STANDARD_INERT_WARNING)
385 } else {
386 None
387 }
388 }
389
390 pub fn normalize_for_solve(&mut self) {
409 if matches!(
420 self.bc_reference_standard,
421 BcReferenceStandard::ArmyStandardMetro
422 ) {
423 self.bc_value *= crate::constants::ASM_TO_ICAO_BC;
424 if let Some(segments) = self.bc_segments.as_mut() {
425 for (_mach, bc) in segments.iter_mut() {
426 *bc *= crate::constants::ASM_TO_ICAO_BC;
427 }
428 }
429 if let Some(segments) = self.bc_segments_data.as_mut() {
430 for segment in segments.iter_mut() {
431 segment.bc_value *= crate::constants::ASM_TO_ICAO_BC;
432 }
433 }
434 self.bc_reference_standard = BcReferenceStandard::Icao;
437 }
438
439 self.caliber_inches = self.bullet_diameter / 0.0254;
444 self.weight_grains = self.bullet_mass / crate::constants::GRAINS_TO_KG;
445
446 self.muzzle_velocity = resolve_powder_adjusted_velocity(
457 self.muzzle_velocity,
458 self.temperature,
459 self.use_powder_sensitivity,
460 self.powder_temp_sensitivity,
461 self.powder_temp,
462 self.powder_temp_curve.as_deref(),
463 self.powder_curve_temp_c,
464 );
465 }
466}
467
468impl Default for BallisticInputs {
469 fn default() -> Self {
470 let mass_kg = 0.01;
471 let diameter_m = 0.00762;
472 let bc = 0.5;
473 let muzzle_angle_rad = 0.0;
474 let bc_type = DragModel::G1;
475
476 Self {
477 bc_value: bc,
479 bc_type,
480 bc_reference_standard: BcReferenceStandard::Icao,
481 bullet_mass: mass_kg,
482 muzzle_velocity: 800.0,
483 bullet_diameter: diameter_m,
484 bullet_length: crate::stability::estimate_bullet_length_m(diameter_m, mass_kg),
488
489 muzzle_angle: muzzle_angle_rad,
491 target_distance: 100.0,
492 azimuth_angle: 0.0,
493 shot_azimuth: 0.0,
494 shooting_angle: 0.0,
495 cant_angle: 0.0,
496 sight_height: 0.05,
497 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,
502 ground_threshold: -100.0, altitude: 0.0,
506 temperature: 15.0,
507 pressure: 1013.25, humidity: 0.5, latitude: None,
510
511 wind_speed: 0.0,
513 wind_angle: 0.0,
514
515 twist_rate: 12.0, is_twist_right: true,
518 caliber_inches: diameter_m / 0.0254, weight_grains: mass_kg / crate::constants::GRAINS_TO_KG, manufacturer: None,
521 bullet_model: None,
522 bullet_id: None,
523 bullet_cluster: None,
524
525 use_rk4: true, use_adaptive_rk45: true, enable_advanced_effects: false,
531 enable_magnus: false,
532 enable_coriolis: false,
533 use_powder_sensitivity: false,
534 powder_temp_sensitivity: 0.0,
535 powder_temp: 15.0,
536 powder_temp_curve: None,
537 powder_curve_temp_c: None,
538 tipoff_yaw: 0.0,
539 tipoff_decay_distance: 50.0,
540 use_bc_segments: false,
541 bc_segments: None,
542 bc_segments_data: None,
543 use_enhanced_spin_drift: false,
544 use_form_factor: false,
545 enable_wind_shear: false,
546 wind_shear_model: "none".to_string(),
547 enable_trajectory_sampling: false,
548 sample_interval: 10.0, drops_reference: DropsReference::Los, enable_pitch_damping: false,
551 enable_precession_nutation: false,
552 enable_aerodynamic_jump: false,
553 use_cluster_bc: false, custom_drag_table: None,
557 cd_scale: 1.0,
558
559 bc_type_str: None,
561 }
562 }
563}
564
565pub fn interpolate_powder_temp_curve(curve: &[(f64, f64)], temp_c: f64) -> f64 {
571 debug_assert!(!curve.is_empty());
572 if curve.is_empty() {
573 return 0.0;
574 }
575 let mut sorted;
578 let pts: &[(f64, f64)] = if curve.windows(2).all(|w| w[0].0 <= w[1].0) {
579 curve
580 } else {
581 sorted = curve.to_vec();
582 sorted.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
583 &sorted
584 };
585 let n = pts.len();
586 if temp_c <= pts[0].0 {
587 return pts[0].1; }
589 if temp_c >= pts[n - 1].0 {
590 return pts[n - 1].1; }
592 for i in 1..n {
593 let (t0, v0) = pts[i - 1];
594 let (t1, v1) = pts[i];
595 if temp_c <= t1 {
596 let span = t1 - t0;
597 if span.abs() < f64::EPSILON {
598 return v1; }
600 let f = (temp_c - t0) / span;
601 return v0 + f * (v1 - v0);
602 }
603 }
604 pts[n - 1].1
605}
606
607pub fn parse_powder_sweep(s: &str) -> Result<Vec<f64>, String> {
613 const MAX_SWEEP_ROWS: usize = 500;
614 let parts: Vec<&str> = s.split(':').collect();
615 if parts.len() != 3 {
616 return Err(format!(
617 "Invalid --sweep '{}': expected START:END:STEP (e.g. \"20:110:10\")",
618 s
619 ));
620 }
621 let parse = |p: &str, name: &str| -> Result<f64, String> {
622 p.trim()
623 .parse::<f64>()
624 .map_err(|_| format!("Invalid --sweep {}: '{}' is not a number", name, p.trim()))
625 };
626 let start = parse(parts[0], "START")?;
627 let end = parse(parts[1], "END")?;
628 let step = parse(parts[2], "STEP")?;
629 if !step.is_finite() || step <= 0.0 {
630 return Err(format!("Invalid --sweep STEP {}: must be positive", step));
631 }
632 if !start.is_finite() || !end.is_finite() || end < start {
633 return Err(format!(
634 "Invalid --sweep range {}:{}: END must be >= START",
635 start, end
636 ));
637 }
638 let n_f = ((end - start) / step + 1e-9).floor();
644 if !n_f.is_finite() || n_f + 1.0 > MAX_SWEEP_ROWS as f64 {
645 return Err(format!(
646 "--sweep would produce more than {} rows; use a larger STEP",
647 MAX_SWEEP_ROWS
648 ));
649 }
650 let n = n_f as usize + 1;
651 Ok((0..n).map(|i| start + step * i as f64).collect())
653}
654
655pub fn resolve_powder_adjusted_velocity(
664 nominal_velocity_mps: f64,
665 ambient_temperature_c: f64,
666 use_powder_sensitivity: bool,
667 powder_temp_sensitivity_mps_per_c: f64,
668 powder_reference_temp_c: f64,
669 powder_temp_curve: Option<&[(f64, f64)]>,
670 powder_curve_temp_c: Option<f64>,
671) -> f64 {
672 if let Some(curve) = powder_temp_curve {
673 if !curve.is_empty() {
674 let lookup_c = powder_curve_temp_c.unwrap_or(ambient_temperature_c);
675 return interpolate_powder_temp_curve(curve, lookup_c);
676 }
677 return nominal_velocity_mps;
680 }
681 if use_powder_sensitivity {
682 let temp_delta_c = ambient_temperature_c - powder_reference_temp_c;
683 return nominal_velocity_mps + powder_temp_sensitivity_mps_per_c * temp_delta_c;
684 }
685 nominal_velocity_mps
686}
687
688#[derive(Debug, Clone)]
690pub struct WindConditions {
691 pub speed: f64, pub direction: f64,
695 pub vertical_speed: f64,
703}
704
705impl Default for WindConditions {
706 fn default() -> Self {
707 Self {
708 speed: 0.0,
709 direction: 0.0,
710 vertical_speed: 0.0,
711 }
712 }
713}
714
715#[derive(Debug, Clone)]
717pub struct AtmosphericConditions {
718 pub temperature: f64, pub pressure: f64, pub humidity: f64,
724 pub altitude: f64, }
726
727impl Default for AtmosphericConditions {
728 fn default() -> Self {
729 Self {
730 temperature: 15.0,
731 pressure: 1013.25,
732 humidity: 50.0,
733 altitude: 0.0,
734 }
735 }
736}
737
738#[derive(Debug, Clone)]
740pub struct TrajectoryPoint {
741 pub time: f64,
742 pub position: Vector3<f64>,
743 pub velocity_magnitude: f64,
744 pub kinetic_energy: f64,
745 pub drag_coefficient: Option<f64>,
752}
753
754impl TrajectoryPoint {
755 pub fn drag_coefficient_json_value(&self, with_drag_coefficient: bool) -> Option<f64> {
763 if with_drag_coefficient {
764 self.drag_coefficient
765 } else {
766 None
767 }
768 }
769}
770
771#[derive(Debug, Clone)]
773pub struct TrajectoryResult {
774 pub max_range: f64,
775 pub max_height: f64,
776 pub time_of_flight: f64,
777 pub impact_velocity: f64,
778 pub impact_energy: f64,
779 pub projectile_mass_kg: f64,
781 pub line_of_sight_height_m: f64,
783 pub station_speed_of_sound_mps: f64,
785 pub termination: TrajectoryTermination,
787 pub points: Vec<TrajectoryPoint>,
788 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>,
797 pub mach_1_2_distance_m: Option<f64>,
802 pub mach_1_0_distance_m: Option<f64>,
806 pub mach_0_9_distance_m: Option<f64>,
814}
815
816const RK45_TOLERANCE: f64 = 1e-6;
817const RK45_SAFETY_FACTOR: f64 = 0.9;
818const RK45_MAX_DT: f64 = 0.01;
819const RK45_MIN_DT: f64 = 1e-6;
820const TRAJECTORY_TIME_LIMIT_S: f64 = 100.0;
821
822pub const MAX_TRAJECTORY_POINTS: usize = 250_000;
828
829fn cli_rk45_error_norm(
831 position: &Vector3<f64>,
832 velocity: &Vector3<f64>,
833 fifth_position: &Vector3<f64>,
834 fifth_velocity: &Vector3<f64>,
835 fourth_position: &Vector3<f64>,
836 fourth_velocity: &Vector3<f64>,
837) -> f64 {
838 let pack_state = |position: &Vector3<f64>, velocity: &Vector3<f64>| {
839 Vector6::new(
840 position.x, position.y, position.z, velocity.x, velocity.y, velocity.z,
841 )
842 };
843 let state = pack_state(position, velocity);
844 let fifth_order = pack_state(fifth_position, fifth_velocity);
845 let fourth_order = pack_state(fourth_position, fourth_velocity);
846
847 crate::trajectory_integration::rk45_error_norm(&state, &fifth_order, &fourth_order)
848}
849
850struct Rk45Trial {
851 position: Vector3<f64>,
852 velocity: Vector3<f64>,
853 suggested_dt: f64,
854 error: f64,
855}
856
857struct Rk45AcceptedStep {
858 position: Vector3<f64>,
859 velocity: Vector3<f64>,
860 used_dt: f64,
861 next_dt: f64,
862 error: f64,
863}
864
865#[derive(Default)]
879struct MachTransitionTracker {
880 previous_mach: Option<f64>,
881 crossed_transonic: bool,
882 crossed_subsonic: bool,
883 crossed_narrow: bool,
884 mach_1_2_distance_m: Option<f64>,
887 mach_1_0_distance_m: Option<f64>,
890 mach_0_9_distance_m: Option<f64>,
893}
894
895impl MachTransitionTracker {
896 fn record_downward_crossings(&mut self, mach: f64, downrange_m: f64, distances: &mut Vec<f64>) {
897 if !mach.is_finite() {
898 self.previous_mach = None;
899 return;
900 }
901
902 if let Some(previous_mach) = self.previous_mach {
903 if !self.crossed_transonic && previous_mach >= 1.2 && mach < 1.2 {
904 self.crossed_transonic = true;
905 distances.push(downrange_m);
906 self.mach_1_2_distance_m = Some(downrange_m);
907 }
908 if !self.crossed_subsonic && previous_mach >= 1.0 && mach < 1.0 {
909 self.crossed_subsonic = true;
910 distances.push(downrange_m);
911 self.mach_1_0_distance_m = Some(downrange_m);
912 }
913 if !self.crossed_narrow && previous_mach >= 0.9 && mach < 0.9 {
914 self.crossed_narrow = true;
915 self.mach_0_9_distance_m = Some(downrange_m);
917 }
918 }
919 self.previous_mach = Some(mach);
920 }
921}
922
923impl TrajectoryResult {
924 pub fn position_at_range(&self, target_range: f64) -> Option<Vector3<f64>> {
928 if self.points.is_empty() {
929 return None;
930 }
931
932 for i in 0..self.points.len() - 1 {
934 let p1 = &self.points[i];
935 let p2 = &self.points[i + 1];
936
937 if p1.position.x <= target_range && p2.position.x >= target_range {
939 let dx = p2.position.x - p1.position.x;
941 if dx.abs() < 1e-10 {
942 return Some(p1.position);
943 }
944 let t = (target_range - p1.position.x) / dx;
945
946 return Some(Vector3::new(
948 target_range,
949 p1.position.y + t * (p2.position.y - p1.position.y),
950 p1.position.z + t * (p2.position.z - p1.position.z),
951 ));
952 }
953 }
954
955 self.points.last().map(|p| p.position)
957 }
958}
959
960#[derive(Debug, Clone, Copy, PartialEq, Eq)]
962enum StationAtmosphereResolution {
963 LegacyDefaultSentinels,
966 Authoritative,
969}
970
971#[derive(Clone)]
972pub struct TrajectorySolver {
973 inputs: BallisticInputs,
974 wind: WindConditions,
975 atmosphere: AtmosphericConditions,
976 station_atmosphere_resolution: StationAtmosphereResolution,
977 max_range: f64,
978 time_step: f64,
979 max_trajectory_points: usize,
980 cluster_bc: Option<ClusterBCDegradation>,
981 precession_nutation_inertias: (f64, f64),
983 wind_sock: Option<crate::wind::WindSock>,
988 atmo_sock: Option<crate::atmosphere::AtmoSock>,
995}
996
997#[derive(Debug, Clone, Copy, PartialEq, Eq)]
1009pub(crate) enum ZeroTargetFrame {
1010 SightLine,
1011 WorldVertical,
1012}
1013
1014#[derive(Debug, Clone, Copy, PartialEq)]
1027pub struct ZeroCrossings {
1028 pub near_m: Option<f64>,
1031 pub far_m: Option<f64>,
1034}
1035
1036impl TrajectorySolver {
1037 pub fn new(
1038 inputs: BallisticInputs,
1039 wind: WindConditions,
1040 atmosphere: AtmosphericConditions,
1041 ) -> Self {
1042 Self::new_with_station_atmosphere_resolution(
1043 inputs,
1044 wind,
1045 atmosphere,
1046 StationAtmosphereResolution::LegacyDefaultSentinels,
1047 )
1048 }
1049
1050 pub fn new_with_resolved_station_atmosphere(
1061 inputs: BallisticInputs,
1062 wind: WindConditions,
1063 atmosphere: AtmosphericConditions,
1064 ) -> Self {
1065 Self::new_with_station_atmosphere_resolution(
1066 inputs,
1067 wind,
1068 atmosphere,
1069 StationAtmosphereResolution::Authoritative,
1070 )
1071 }
1072
1073 fn new_with_station_atmosphere_resolution(
1074 mut inputs: BallisticInputs,
1075 wind: WindConditions,
1076 atmosphere: AtmosphericConditions,
1077 station_atmosphere_resolution: StationAtmosphereResolution,
1078 ) -> Self {
1079 inputs.normalize_for_solve();
1084
1085 let cluster_bc = if inputs.use_cluster_bc {
1087 Some(ClusterBCDegradation::new())
1088 } else {
1089 None
1090 };
1091 let precession_nutation_inertias = projectile_moments_of_inertia(
1092 inputs.bullet_mass,
1093 inputs.bullet_diameter,
1094 inputs.bullet_length,
1095 );
1096
1097 Self {
1098 inputs,
1099 wind,
1100 atmosphere,
1101 station_atmosphere_resolution,
1102 max_range: 1000.0,
1103 time_step: 0.001,
1104 max_trajectory_points: MAX_TRAJECTORY_POINTS,
1105 cluster_bc,
1106 precession_nutation_inertias,
1107 wind_sock: None,
1108 atmo_sock: None,
1109 }
1110 }
1111
1112 pub fn set_max_range(&mut self, range: f64) {
1113 self.max_range = range;
1114 }
1115
1116 pub fn set_time_step(&mut self, step: f64) {
1117 self.time_step = step;
1118 }
1119
1120 pub(crate) fn calculate_and_set_zero_angle(
1124 &mut self,
1125 target_distance_m: f64,
1126 target_height_m: f64,
1127 frame: ZeroTargetFrame,
1128 ) -> Result<f64, BallisticsError> {
1129 let angle = self.find_zero_angle(target_distance_m, target_height_m, frame)?;
1130 let angle = if target_distance_m > 0.0 {
1139 angle + self.inputs.zero_poi_vertical_m / target_distance_m
1140 } else {
1141 angle
1142 };
1143 self.inputs.muzzle_angle = angle;
1144 self.inputs.azimuth_angle += self.inputs.windage_zero_bias_rad(target_distance_m);
1145 Ok(angle)
1146 }
1147
1148 fn find_zero_angle(
1149 &self,
1150 target_distance_m: f64,
1151 target_height_m: f64,
1152 frame: ZeroTargetFrame,
1153 ) -> Result<f64, BallisticsError> {
1154 let mut low_angle = 0.0;
1157 let mut high_angle = 0.2; let tolerance = 1e-7;
1159 let max_iterations = 60;
1160
1161 let low_height = self.zero_trial_height_at(low_angle, target_distance_m, frame)?;
1163 let high_height = self.zero_trial_height_at(high_angle, target_distance_m, frame)?;
1164
1165 match (low_height, high_height) {
1166 (Some(low_height), Some(high_height)) => {
1167 let low_error = low_height - target_height_m;
1168 let high_error = high_height - target_height_m;
1169
1170 if low_error > 0.0 && high_error > 0.0 {
1171 } else if low_error < 0.0 && high_error < 0.0 {
1174 let mut expanded = false;
1176 for multiplier in [2.0, 3.0, 4.0] {
1177 let new_high = (high_angle * multiplier).min(0.785);
1178 if let Ok(Some(height)) =
1179 self.zero_trial_height_at(new_high, target_distance_m, frame)
1180 {
1181 if height - target_height_m > 0.0 {
1182 high_angle = new_high;
1183 expanded = true;
1184 break;
1185 }
1186 }
1187 if new_high >= 0.785 {
1188 break;
1189 }
1190 }
1191 if !expanded {
1192 return Err("Cannot find zero angle: target beyond effective range even at maximum angle".into());
1193 }
1194 }
1195 }
1196 (None, Some(_)) => {
1197 }
1200 (Some(_), None) => {
1201 return Err(
1202 "Cannot find zero angle: high angle trajectory doesn't reach target distance"
1203 .into(),
1204 );
1205 }
1206 (None, None) => {
1207 return Err(
1208 "Cannot find zero angle: trajectory cannot reach target distance at any angle"
1209 .into(),
1210 );
1211 }
1212 }
1213
1214 for _ in 0..max_iterations {
1215 let mid_angle = (low_angle + high_angle) / 2.0;
1216 match self.zero_trial_height_at(mid_angle, target_distance_m, frame)? {
1217 Some(height) => {
1218 let error = height - target_height_m;
1219 if error.abs() < 0.0001 {
1222 return Ok(mid_angle);
1223 }
1224
1225 if (high_angle - low_angle).abs() < tolerance {
1228 if error.abs() < 0.01 {
1229 return Ok(mid_angle);
1230 }
1231 return Err("Zero angle did not converge: residual height error too large (target not reachable / not bracketed)".into());
1232 }
1233
1234 if error > 0.0 {
1235 high_angle = mid_angle;
1236 } else {
1237 low_angle = mid_angle;
1238 }
1239 }
1240 None => {
1241 low_angle = mid_angle;
1242 if (high_angle - low_angle).abs() < tolerance {
1243 return Err("Trajectory cannot reach target distance - angle converged without valid solution".into());
1244 }
1245 }
1246 }
1247 }
1248
1249 Err("Failed to find zero angle".into())
1250 }
1251
1252 fn zero_trial_height_at(
1255 &self,
1256 angle_rad: f64,
1257 target_distance_m: f64,
1258 frame: ZeroTargetFrame,
1259 ) -> Result<Option<f64>, BallisticsError> {
1260 let mut trial = self.clone();
1261 trial.inputs.muzzle_angle = angle_rad;
1262 trial.inputs.enable_aerodynamic_jump = false;
1265 trial.inputs.cant_angle = 0.0;
1268 if frame == ZeroTargetFrame::SightLine {
1275 trial.inputs.shooting_angle = 0.0;
1276 }
1277 trial.set_max_range(target_distance_m * 2.0);
1278 let result = trial.solve()?;
1279
1280 for (index, point) in result.points.iter().enumerate() {
1281 if point.position.x >= target_distance_m {
1282 let shot_y_m = if index == 0 {
1283 point.position.y
1284 } else {
1285 let previous = &result.points[index - 1];
1286 let span = point.position.x - previous.position.x;
1287 let fraction = (target_distance_m - previous.position.x) / span;
1288 previous.position.y + fraction * (point.position.y - previous.position.y)
1289 };
1290 return Ok(Some(crate::atmosphere::shot_frame_altitude(
1291 0.0,
1292 target_distance_m,
1293 shot_y_m,
1294 trial.inputs.shooting_angle,
1295 )));
1296 }
1297 }
1298 Ok(None)
1299 }
1300
1301 fn find_zero_range(
1328 &self,
1329 angle_rad: f64,
1330 target_height_m: f64,
1331 frame: ZeroTargetFrame,
1332 ) -> Result<ZeroCrossings, BallisticsError> {
1333 let mut trial = self.clone();
1334 trial.inputs.muzzle_angle = angle_rad;
1335 trial.inputs.enable_aerodynamic_jump = false;
1338 trial.inputs.cant_angle = 0.0;
1339 if frame == ZeroTargetFrame::SightLine {
1340 trial.inputs.shooting_angle = 0.0;
1341 }
1342 let result = trial.solve()?;
1343
1344 let mut near_crossing: Option<f64> = None;
1352 let mut far_crossing: Option<f64> = None;
1353 let mut previous: Option<(f64, f64)> = None; for point in &result.points {
1355 let height = crate::atmosphere::shot_frame_altitude(
1356 0.0,
1357 point.position.x,
1358 point.position.y,
1359 trial.inputs.shooting_angle,
1360 );
1361 let error = height - target_height_m;
1362 if let Some((prev_x, prev_error)) = previous {
1363 if prev_error == 0.0 {
1364 if near_crossing.is_none() {
1367 near_crossing = Some(prev_x);
1368 } else {
1369 far_crossing = Some(prev_x);
1370 }
1371 }
1372 if prev_error * error < 0.0 {
1373 let fraction = prev_error / (prev_error - error);
1374 let crossing = prev_x + fraction * (point.position.x - prev_x);
1375 if prev_error < 0.0 && error > 0.0 {
1376 if near_crossing.is_none() {
1378 near_crossing = Some(crossing);
1379 }
1380 } else {
1381 far_crossing = Some(crossing);
1383 }
1384 }
1385 }
1386 previous = Some((point.position.x, error));
1387 }
1388 if let Some((last_x, last_error)) = previous {
1391 if last_error == 0.0 {
1392 if near_crossing.is_none() {
1393 near_crossing = Some(last_x);
1394 } else {
1395 far_crossing = Some(last_x);
1396 }
1397 }
1398 }
1399
1400 if near_crossing.is_none() && far_crossing.is_none() {
1401 return Err(BallisticsError::from(
1402 "Cannot find zero range: this angle never crosses the target height within the \
1403 solved range (angle too shallow to reach it, or both crossings lie beyond \
1404 the solver's max range)."
1405 .to_string(),
1406 ));
1407 }
1408
1409 Ok(ZeroCrossings {
1410 near_m: near_crossing,
1411 far_m: far_crossing,
1412 })
1413 }
1414
1415 pub fn equivalent_horizontal_range(
1439 &self,
1440 target_range_m: f64,
1441 zero_distance_m: f64,
1442 ) -> Option<f64> {
1443 if !target_range_m.is_finite() || !zero_distance_m.is_finite() {
1444 return None;
1445 }
1446 if target_range_m <= zero_distance_m || target_range_m <= 0.0 {
1447 return None;
1448 }
1449
1450 fn path_y_at(points: &[TrajectoryPoint], distance_m: f64) -> Option<f64> {
1455 let mut previous: Option<(f64, f64)> = None;
1456 for point in points {
1457 if point.position.x >= distance_m {
1458 return Some(match previous {
1459 None => point.position.y,
1460 Some((prev_x, prev_y)) => {
1461 let span = point.position.x - prev_x;
1462 if span <= 0.0 {
1463 point.position.y
1464 } else {
1465 let fraction = (distance_m - prev_x) / span;
1466 prev_y + fraction * (point.position.y - prev_y)
1467 }
1468 }
1469 });
1470 }
1471 previous = Some((point.position.x, point.position.y));
1472 }
1473 None
1474 }
1475
1476 let mut inclined = self.clone();
1480 inclined.inputs.enable_trajectory_sampling = false;
1481 let inclined_result = inclined.solve().ok()?;
1482 let los_height = inclined_result.line_of_sight_height_m;
1483 let inclined_drop = los_height - path_y_at(&inclined_result.points, target_range_m)?;
1484 let correction = inclined_drop / target_range_m;
1485 if correction <= 0.0 {
1486 return None;
1487 }
1488
1489 let mut flat = self.clone();
1491 flat.inputs.enable_trajectory_sampling = false;
1492 flat.inputs.shooting_angle = 0.0;
1493 let flat_result = flat.solve().ok()?;
1494 let flat_correction_at = |range_m: f64| -> Option<f64> {
1495 Some((los_height - path_y_at(&flat_result.points, range_m)?) / range_m)
1496 };
1497
1498 let flat_terminal_m = flat_result.points.last().map(|p| p.position.x)?;
1506 let mut low = zero_distance_m.max(1.0);
1507 let mut high = target_range_m.min(flat_terminal_m);
1508 if high <= low {
1509 return None;
1510 }
1511 if flat_correction_at(low)? - correction > 0.0 {
1512 return None; }
1514 if flat_correction_at(high)? - correction < 0.0 {
1515 return None; }
1517 for _ in 0..60 {
1518 let mid = 0.5 * (low + high);
1519 let error = flat_correction_at(mid)? - correction;
1520 if error.abs() == 0.0 {
1521 return Some(mid);
1522 }
1523 if error < 0.0 {
1524 low = mid;
1525 } else {
1526 high = mid;
1527 }
1528 if high - low < 0.01 {
1529 break;
1530 }
1531 }
1532 Some(0.5 * (low + high))
1533 }
1534
1535 fn validate_for_solve(&self) -> Result<(), BallisticsError> {
1541 let require_finite = |name: &str, value: f64| {
1542 if value.is_finite() {
1543 Ok(())
1544 } else {
1545 Err(BallisticsError::from(format!("{name} must be finite")))
1546 }
1547 };
1548 let require_positive = |name: &str, value: f64| {
1549 if value.is_finite() && value > 0.0 {
1550 Ok(())
1551 } else {
1552 Err(BallisticsError::from(format!(
1553 "{name} must be finite and greater than zero"
1554 )))
1555 }
1556 };
1557
1558 if self.inputs.custom_drag_table.is_none() {
1564 require_positive("bc_value", self.inputs.bc_value)?;
1565 }
1566 require_positive("bullet_mass", self.inputs.bullet_mass)?;
1567 require_positive("bullet_diameter", self.inputs.bullet_diameter)?;
1568 require_positive("muzzle_velocity", self.inputs.muzzle_velocity)?;
1569 require_positive("cd_scale", self.inputs.cd_scale)?;
1575
1576 require_finite("muzzle_angle", self.inputs.muzzle_angle)?;
1577 require_finite("azimuth_angle", self.inputs.azimuth_angle)?;
1578 require_finite("shooting_angle", self.inputs.shooting_angle)?;
1579 require_finite("cant_angle", self.inputs.cant_angle)?;
1580 require_finite("muzzle_height", self.inputs.muzzle_height)?;
1581
1582 for (name, value) in [
1586 ("zero_poi_vertical_m", self.inputs.zero_poi_vertical_m),
1587 ("zero_poi_horizontal_m", self.inputs.zero_poi_horizontal_m),
1588 ] {
1589 require_finite(name, value)?;
1590 if value.abs() >= 1.0 {
1591 return Err(BallisticsError::from(format!(
1592 "{name} must be smaller than 1.0 m in magnitude (it is a linear POI \
1593 offset at the zero range, in meters)"
1594 )));
1595 }
1596 }
1597
1598 require_finite(
1602 "sight_offset_lateral_m",
1603 self.inputs.sight_offset_lateral_m,
1604 )?;
1605 if self.inputs.sight_offset_lateral_m.abs() >= 0.5 {
1606 return Err(BallisticsError::from(
1607 "sight_offset_lateral_m must be smaller than 0.5 m in magnitude (it is \
1608 the lateral sight-to-bore mount offset, in meters)",
1609 ));
1610 }
1611
1612 if !(self.inputs.ground_threshold.is_finite()
1615 || self.inputs.ground_threshold == f64::NEG_INFINITY)
1616 {
1617 return Err(BallisticsError::from(
1618 "ground_threshold must be finite or negative infinity",
1619 ));
1620 }
1621
1622 match &self.wind_sock {
1623 Some(wind_sock) => wind_sock
1624 .validate_segments()
1625 .map_err(BallisticsError::from)?,
1626 None => {
1627 require_finite("wind.speed", self.wind.speed)?;
1628 require_finite("wind.direction", self.wind.direction)?;
1629 require_finite("wind.vertical_speed", self.wind.vertical_speed)?;
1630 }
1631 }
1632
1633 require_finite("atmosphere.temperature", self.atmosphere.temperature)?;
1634 require_finite("atmosphere.pressure", self.atmosphere.pressure)?;
1635 require_finite("atmosphere.humidity", self.atmosphere.humidity)?;
1636 require_finite("atmosphere.altitude", self.atmosphere.altitude)?;
1637
1638 require_positive("max_range", self.max_range)?;
1639 if !self.inputs.use_rk4 || !self.inputs.use_adaptive_rk45 {
1642 require_positive("time_step", self.time_step)?;
1643 }
1644
1645 if self.inputs.enable_trajectory_sampling {
1646 require_finite("sight_height", self.inputs.sight_height)?;
1647 require_positive("sample_interval", self.inputs.sample_interval)?;
1648 projected_sample_count(self.max_range, self.inputs.sample_interval)?;
1649 }
1650
1651 if self.inputs.drops_reference == DropsReference::Target {
1655 require_finite("target_height", self.inputs.target_height)?;
1656 if self.inputs.shooting_angle.cos() <= 1e-9 {
1657 return Err(BallisticsError::from(
1658 "drops reference 'target' is undefined for shooting angles at or beyond 90 degrees",
1659 ));
1660 }
1661 }
1662
1663 if self.inputs.enable_coriolis {
1664 require_finite("shot_azimuth", self.inputs.shot_azimuth)?;
1665 if let Some(latitude) = self.inputs.latitude {
1666 require_finite("latitude", latitude)?;
1667 }
1668 }
1669
1670 Ok(())
1671 }
1672
1673 fn validate_result_sanity(&self, result: &TrajectoryResult) -> Result<(), BallisticsError> {
1680 let require_finite = |name: &str, value: f64| {
1681 if value.is_finite() {
1682 Ok(())
1683 } else {
1684 Err(BallisticsError::from(format!(
1685 "trajectory result contains non-finite {name}"
1686 )))
1687 }
1688 };
1689 let require_non_negative = |name: &str, value: f64| {
1690 if value >= 0.0 {
1691 Ok(())
1692 } else {
1693 Err(BallisticsError::from(format!(
1694 "trajectory result contains non-physical negative {name} ({value})"
1695 )))
1696 }
1697 };
1698 let require_indexed_finite = |collection: &str, index: usize, field: &str, value: f64| {
1699 if value.is_finite() {
1700 Ok(())
1701 } else {
1702 Err(BallisticsError::from(format!(
1703 "trajectory result contains non-finite {collection}[{index}].{field}"
1704 )))
1705 }
1706 };
1707 let require_indexed_non_negative =
1708 |collection: &str, index: usize, field: &str, value: f64| {
1709 if value >= 0.0 {
1710 Ok(())
1711 } else {
1712 Err(BallisticsError::from(format!(
1713 "trajectory result contains non-physical negative {collection}[{index}].{field} ({value})"
1714 )))
1715 }
1716 };
1717
1718 require_finite("max_range", result.max_range)?;
1719 require_finite("max_height", result.max_height)?;
1720 require_finite("time_of_flight", result.time_of_flight)?;
1721 require_finite("impact_velocity", result.impact_velocity)?;
1722 require_finite("impact_energy", result.impact_energy)?;
1723 require_finite("projectile_mass_kg", result.projectile_mass_kg)?;
1724 require_finite(
1725 "line_of_sight_height_m",
1726 result.line_of_sight_height_m,
1727 )?;
1728 require_finite(
1729 "station_speed_of_sound_mps",
1730 result.station_speed_of_sound_mps,
1731 )?;
1732
1733 require_non_negative("max_range", result.max_range)?;
1737 require_non_negative("time_of_flight", result.time_of_flight)?;
1738 require_non_negative("impact_velocity", result.impact_velocity)?;
1739 require_non_negative("impact_energy", result.impact_energy)?;
1740 require_non_negative("projectile_mass_kg", result.projectile_mass_kg)?;
1741 require_non_negative(
1742 "station_speed_of_sound_mps",
1743 result.station_speed_of_sound_mps,
1744 )?;
1745
1746 for (index, point) in result.points.iter().enumerate() {
1747 require_indexed_finite("points", index, "time", point.time)?;
1748 require_indexed_finite("points", index, "position.x", point.position.x)?;
1749 require_indexed_finite("points", index, "position.y", point.position.y)?;
1750 require_indexed_finite("points", index, "position.z", point.position.z)?;
1751 require_indexed_finite(
1752 "points",
1753 index,
1754 "velocity_magnitude",
1755 point.velocity_magnitude,
1756 )?;
1757 require_indexed_finite("points", index, "kinetic_energy", point.kinetic_energy)?;
1758 require_indexed_non_negative("points", index, "time", point.time)?;
1759 require_indexed_non_negative(
1760 "points",
1761 index,
1762 "velocity_magnitude",
1763 point.velocity_magnitude,
1764 )?;
1765 require_indexed_non_negative("points", index, "kinetic_energy", point.kinetic_energy)?;
1766 }
1767
1768 if let Some(samples) = &result.sampled_points {
1769 for (index, sample) in samples.iter().enumerate() {
1770 require_indexed_finite("sampled_points", index, "distance_m", sample.distance_m)?;
1771 require_indexed_finite("sampled_points", index, "drop_m", sample.drop_m)?;
1772 require_indexed_finite(
1773 "sampled_points",
1774 index,
1775 "wind_drift_m",
1776 sample.wind_drift_m,
1777 )?;
1778 require_indexed_finite(
1779 "sampled_points",
1780 index,
1781 "velocity_mps",
1782 sample.velocity_mps,
1783 )?;
1784 require_indexed_finite("sampled_points", index, "energy_j", sample.energy_j)?;
1785 require_indexed_finite("sampled_points", index, "time_s", sample.time_s)?;
1786 }
1787 }
1788
1789 for (name, value) in [
1790 ("min_pitch_damping", result.min_pitch_damping),
1791 ("transonic_mach", result.transonic_mach),
1792 ("max_yaw_angle", result.max_yaw_angle),
1793 ("max_precession_angle", result.max_precession_angle),
1794 ] {
1795 if let Some(value) = value {
1796 require_finite(name, value)?;
1797 }
1798 }
1799
1800 if let Some(state) = result.angular_state {
1801 for (name, value) in [
1802 ("angular_state.pitch_angle", state.pitch_angle),
1803 ("angular_state.yaw_angle", state.yaw_angle),
1804 ("angular_state.pitch_rate", state.pitch_rate),
1805 ("angular_state.yaw_rate", state.yaw_rate),
1806 ("angular_state.precession_angle", state.precession_angle),
1807 ("angular_state.nutation_phase", state.nutation_phase),
1808 ] {
1809 require_finite(name, value)?;
1810 }
1811 }
1812
1813 if let Some(jump) = result.aerodynamic_jump {
1814 for (name, value) in [
1815 ("aerodynamic_jump.vertical_jump_moa", jump.vertical_jump_moa),
1816 (
1817 "aerodynamic_jump.horizontal_jump_moa",
1818 jump.horizontal_jump_moa,
1819 ),
1820 ("aerodynamic_jump.jump_angle_rad", jump.jump_angle_rad),
1821 (
1822 "aerodynamic_jump.magnus_component_moa",
1823 jump.magnus_component_moa,
1824 ),
1825 ("aerodynamic_jump.yaw_component_moa", jump.yaw_component_moa),
1826 (
1827 "aerodynamic_jump.stabilization_factor",
1828 jump.stabilization_factor,
1829 ),
1830 ] {
1831 require_finite(name, value)?;
1832 }
1833 }
1834
1835 Ok(())
1836 }
1837
1838 fn validate_integration_state(
1851 &self,
1852 position: &Vector3<f64>,
1853 velocity: &Vector3<f64>,
1854 time: f64,
1855 ) -> Result<(), BallisticsError> {
1856 if !(position.iter().all(|value| value.is_finite())
1857 && velocity.iter().all(|value| value.is_finite())
1858 && time.is_finite())
1859 {
1860 return Err(BallisticsError::from(
1861 "trajectory integration produced a non-finite state (often from physically \
1862 extreme inputs — e.g. an absurd bore/muzzle height placing the launch far \
1863 from sea level, or a degenerate atmosphere; check those inputs, or set \
1864 --altitude explicitly)",
1865 ));
1866 }
1867
1868 let speed = velocity.magnitude();
1869 let budget = self.speed_budget(time);
1870 if speed > budget {
1871 return Err(BallisticsError::from(format!(
1872 "trajectory integration diverged: speed {speed:.3e} m/s at t={time:.6}s exceeds \
1873 the physical budget of {budget:.3e} m/s"
1874 )));
1875 }
1876 Ok(())
1877 }
1878
1879 fn speed_budget(&self, time: f64) -> f64 {
1884 let scalar_wind = self.wind.speed.abs() + self.wind.vertical_speed.abs();
1885 let wind_bound = match &self.wind_sock {
1886 Some(sock) => scalar_wind.max(sock.max_speed_mps()),
1887 None => scalar_wind,
1888 };
1889 2.0 * (self.inputs.muzzle_velocity + wind_bound + 10.0)
1890 + crate::constants::G_ACCEL_MPS2 * time
1891 }
1892
1893 fn push_trajectory_point(
1895 &self,
1896 points: &mut Vec<TrajectoryPoint>,
1897 point: TrajectoryPoint,
1898 ) -> Result<(), BallisticsError> {
1899 if points.len() >= self.max_trajectory_points {
1900 return Err(BallisticsError::from(format!(
1901 "trajectory point limit of {} exceeded",
1902 self.max_trajectory_points
1903 )));
1904 }
1905 points.push(point);
1906 Ok(())
1907 }
1908
1909 pub fn set_wind_segments(&mut self, segments: Vec<crate::wind::WindSegment>) {
1916 self.wind_sock = if segments.is_empty() {
1917 None
1918 } else {
1919 Some(crate::wind::WindSock::new(segments))
1920 };
1921 }
1922
1923 pub fn set_atmo_segments(&mut self, segments: Vec<crate::atmosphere::AtmoSegment>) {
1932 self.atmo_sock = if segments.is_empty() {
1933 None
1934 } else {
1935 Some(crate::atmosphere::AtmoSock::new(segments))
1936 };
1937 }
1938
1939 fn launch_angles_from(
1948 &self,
1949 aj: Option<&crate::aerodynamic_jump::AerodynamicJumpComponents>,
1950 ) -> (f64, f64) {
1951 let (mut elev, mut azim) = (self.inputs.muzzle_angle, self.inputs.azimuth_angle);
1952 if self.inputs.cant_angle != 0.0 {
1959 let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
1960 let (e0, a0) = (elev, azim);
1961 elev = e0 * cos_c - a0 * sin_c;
1962 azim = a0 * cos_c + e0 * sin_c;
1963 }
1964 match aj {
1965 Some(c) => {
1966 const MOA_PER_RAD: f64 = 3437.7467707849;
1968 (
1969 elev + c.vertical_jump_moa / MOA_PER_RAD,
1970 azim + c.horizontal_jump_moa / MOA_PER_RAD,
1971 )
1972 }
1973 None => (elev, azim),
1974 }
1975 }
1976
1977 fn aerodynamic_jump_components(
1985 &self,
1986 ) -> Option<crate::aerodynamic_jump::AerodynamicJumpComponents> {
1987 if !self.inputs.enable_aerodynamic_jump {
1988 return None;
1989 }
1990 let diameter_m = self.inputs.bullet_diameter;
1994 if !(self.inputs.twist_rate.is_finite()
1995 && self.inputs.twist_rate != 0.0
1996 && diameter_m.is_finite()
1997 && diameter_m > 0.0
1998 && self.inputs.bullet_length.is_finite()
1999 && self.inputs.bullet_length > 0.0
2000 && self.inputs.muzzle_velocity.is_finite())
2001 {
2002 return None;
2003 }
2004
2005 let (_, _, temp_c, pressure_hpa) = self.resolved_atmosphere();
2007 let sg = crate::stability::compute_stability_coefficient(
2008 &self.inputs,
2009 (self.atmosphere.altitude, temp_c, pressure_hpa, 0.0),
2010 );
2011 if !(sg.is_finite() && sg > 0.0) {
2012 return None;
2013 }
2014 let length_calibers = self.inputs.bullet_length / diameter_m;
2015
2016 const MS_TO_MPH: f64 = 2.236_936_292_054_4;
2022 let crosswind_from_right_mps = if let Some(sock) = &self.wind_sock {
2023 -sock.vector_for_range_stateless(0.0)[2]
2024 } else {
2025 self.wind.speed * self.wind.direction.sin()
2026 };
2027 let crosswind_from_right_mph = crosswind_from_right_mps * MS_TO_MPH;
2028
2029 let vertical_jump_moa = crate::aerodynamic_jump::litz_crosswind_jump_moa(
2030 sg,
2031 length_calibers,
2032 crosswind_from_right_mph,
2033 self.inputs.is_twist_right,
2034 );
2035 if !vertical_jump_moa.is_finite() {
2036 return None;
2037 }
2038
2039 const MOA_PER_RAD: f64 = 3437.7467707849;
2040 Some(crate::aerodynamic_jump::AerodynamicJumpComponents {
2041 vertical_jump_moa,
2042 horizontal_jump_moa: 0.0,
2044 jump_angle_rad: vertical_jump_moa.abs() / MOA_PER_RAD,
2045 magnus_component_moa: 0.0,
2046 yaw_component_moa: 0.0,
2047 stabilization_factor: (sg / 1.5).clamp(0.0, 1.0),
2048 })
2049 }
2050
2051 fn resolved_atmosphere(&self) -> (f64, f64, f64, f64) {
2052 let (temp_c, pressure_hpa) = match self.station_atmosphere_resolution {
2053 StationAtmosphereResolution::LegacyDefaultSentinels => {
2054 crate::atmosphere::resolve_station_conditions(
2055 self.atmosphere.temperature,
2056 self.atmosphere.pressure,
2057 self.atmosphere.altitude,
2058 )
2059 }
2060 StationAtmosphereResolution::Authoritative => {
2061 (self.atmosphere.temperature, self.atmosphere.pressure)
2062 }
2063 };
2064 let (density, speed_of_sound) = crate::atmosphere::calculate_atmosphere(
2065 self.atmosphere.altitude,
2066 Some(temp_c),
2067 Some(pressure_hpa),
2068 self.atmosphere.humidity,
2069 );
2070 (density, speed_of_sound, temp_c, pressure_hpa)
2071 }
2072
2073 fn precession_nutation_params(
2074 &self,
2075 velocity_mps: f64,
2076 air_density_kg_m3: f64,
2077 speed_of_sound_mps: f64,
2078 ) -> PrecessionNutationParams {
2079 let (spin_inertia, transverse_inertia) = self.precession_nutation_inertias;
2080 let spin_rate_rad_s = if self.inputs.twist_rate > 0.0 {
2081 let velocity_fps = velocity_mps * 3.28084;
2082 let twist_rate_ft = self.inputs.twist_rate / 12.0;
2083 (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI
2084 } else {
2085 0.0
2086 };
2087
2088 PrecessionNutationParams {
2089 mass_kg: self.inputs.bullet_mass,
2090 caliber_m: self.inputs.bullet_diameter,
2091 length_m: self.inputs.bullet_length,
2092 spin_rate_rad_s,
2093 spin_inertia,
2094 transverse_inertia,
2095 velocity_mps,
2096 air_density_kg_m3,
2097 mach: velocity_mps / speed_of_sound_mps,
2098 pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
2099 nutation_damping_factor: 0.05,
2100 }
2101 }
2102
2103 fn append_terminal_endpoint(
2110 &self,
2111 points: &mut Vec<TrajectoryPoint>,
2112 post_position: Vector3<f64>,
2113 post_velocity: Vector3<f64>,
2114 post_time: f64,
2115 max_height: &mut f64,
2116 ) -> Result<TrajectoryTermination, BallisticsError> {
2117 let previous = points
2118 .last()
2119 .cloned()
2120 .ok_or_else(|| BallisticsError::from("No trajectory points generated"))?;
2121
2122 let mut crossings = Vec::with_capacity(3);
2123 if previous.position.x < self.max_range && post_position.x >= self.max_range {
2124 let span = post_position.x - previous.position.x;
2125 if span.is_finite() && span > 0.0 {
2126 crossings.push((
2127 (self.max_range - previous.position.x) / span,
2128 TrajectoryTermination::MaxRange,
2129 ));
2130 }
2131 }
2132 if self.inputs.ground_threshold.is_finite()
2133 && previous.position.y > self.inputs.ground_threshold
2134 && post_position.y <= self.inputs.ground_threshold
2135 {
2136 let span = post_position.y - previous.position.y;
2137 if span.is_finite() && span < 0.0 {
2138 crossings.push((
2139 (self.inputs.ground_threshold - previous.position.y) / span,
2140 TrajectoryTermination::GroundThreshold,
2141 ));
2142 }
2143 }
2144 if previous.time < TRAJECTORY_TIME_LIMIT_S && post_time >= TRAJECTORY_TIME_LIMIT_S {
2145 let span = post_time - previous.time;
2146 if span.is_finite() && span > 0.0 {
2147 crossings.push((
2148 (TRAJECTORY_TIME_LIMIT_S - previous.time) / span,
2149 TrajectoryTermination::TimeLimit,
2150 ));
2151 }
2152 }
2153
2154 let (fraction, termination) = crossings
2155 .into_iter()
2156 .filter(|(fraction, _)| fraction.is_finite() && (0.0..=1.0).contains(fraction))
2157 .min_by(|left, right| {
2158 let priority = |termination: TrajectoryTermination| match termination {
2159 TrajectoryTermination::GroundThreshold => 0,
2160 TrajectoryTermination::MaxRange => 1,
2161 TrajectoryTermination::TimeLimit => 2,
2162 TrajectoryTermination::VelocityFloor => 3,
2163 };
2164 left.0
2165 .total_cmp(&right.0)
2166 .then_with(|| priority(left.1).cmp(&priority(right.1)))
2167 })
2168 .ok_or_else(|| {
2169 BallisticsError::from(
2170 "trajectory integration stopped without crossing a supported boundary",
2171 )
2172 })?;
2173
2174 let mut position = previous.position + (post_position - previous.position) * fraction;
2175 match termination {
2176 TrajectoryTermination::MaxRange => position.x = self.max_range,
2177 TrajectoryTermination::GroundThreshold => {
2178 position.y = self.inputs.ground_threshold;
2179 }
2180 TrajectoryTermination::TimeLimit | TrajectoryTermination::VelocityFloor => {}
2181 }
2182 let velocity_magnitude = previous.velocity_magnitude
2183 + (post_velocity.magnitude() - previous.velocity_magnitude) * fraction;
2184 let mut time = previous.time + (post_time - previous.time) * fraction;
2185 if termination == TrajectoryTermination::TimeLimit {
2186 time = TRAJECTORY_TIME_LIMIT_S;
2187 }
2188 let kinetic_energy =
2189 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2190
2191 if position.y > *max_height {
2192 *max_height = position.y;
2193 }
2194 let terminal_point = TrajectoryPoint {
2195 time,
2196 position,
2197 velocity_magnitude,
2198 kinetic_energy,
2199 drag_coefficient: None,
2200 };
2201 if terminal_point.position.x < previous.position.x {
2202 return Err(BallisticsError::from(
2203 "trajectory terminal state reversed downrange before the crossed boundary",
2204 ));
2205 }
2206 if terminal_point.position.x == previous.position.x {
2207 let last = points.last_mut().ok_or_else(|| {
2212 BallisticsError::from("trajectory points disappeared during terminal finalization")
2213 })?;
2214 *last = terminal_point;
2215 } else {
2216 self.push_trajectory_point(points, terminal_point)?;
2217 }
2218 Ok(termination)
2219 }
2220
2221 fn gravity_acceleration(&self) -> Vector3<f64> {
2222 let theta = self.inputs.shooting_angle;
2223 Vector3::new(
2224 -crate::constants::G_ACCEL_MPS2 * theta.sin(),
2225 -crate::constants::G_ACCEL_MPS2 * theta.cos(),
2226 0.0,
2227 )
2228 }
2229
2230 fn get_wind_at_altitude(&self, altitude_m: f64) -> Vector3<f64> {
2231 let model = match self.inputs.wind_shear_model.as_str() {
2246 "logarithmic" => WindShearModel::Logarithmic,
2247 "power_law" | "powerlaw" | "exponential" => WindShearModel::PowerLaw,
2248 "ekman_spiral" | "ekman" => WindShearModel::EkmanSpiral,
2249 "custom_layers" | "custom" => WindShearModel::CustomLayers,
2250 _ => WindShearModel::PowerLaw,
2251 };
2252 let speed_ratio = crate::wind_shear::boundary_layer_speed_ratio(altitude_m, model);
2253
2254 crate::wind::wind_vector(self.wind.speed, self.wind.direction, 0.0) * speed_ratio
2261 + Vector3::new(0.0, self.wind.vertical_speed, 0.0)
2262 }
2263
2264 pub fn solve(&self) -> Result<TrajectoryResult, BallisticsError> {
2265 self.validate_for_solve()?;
2266 let mut result = if self.inputs.use_rk4 {
2267 if self.inputs.use_adaptive_rk45 {
2268 self.solve_rk45()?
2269 } else {
2270 self.solve_rk4()?
2271 }
2272 } else {
2273 self.solve_euler()?
2274 };
2275 self.apply_spin_drift(&mut result);
2276 self.validate_result_sanity(&result)?;
2277 Ok(result)
2278 }
2279
2280 fn apply_spin_drift(&self, result: &mut TrajectoryResult) {
2286 if !self.inputs.use_enhanced_spin_drift {
2287 return;
2288 }
2289 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 {
2293 return;
2294 }
2295
2296 let sg = self.effective_spin_drift_sg();
2303
2304 for p in result.points.iter_mut() {
2305 if p.time <= 0.0 {
2306 continue;
2307 }
2308 p.position.z +=
2310 crate::spin_drift::litz_drift_meters(sg, p.time, self.inputs.is_twist_right);
2311 }
2312
2313 if let Some(samples) = result.sampled_points.as_mut() {
2317 for s in samples.iter_mut() {
2318 if s.time_s <= 0.0 {
2319 continue;
2320 }
2321 s.wind_drift_m +=
2322 crate::spin_drift::litz_drift_meters(sg, s.time_s, self.inputs.is_twist_right);
2323 }
2324 }
2325 }
2326
2327 fn effective_spin_drift_sg(&self) -> f64 {
2332 let (_, _, temp_c, press_hpa) = self.resolved_atmosphere();
2333 crate::spin_drift::effective_sg_from_inputs(&self.inputs, temp_c, press_hpa)
2334 }
2335
2336 fn initial_position(&self) -> Vector3<f64> {
2354 if self.inputs.cant_angle == 0.0 && self.inputs.sight_offset_lateral_m == 0.0 {
2355 return Vector3::new(0.0, self.inputs.muzzle_height, 0.0);
2356 }
2357 let (sin_c, cos_c) = self.inputs.cant_angle.sin_cos();
2358 let sh = self.inputs.sight_height;
2359 let off = self.inputs.sight_offset_lateral_m;
2360 Vector3::new(
2361 0.0,
2362 self.inputs.muzzle_height + sh * (1.0 - cos_c) + off * sin_c,
2363 -sh * sin_c - off * cos_c,
2364 )
2365 }
2366
2367 fn build_sampled_points(
2374 &self,
2375 points: &[TrajectoryPoint],
2376 max_height: f64,
2377 transonic_distances: Vec<f64>,
2378 mach_transitions: &MachTransitionTracker,
2379 ) -> Result<Option<Vec<TrajectorySample>>, BallisticsError> {
2380 if !self.inputs.enable_trajectory_sampling {
2381 return Ok(None);
2382 }
2383
2384 let last_point = points.last().ok_or("No trajectory points generated")?;
2385 let trajectory_data = TrajectoryData {
2386 times: points.iter().map(|p| p.time).collect(),
2387 positions: points.iter().map(|p| p.position).collect(),
2388 velocities: points
2389 .iter()
2390 .map(|p| {
2391 Vector3::new(0.0, 0.0, p.velocity_magnitude)
2393 })
2394 .collect(),
2395 transonic_distances, mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2397 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2398 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2399 };
2400
2401 let sight_position_m = self.inputs.muzzle_height + self.inputs.sight_height;
2406 let target_reference = self.inputs.drops_reference == DropsReference::Target;
2407 let target_vertical_height_m = if target_reference && self.inputs.target_height != 0.0 {
2414 self.inputs.target_height
2415 } else {
2416 sight_position_m
2417 };
2418 let outputs = TrajectoryOutputs {
2419 target_distance_horiz_m: last_point.position.x, target_vertical_height_m,
2421 time_of_flight_s: last_point.time,
2422 max_ord_dist_horiz_m: max_height,
2423 sight_height_m: sight_position_m,
2424 };
2425
2426 let mut samples = sample_trajectory(
2428 &trajectory_data,
2429 &outputs,
2430 self.inputs.sample_interval,
2431 self.inputs.bullet_mass,
2432 )?;
2433
2434 if target_reference {
2441 let cos_theta = self.inputs.shooting_angle.cos();
2442 for sample in &mut samples {
2443 sample.drop_m /= cos_theta;
2444 }
2445 }
2446 Ok(Some(samples))
2447 }
2448
2449 fn solve_euler(&self) -> Result<TrajectoryResult, BallisticsError> {
2450 let mut time = 0.0;
2452 let mut position = self.initial_position();
2456 let aj_components = self.aerodynamic_jump_components();
2462 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2463 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2464 let mut velocity = Vector3::new(
2465 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2469
2470 let mut points = Vec::new();
2471 let mut max_height = position.y;
2472 let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
2478 let mut mach_transitions = MachTransitionTracker::default();
2479
2480 let mut angular_state = if self.inputs.enable_precession_nutation {
2482 Some(AngularState {
2483 pitch_angle: 0.001, yaw_angle: 0.001,
2485 pitch_rate: 0.0,
2486 yaw_rate: 0.0,
2487 precession_angle: 0.0,
2488 nutation_phase: 0.0,
2489 })
2490 } else {
2491 None
2492 };
2493 let mut max_yaw_angle = 0.0;
2494 let mut max_precession_angle = 0.0;
2495
2496 let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2498 self.resolved_atmosphere();
2499 let base_ratio = air_density / 1.225;
2504
2505 let wind_vector =
2511 crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
2512
2513 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2516 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2517 );
2518
2519 while position.x < self.max_range
2521 && position.y > self.inputs.ground_threshold
2522 && time < TRAJECTORY_TIME_LIMIT_S
2523 {
2524 let velocity_magnitude = velocity.magnitude();
2526 let kinetic_energy =
2527 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2528
2529 self.push_trajectory_point(
2530 &mut points,
2531 TrajectoryPoint {
2532 time,
2533 position,
2534 velocity_magnitude,
2535 kinetic_energy,
2536 drag_coefficient: None,
2537 },
2538 )?;
2539
2540 {
2543 let mach_here = if speed_of_sound > 0.0 {
2544 velocity_magnitude / speed_of_sound
2545 } else {
2546 0.0
2547 };
2548 mach_transitions.record_downward_crossings(
2549 mach_here,
2550 position.x,
2551 &mut transonic_distances,
2552 );
2553 }
2554
2555 if position.y > max_height {
2557 max_height = position.y;
2558 }
2559
2560 if self.inputs.enable_pitch_damping {
2562 let mach = velocity_magnitude / speed_of_sound;
2563
2564 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
2566 transonic_mach = Some(mach);
2567 }
2568
2569 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
2571
2572 if pitch_damping < min_pitch_damping {
2574 min_pitch_damping = pitch_damping;
2575 }
2576 }
2577
2578 if self.inputs.enable_precession_nutation {
2580 if let Some(ref mut state) = angular_state {
2581 let velocity_magnitude = velocity.magnitude();
2582 let params = self.precession_nutation_params(
2583 velocity_magnitude,
2584 air_density,
2585 speed_of_sound,
2586 );
2587
2588 *state = calculate_combined_angular_motion(
2590 ¶ms,
2591 state,
2592 time,
2593 self.time_step,
2594 0.001, );
2596
2597 if state.yaw_angle.abs() > max_yaw_angle {
2599 max_yaw_angle = state.yaw_angle.abs();
2600 }
2601 if state.precession_angle.abs() > max_precession_angle {
2602 max_precession_angle = state.precession_angle.abs();
2603 }
2604 }
2605 }
2606
2607 let acceleration = self.calculate_acceleration(
2614 &position,
2615 &velocity,
2616 &wind_vector,
2617 (resolved_temp_c, resolved_press_hpa, base_ratio),
2618 );
2619
2620 velocity += acceleration * self.time_step;
2622 position += velocity * self.time_step;
2623 time += self.time_step;
2624 self.validate_integration_state(&position, &velocity, time)?;
2625 }
2626
2627 let termination =
2628 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
2629
2630 self.annotate_drag_coefficients(&mut points, speed_of_sound);
2635
2636 let last_point = points.last().ok_or("No trajectory points generated")?;
2637
2638 let sampled_points = self.build_sampled_points(
2640 &points,
2641 max_height,
2642 transonic_distances,
2643 &mach_transitions,
2644 )?;
2645
2646 Ok(TrajectoryResult {
2647 max_range: last_point.position.x, max_height,
2649 time_of_flight: last_point.time,
2650 impact_velocity: last_point.velocity_magnitude,
2651 impact_energy: last_point.kinetic_energy,
2652 projectile_mass_kg: self.inputs.bullet_mass,
2653 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
2654 station_speed_of_sound_mps: speed_of_sound,
2655 termination,
2656 points,
2657 sampled_points,
2658 min_pitch_damping: if self.inputs.enable_pitch_damping {
2659 Some(min_pitch_damping)
2660 } else {
2661 None
2662 },
2663 transonic_mach,
2664 angular_state,
2665 max_yaw_angle: if self.inputs.enable_precession_nutation {
2666 Some(max_yaw_angle)
2667 } else {
2668 None
2669 },
2670 max_precession_angle: if self.inputs.enable_precession_nutation {
2671 Some(max_precession_angle)
2672 } else {
2673 None
2674 },
2675 aerodynamic_jump: aj_components,
2676 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2677 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2678 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2679 })
2680 }
2681
2682 fn solve_rk4(&self) -> Result<TrajectoryResult, BallisticsError> {
2683 let mut time = 0.0;
2685 let mut position = self.initial_position();
2690
2691 let aj_components = self.aerodynamic_jump_components();
2697 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2698 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2699 let mut velocity = Vector3::new(
2700 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2704
2705 let mut points = Vec::new();
2706 let mut max_height = position.y;
2707 let mut min_pitch_damping = f64::INFINITY; let mut transonic_mach = None; let mut transonic_distances: Vec<f64> = Vec::new();
2713 let mut mach_transitions = MachTransitionTracker::default();
2714
2715 let mut angular_state = if self.inputs.enable_precession_nutation {
2717 Some(AngularState {
2718 pitch_angle: 0.001, yaw_angle: 0.001,
2720 pitch_rate: 0.0,
2721 yaw_rate: 0.0,
2722 precession_angle: 0.0,
2723 nutation_phase: 0.0,
2724 })
2725 } else {
2726 None
2727 };
2728 let mut max_yaw_angle = 0.0;
2729 let mut max_precession_angle = 0.0;
2730
2731 let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2733 self.resolved_atmosphere();
2734 let base_ratio = air_density / 1.225;
2739
2740 let wind_vector =
2746 crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
2747
2748 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2751 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2752 );
2753
2754 while position.x < self.max_range
2756 && position.y > self.inputs.ground_threshold
2757 && time < TRAJECTORY_TIME_LIMIT_S
2758 {
2759 let velocity_magnitude = velocity.magnitude();
2761 let kinetic_energy =
2762 0.5 * self.inputs.bullet_mass * velocity_magnitude * velocity_magnitude;
2763
2764 self.push_trajectory_point(
2765 &mut points,
2766 TrajectoryPoint {
2767 time,
2768 position,
2769 velocity_magnitude,
2770 kinetic_energy,
2771 drag_coefficient: None,
2772 },
2773 )?;
2774
2775 {
2778 let mach_here = if speed_of_sound > 0.0 {
2779 velocity_magnitude / speed_of_sound
2780 } else {
2781 0.0
2782 };
2783 mach_transitions.record_downward_crossings(
2784 mach_here,
2785 position.x,
2786 &mut transonic_distances,
2787 );
2788 }
2789
2790 if position.y > max_height {
2791 max_height = position.y;
2792 }
2793
2794 if self.inputs.enable_pitch_damping {
2796 let mach = velocity_magnitude / speed_of_sound;
2797
2798 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
2800 transonic_mach = Some(mach);
2801 }
2802
2803 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
2805
2806 if pitch_damping < min_pitch_damping {
2808 min_pitch_damping = pitch_damping;
2809 }
2810 }
2811
2812 if self.inputs.enable_precession_nutation {
2814 if let Some(ref mut state) = angular_state {
2815 let velocity_magnitude = velocity.magnitude();
2816 let params = self.precession_nutation_params(
2817 velocity_magnitude,
2818 air_density,
2819 speed_of_sound,
2820 );
2821
2822 *state = calculate_combined_angular_motion(
2824 ¶ms,
2825 state,
2826 time,
2827 self.time_step,
2828 0.001, );
2830
2831 if state.yaw_angle.abs() > max_yaw_angle {
2833 max_yaw_angle = state.yaw_angle.abs();
2834 }
2835 if state.precession_angle.abs() > max_precession_angle {
2836 max_precession_angle = state.precession_angle.abs();
2837 }
2838 }
2839 }
2840
2841 let dt = self.time_step;
2843
2844 let acc1 = self.calculate_acceleration(
2846 &position,
2847 &velocity,
2848 &wind_vector,
2849 (resolved_temp_c, resolved_press_hpa, base_ratio),
2850 );
2851
2852 let pos2 = position + velocity * (dt * 0.5);
2854 let vel2 = velocity + acc1 * (dt * 0.5);
2855 let acc2 = self.calculate_acceleration(
2856 &pos2,
2857 &vel2,
2858 &wind_vector,
2859 (resolved_temp_c, resolved_press_hpa, base_ratio),
2860 );
2861
2862 let pos3 = position + vel2 * (dt * 0.5);
2864 let vel3 = velocity + acc2 * (dt * 0.5);
2865 let acc3 = self.calculate_acceleration(
2866 &pos3,
2867 &vel3,
2868 &wind_vector,
2869 (resolved_temp_c, resolved_press_hpa, base_ratio),
2870 );
2871
2872 let pos4 = position + vel3 * dt;
2874 let vel4 = velocity + acc3 * dt;
2875 let acc4 = self.calculate_acceleration(
2876 &pos4,
2877 &vel4,
2878 &wind_vector,
2879 (resolved_temp_c, resolved_press_hpa, base_ratio),
2880 );
2881
2882 position += (velocity + vel2 * 2.0 + vel3 * 2.0 + vel4) * (dt / 6.0);
2884 velocity += (acc1 + acc2 * 2.0 + acc3 * 2.0 + acc4) * (dt / 6.0);
2885 time += dt;
2886 self.validate_integration_state(&position, &velocity, time)?;
2887 }
2888
2889 let termination =
2890 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
2891
2892 self.annotate_drag_coefficients(&mut points, speed_of_sound);
2897
2898 let last_point = points.last().ok_or("No trajectory points generated")?;
2899
2900 let sampled_points = self.build_sampled_points(
2902 &points,
2903 max_height,
2904 transonic_distances,
2905 &mach_transitions,
2906 )?;
2907
2908 Ok(TrajectoryResult {
2909 max_range: last_point.position.x, max_height,
2911 time_of_flight: last_point.time,
2912 impact_velocity: last_point.velocity_magnitude,
2913 impact_energy: last_point.kinetic_energy,
2914 projectile_mass_kg: self.inputs.bullet_mass,
2915 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
2916 station_speed_of_sound_mps: speed_of_sound,
2917 termination,
2918 points,
2919 sampled_points,
2920 min_pitch_damping: if self.inputs.enable_pitch_damping {
2921 Some(min_pitch_damping)
2922 } else {
2923 None
2924 },
2925 transonic_mach,
2926 angular_state,
2927 max_yaw_angle: if self.inputs.enable_precession_nutation {
2928 Some(max_yaw_angle)
2929 } else {
2930 None
2931 },
2932 max_precession_angle: if self.inputs.enable_precession_nutation {
2933 Some(max_precession_angle)
2934 } else {
2935 None
2936 },
2937 aerodynamic_jump: aj_components,
2938 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
2939 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
2940 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
2941 })
2942 }
2943
2944 fn solve_rk45(&self) -> Result<TrajectoryResult, BallisticsError> {
2945 let mut time = 0.0;
2947 let mut position = self.initial_position();
2951
2952 let aj_components = self.aerodynamic_jump_components();
2958 let (launch_elev, launch_azim) = self.launch_angles_from(aj_components.as_ref());
2959 let horizontal_velocity = self.inputs.muzzle_velocity * launch_elev.cos();
2960 let mut velocity = Vector3::new(
2961 horizontal_velocity * launch_azim.cos(), self.inputs.muzzle_velocity * launch_elev.sin(), horizontal_velocity * launch_azim.sin(), );
2965
2966 let mut points = Vec::new();
2967 let mut max_height = position.y;
2968 let mut dt = 0.001; let (air_density, speed_of_sound, resolved_temp_c, resolved_press_hpa) =
2973 self.resolved_atmosphere();
2974 let base_ratio = air_density / 1.225;
2979 let wind_vector =
2984 crate::wind::wind_vector(self.wind.speed, self.wind.direction, self.wind.vertical_speed);
2985
2986 let mut transonic_distances: Vec<f64> = Vec::new();
2988 let mut mach_transitions = MachTransitionTracker::default();
2989
2990 let mut min_pitch_damping = f64::INFINITY;
2995 let mut transonic_mach: Option<f64> = None;
2996 let pitch_coeffs = PitchDampingCoefficients::from_bullet_type(
2997 self.inputs.bullet_model.as_deref().unwrap_or("default"),
2998 );
2999 let mut angular_state = if self.inputs.enable_precession_nutation {
3000 Some(AngularState {
3001 pitch_angle: 0.001,
3002 yaw_angle: 0.001,
3003 pitch_rate: 0.0,
3004 yaw_rate: 0.0,
3005 precession_angle: 0.0,
3006 nutation_phase: 0.0,
3007 })
3008 } else {
3009 None
3010 };
3011 let mut max_yaw_angle = 0.0;
3012 let mut max_precession_angle = 0.0;
3013
3014 while position.x < self.max_range
3015 && position.y > self.inputs.ground_threshold
3016 && time < TRAJECTORY_TIME_LIMIT_S
3017 {
3018 let velocity_magnitude = velocity.magnitude();
3020 let kinetic_energy = 0.5 * self.inputs.bullet_mass * velocity_magnitude.powi(2);
3021
3022 self.push_trajectory_point(
3023 &mut points,
3024 TrajectoryPoint {
3025 time,
3026 position,
3027 velocity_magnitude,
3028 kinetic_energy,
3029 drag_coefficient: None,
3030 },
3031 )?;
3032
3033 {
3036 let mach_here = if speed_of_sound > 0.0 {
3037 velocity_magnitude / speed_of_sound
3038 } else {
3039 0.0
3040 };
3041 mach_transitions.record_downward_crossings(
3042 mach_here,
3043 position.x,
3044 &mut transonic_distances,
3045 );
3046 }
3047
3048 if position.y > max_height {
3049 max_height = position.y;
3050 }
3051
3052 if self.inputs.enable_pitch_damping {
3055 let mach = velocity_magnitude / speed_of_sound;
3056 if transonic_mach.is_none() && mach < 1.2 && mach > 0.8 {
3057 transonic_mach = Some(mach);
3058 }
3059 let pitch_damping = calculate_pitch_damping_coefficient(mach, &pitch_coeffs);
3060 if pitch_damping < min_pitch_damping {
3061 min_pitch_damping = pitch_damping;
3062 }
3063 }
3064
3065 let accepted_step = self.adaptive_rk45_step(
3068 &position,
3069 &velocity,
3070 dt,
3071 &wind_vector,
3072 (resolved_temp_c, resolved_press_hpa, base_ratio),
3073 );
3074 debug_assert!(
3075 accepted_step.error <= RK45_TOLERANCE || accepted_step.used_dt <= RK45_MIN_DT
3076 );
3077
3078 if self.inputs.enable_precession_nutation {
3082 if let Some(ref mut state) = angular_state {
3083 let params = self.precession_nutation_params(
3084 velocity_magnitude,
3085 air_density,
3086 speed_of_sound,
3087 );
3088
3089 *state = calculate_combined_angular_motion(
3090 ¶ms,
3091 state,
3092 time,
3093 accepted_step.used_dt,
3094 0.001,
3095 );
3096
3097 if state.yaw_angle.abs() > max_yaw_angle {
3098 max_yaw_angle = state.yaw_angle.abs();
3099 }
3100 if state.precession_angle.abs() > max_precession_angle {
3101 max_precession_angle = state.precession_angle.abs();
3102 }
3103 }
3104 }
3105
3106 position = accepted_step.position;
3107 velocity = accepted_step.velocity;
3108 time += accepted_step.used_dt;
3109 self.validate_integration_state(&position, &velocity, time)?;
3110
3111 dt = accepted_step.next_dt;
3113 }
3114
3115 if points.is_empty() {
3117 return Err(BallisticsError::from("No trajectory points calculated"));
3118 }
3119
3120 let termination =
3122 self.append_terminal_endpoint(&mut points, position, velocity, time, &mut max_height)?;
3123
3124 self.annotate_drag_coefficients(&mut points, speed_of_sound);
3126
3127 let last_point = points.last().unwrap();
3128
3129 let sampled_points = self.build_sampled_points(
3131 &points,
3132 max_height,
3133 transonic_distances,
3134 &mach_transitions,
3135 )?;
3136
3137 Ok(TrajectoryResult {
3138 max_range: last_point.position.x, max_height,
3140 time_of_flight: last_point.time,
3141 impact_velocity: last_point.velocity_magnitude,
3142 impact_energy: last_point.kinetic_energy,
3143 projectile_mass_kg: self.inputs.bullet_mass,
3144 line_of_sight_height_m: self.inputs.muzzle_height + self.inputs.sight_height,
3145 station_speed_of_sound_mps: speed_of_sound,
3146 termination,
3147 points,
3148 sampled_points,
3149 min_pitch_damping: if self.inputs.enable_pitch_damping {
3150 Some(min_pitch_damping)
3151 } else {
3152 None
3153 },
3154 transonic_mach,
3155 angular_state,
3156 max_yaw_angle: if self.inputs.enable_precession_nutation {
3157 Some(max_yaw_angle)
3158 } else {
3159 None
3160 },
3161 max_precession_angle: if self.inputs.enable_precession_nutation {
3162 Some(max_precession_angle)
3163 } else {
3164 None
3165 },
3166 aerodynamic_jump: aj_components,
3167 mach_1_2_distance_m: mach_transitions.mach_1_2_distance_m,
3168 mach_1_0_distance_m: mach_transitions.mach_1_0_distance_m,
3169 mach_0_9_distance_m: mach_transitions.mach_0_9_distance_m,
3170 })
3171 }
3172
3173 fn adaptive_rk45_step(
3174 &self,
3175 position: &Vector3<f64>,
3176 velocity: &Vector3<f64>,
3177 initial_dt: f64,
3178 wind_vector: &Vector3<f64>,
3179 resolved_atmo: (f64, f64, f64),
3180 ) -> Rk45AcceptedStep {
3181 let mut trial_dt = initial_dt;
3182
3183 loop {
3184 let trial = self.rk45_step(
3185 position,
3186 velocity,
3187 trial_dt,
3188 wind_vector,
3189 RK45_TOLERANCE,
3190 resolved_atmo,
3191 );
3192 let next_dt = if trial.suggested_dt.is_finite() {
3197 (RK45_SAFETY_FACTOR * trial.suggested_dt).clamp(RK45_MIN_DT, RK45_MAX_DT)
3198 } else {
3199 RK45_MIN_DT
3200 };
3201
3202 if trial.error <= RK45_TOLERANCE || trial_dt <= RK45_MIN_DT {
3203 return Rk45AcceptedStep {
3204 position: trial.position,
3205 velocity: trial.velocity,
3206 used_dt: trial_dt,
3207 next_dt,
3208 error: trial.error,
3209 };
3210 }
3211
3212 trial_dt = next_dt;
3213 }
3214 }
3215
3216 fn rk45_step(
3217 &self,
3218 position: &Vector3<f64>,
3219 velocity: &Vector3<f64>,
3220 dt: f64,
3221 wind_vector: &Vector3<f64>,
3222 tolerance: f64,
3223 resolved_atmo: (f64, f64, f64), ) -> Rk45Trial {
3225 const A21: f64 = 1.0 / 5.0;
3227 const A31: f64 = 3.0 / 40.0;
3228 const A32: f64 = 9.0 / 40.0;
3229 const A41: f64 = 44.0 / 45.0;
3230 const A42: f64 = -56.0 / 15.0;
3231 const A43: f64 = 32.0 / 9.0;
3232 const A51: f64 = 19372.0 / 6561.0;
3233 const A52: f64 = -25360.0 / 2187.0;
3234 const A53: f64 = 64448.0 / 6561.0;
3235 const A54: f64 = -212.0 / 729.0;
3236 const A61: f64 = 9017.0 / 3168.0;
3237 const A62: f64 = -355.0 / 33.0;
3238 const A63: f64 = 46732.0 / 5247.0;
3239 const A64: f64 = 49.0 / 176.0;
3240 const A65: f64 = -5103.0 / 18656.0;
3241 const A71: f64 = 35.0 / 384.0;
3242 const A73: f64 = 500.0 / 1113.0;
3243 const A74: f64 = 125.0 / 192.0;
3244 const A75: f64 = -2187.0 / 6784.0;
3245 const A76: f64 = 11.0 / 84.0;
3246
3247 const B1: f64 = 35.0 / 384.0;
3249 const B3: f64 = 500.0 / 1113.0;
3250 const B4: f64 = 125.0 / 192.0;
3251 const B5: f64 = -2187.0 / 6784.0;
3252 const B6: f64 = 11.0 / 84.0;
3253
3254 const B1_ERR: f64 = 5179.0 / 57600.0;
3256 const B3_ERR: f64 = 7571.0 / 16695.0;
3257 const B4_ERR: f64 = 393.0 / 640.0;
3258 const B5_ERR: f64 = -92097.0 / 339200.0;
3259 const B6_ERR: f64 = 187.0 / 2100.0;
3260 const B7_ERR: f64 = 1.0 / 40.0;
3261
3262 let k1_v = self.calculate_acceleration(position, velocity, wind_vector, resolved_atmo);
3264 let k1_p = *velocity;
3265
3266 let p2 = position + dt * A21 * k1_p;
3267 let v2 = velocity + dt * A21 * k1_v;
3268 let k2_v = self.calculate_acceleration(&p2, &v2, wind_vector, resolved_atmo);
3269 let k2_p = v2;
3270
3271 let p3 = position + dt * (A31 * k1_p + A32 * k2_p);
3272 let v3 = velocity + dt * (A31 * k1_v + A32 * k2_v);
3273 let k3_v = self.calculate_acceleration(&p3, &v3, wind_vector, resolved_atmo);
3274 let k3_p = v3;
3275
3276 let p4 = position + dt * (A41 * k1_p + A42 * k2_p + A43 * k3_p);
3277 let v4 = velocity + dt * (A41 * k1_v + A42 * k2_v + A43 * k3_v);
3278 let k4_v = self.calculate_acceleration(&p4, &v4, wind_vector, resolved_atmo);
3279 let k4_p = v4;
3280
3281 let p5 = position + dt * (A51 * k1_p + A52 * k2_p + A53 * k3_p + A54 * k4_p);
3282 let v5 = velocity + dt * (A51 * k1_v + A52 * k2_v + A53 * k3_v + A54 * k4_v);
3283 let k5_v = self.calculate_acceleration(&p5, &v5, wind_vector, resolved_atmo);
3284 let k5_p = v5;
3285
3286 let p6 = position + dt * (A61 * k1_p + A62 * k2_p + A63 * k3_p + A64 * k4_p + A65 * k5_p);
3287 let v6 = velocity + dt * (A61 * k1_v + A62 * k2_v + A63 * k3_v + A64 * k4_v + A65 * k5_v);
3288 let k6_v = self.calculate_acceleration(&p6, &v6, wind_vector, resolved_atmo);
3289 let k6_p = v6;
3290
3291 let p7 = position + dt * (A71 * k1_p + A73 * k3_p + A74 * k4_p + A75 * k5_p + A76 * k6_p);
3292 let v7 = velocity + dt * (A71 * k1_v + A73 * k3_v + A74 * k4_v + A75 * k5_v + A76 * k6_v);
3293 let k7_v = self.calculate_acceleration(&p7, &v7, wind_vector, resolved_atmo);
3294 let k7_p = v7;
3295
3296 let new_pos = position + dt * (B1 * k1_p + B3 * k3_p + B4 * k4_p + B5 * k5_p + B6 * k6_p);
3298 let new_vel = velocity + dt * (B1 * k1_v + B3 * k3_v + B4 * k4_v + B5 * k5_v + B6 * k6_v);
3299
3300 let pos_err = position
3302 + dt * (B1_ERR * k1_p
3303 + B3_ERR * k3_p
3304 + B4_ERR * k4_p
3305 + B5_ERR * k5_p
3306 + B6_ERR * k6_p
3307 + B7_ERR * k7_p);
3308 let vel_err = velocity
3309 + dt * (B1_ERR * k1_v
3310 + B3_ERR * k3_v
3311 + B4_ERR * k4_v
3312 + B5_ERR * k5_v
3313 + B6_ERR * k6_v
3314 + B7_ERR * k7_v);
3315
3316 let error = cli_rk45_error_norm(position, velocity, &new_pos, &new_vel, &pos_err, &vel_err);
3318
3319 let dt_new = if error < tolerance {
3321 dt * (tolerance / error).powf(0.2).min(2.0)
3322 } else {
3323 dt * (tolerance / error).powf(0.25).max(0.1)
3324 };
3325
3326 Rk45Trial {
3327 position: new_pos,
3328 velocity: new_vel,
3329 suggested_dt: dt_new,
3330 error,
3331 }
3332 }
3333
3334 fn apply_cluster_bc_correction(&self, base_bc: f64, velocity_fps: f64) -> f64 {
3335 if let Some(ref cluster_bc) = self.cluster_bc {
3336 cluster_bc.apply_correction_for_drag_model(
3337 base_bc,
3338 self.inputs.caliber_inches,
3339 self.inputs.weight_grains,
3340 velocity_fps,
3341 self.inputs.bc_type,
3342 )
3343 } else {
3344 base_bc
3345 }
3346 }
3347
3348 fn calculate_acceleration(
3349 &self,
3350 position: &Vector3<f64>,
3351 velocity: &Vector3<f64>,
3352 wind_vector: &Vector3<f64>,
3353 resolved_atmo: (f64, f64, f64), ) -> Vector3<f64> {
3355 let actual_wind = if let Some(ref sock) = self.wind_sock {
3361 sock.vector_for_range_stateless(position.x)
3362 } else if self.inputs.enable_wind_shear {
3363 self.get_wind_at_altitude(position.y)
3364 } else {
3365 *wind_vector
3366 };
3367 let actual_wind =
3368 crate::derivatives::level_vector_to_shot_frame(actual_wind, self.inputs.shooting_angle);
3369
3370 let relative_velocity = velocity - actual_wind;
3371 let velocity_magnitude = relative_velocity.magnitude();
3372
3373 if velocity_magnitude < 0.001 {
3374 return self.gravity_acceleration();
3375 }
3376
3377 let (base_temp_c, base_press_hpa, station_ratio) = resolved_atmo;
3388
3389 let (drag_base_temp_c, drag_base_press_hpa, drag_base_ratio, drag_humidity_percent) =
3397 if let Some(ref sock) = self.atmo_sock {
3398 let (zone_temp_c, zone_press_hpa, zone_humidity) = sock.atmo_for_range(position.x);
3399 let zone_base_ratio = crate::atmosphere::calculate_air_density_cimp(
3400 zone_temp_c,
3401 zone_press_hpa,
3402 zone_humidity,
3403 ) / 1.225;
3404 (zone_temp_c, zone_press_hpa, zone_base_ratio, zone_humidity)
3405 } else {
3406 (
3407 base_temp_c,
3408 base_press_hpa,
3409 station_ratio,
3410 self.atmosphere.humidity,
3411 )
3412 };
3413 let local_alt = crate::atmosphere::shot_frame_altitude(
3414 self.atmosphere.altitude,
3415 position.x,
3416 position.y,
3417 self.inputs.shooting_angle,
3418 );
3419 let (air_density, speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
3420 local_alt,
3421 self.atmosphere.altitude,
3422 drag_base_temp_c,
3423 drag_base_press_hpa,
3424 drag_base_ratio,
3425 drag_humidity_percent,
3426 );
3427
3428 let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
3432
3433 let velocity_fps = velocity_magnitude * 3.28084;
3435
3436 let cd_to_retard = crate::constants::CD_TO_RETARD;
3441 let standard_factor = cd * cd_to_retard;
3442 let density_scale = air_density / 1.225; let a_drag_ft_s2 =
3446 (velocity_fps * velocity_fps) * standard_factor * density_scale / retard_denom;
3447 let a_drag_m_s2 = a_drag_ft_s2 * 0.3048; let drag_acceleration = -a_drag_m_s2 * (relative_velocity / velocity_magnitude);
3451
3452 let mut accel = drag_acceleration + self.gravity_acceleration();
3455
3456 if self.inputs.enable_coriolis {
3459 if let Some(lat_deg) = self.inputs.latitude {
3460 let omega_earth = 7.2921159e-5_f64; let lat = lat_deg.to_radians();
3462 let az = self.inputs.shot_azimuth; let omega = Vector3::new(
3469 omega_earth * lat.cos() * az.cos(), omega_earth * lat.sin(), -omega_earth * lat.cos() * az.sin(), );
3473 let omega = crate::derivatives::level_vector_to_shot_frame(
3474 omega,
3475 self.inputs.shooting_angle,
3476 );
3477 accel += -2.0 * omega.cross(velocity);
3482 }
3483 }
3484
3485 if self.inputs.enable_magnus
3492 && !self.inputs.use_enhanced_spin_drift
3493 && self.inputs.bullet_diameter > 0.0
3494 && self.inputs.twist_rate > 0.0
3495 {
3496 let diameter_m = self.inputs.bullet_diameter;
3497 let (spin_rad_s, spin_param) = crate::spin_drift::calculate_magnus_spin_state(
3498 self.inputs.muzzle_velocity,
3499 velocity_magnitude,
3500 self.inputs.twist_rate,
3501 diameter_m,
3502 );
3503 let mach = velocity_magnitude / speed_of_sound;
3505
3506 let d_in = self.inputs.bullet_diameter / 0.0254;
3508 let m_gr = self.inputs.bullet_mass / crate::constants::GRAINS_TO_KG;
3509 let l_in = if self.inputs.bullet_length > 0.0 {
3510 self.inputs.bullet_length / 0.0254
3511 } else {
3512 let est_m = crate::stability::estimate_bullet_length_m(
3514 self.inputs.bullet_diameter,
3515 self.inputs.bullet_mass,
3516 );
3517 if est_m > 0.0 {
3518 est_m / 0.0254
3519 } else {
3520 4.5 * d_in
3521 }
3522 };
3523 let sg = crate::spin_drift::calculate_dynamic_stability(
3527 m_gr,
3528 velocity_magnitude,
3529 spin_rad_s,
3530 d_in,
3531 l_in,
3532 air_density,
3533 );
3534
3535 let (yaw_rad, _) = crate::spin_drift::calculate_yaw_of_repose(
3537 sg,
3538 velocity_magnitude,
3539 spin_rad_s,
3540 0.0, 0.0, air_density,
3543 d_in,
3544 l_in,
3545 m_gr,
3546 mach,
3547 "match",
3548 false,
3549 );
3550
3551 let c_np = crate::derivatives::calculate_magnus_moment_coefficient(mach);
3553 let area = std::f64::consts::PI * (diameter_m / 2.0).powi(2);
3554 let magnus_force = 0.5
3555 * air_density
3556 * velocity_magnitude.powi(2)
3557 * area
3558 * c_np
3559 * spin_param
3560 * yaw_rad.sin();
3561
3562 if magnus_force.abs() > 1e-12 {
3566 if let Some(dir) = crate::derivatives::yaw_of_repose_magnus_direction(
3567 relative_velocity,
3568 self.gravity_acceleration(),
3569 self.inputs.is_twist_right,
3570 ) {
3571 accel += (magnus_force / self.inputs.bullet_mass) * dir;
3572 }
3573 }
3574 }
3575
3576 accel
3577 }
3578
3579 fn drag_terms(&self, velocity_magnitude: f64, speed_of_sound: f64) -> (f64, f64) {
3591 let cd = self.calculate_drag_coefficient(velocity_magnitude, speed_of_sound);
3592
3593 let velocity_fps = velocity_magnitude * 3.28084;
3594
3595 let (base_bc, bc_from_segments) = if let Some(segments) = self
3600 .inputs
3601 .bc_segments_data
3602 .as_ref()
3603 .filter(|segments| self.inputs.use_bc_segments && !segments.is_empty())
3604 {
3605 (
3607 crate::bc_estimation::velocity_segment_bc(
3608 velocity_fps,
3609 segments,
3610 self.inputs.bc_value,
3611 ),
3612 true,
3613 )
3614 } else if let Some(segments) = self
3615 .inputs
3616 .bc_segments
3617 .as_ref()
3618 .filter(|segments| !segments.is_empty())
3619 {
3620 (
3621 crate::derivatives::interpolated_bc(
3622 velocity_magnitude / speed_of_sound,
3623 segments,
3624 Some(&self.inputs),
3625 ),
3626 true,
3627 )
3628 } else {
3629 (self.inputs.bc_value, false)
3630 };
3631
3632 let effective_bc = if bc_from_segments {
3637 base_bc
3638 } else {
3639 self.apply_cluster_bc_correction(base_bc, velocity_fps)
3640 };
3641 let effective_bc = effective_bc.max(1e-6);
3644
3645 let retard_denom = if self.inputs.custom_drag_table.is_some() {
3650 self.inputs.custom_drag_denominator(effective_bc)
3651 } else {
3652 effective_bc
3653 };
3654
3655 (cd, retard_denom)
3656 }
3657
3658 pub fn effective_drag_coefficient(
3679 &self,
3680 velocity_magnitude: f64,
3681 speed_of_sound: f64,
3682 ) -> Option<f64> {
3683 if !velocity_magnitude.is_finite() || speed_of_sound <= 1e-9 {
3684 return None;
3685 }
3686 let sectional_density = self.inputs.sectional_density_lb_in2()?;
3687 let (cd, retard_denom) = self.drag_terms(velocity_magnitude, speed_of_sound);
3688 if retard_denom <= 0.0 {
3689 return None;
3690 }
3691 let effective = cd * sectional_density / retard_denom;
3692 effective.is_finite().then_some(effective)
3693 }
3694
3695 fn annotate_drag_coefficients(&self, points: &mut [TrajectoryPoint], speed_of_sound: f64) {
3705 for point in points.iter_mut() {
3706 point.drag_coefficient =
3707 self.effective_drag_coefficient(point.velocity_magnitude, speed_of_sound);
3708 }
3709 }
3710
3711 fn calculate_drag_coefficient(&self, velocity: f64, speed_of_sound: f64) -> f64 {
3712 let mach = velocity / speed_of_sound;
3713
3714 if let Some(ref table) = self.inputs.custom_drag_table {
3718 return table.interpolate(mach) * self.inputs.cd_scale;
3723 }
3724
3725 crate::drag::get_drag_coefficient(mach, &self.inputs.bc_type)
3728 }
3729}
3730
3731#[derive(Debug, Clone)]
3733pub struct MonteCarloParams {
3734 pub num_simulations: usize,
3735 pub velocity_std_dev: f64,
3736 pub angle_std_dev: f64,
3737 pub bc_std_dev: f64,
3738 pub wind_speed_std_dev: f64,
3739 pub target_distance: Option<f64>,
3740 pub base_wind_speed: f64,
3741 pub base_wind_direction: f64,
3742 pub azimuth_std_dev: f64, }
3744
3745impl Default for MonteCarloParams {
3746 fn default() -> Self {
3747 Self {
3748 num_simulations: 1000,
3749 velocity_std_dev: 1.0,
3750 angle_std_dev: 0.001,
3751 bc_std_dev: 0.01,
3752 wind_speed_std_dev: 1.0,
3753 target_distance: None,
3754 base_wind_speed: 0.0,
3755 base_wind_direction: 0.0,
3756 azimuth_std_dev: 0.001, }
3758 }
3759}
3760
3761#[derive(Debug, Clone)]
3763pub struct MonteCarloResults {
3764 pub ranges: Vec<f64>,
3765 pub impact_velocities: Vec<f64>,
3766 pub impact_positions: Vec<Vector3<f64>>,
3772}
3773
3774pub const DEFAULT_HIT_RADIUS_M: f64 = 0.3;
3777
3778pub const TARGET_NOT_REACHED_SENTINEL_M: f64 = -1.0e9;
3784
3785impl MonteCarloResults {
3786 pub fn position_reached_target(position: &Vector3<f64>) -> bool {
3788 position.iter().all(|component| component.is_finite())
3789 && position.y != TARGET_NOT_REACHED_SENTINEL_M
3790 }
3791
3792 pub fn target_arrival_count(&self) -> usize {
3794 self.impact_positions
3795 .iter()
3796 .filter(|position| Self::position_reached_target(position))
3797 .count()
3798 }
3799
3800 pub fn target_shortfall_fraction(&self) -> f64 {
3803 if self.impact_positions.is_empty() {
3804 return 0.0;
3805 }
3806 (self.impact_positions.len() - self.target_arrival_count()) as f64
3807 / self.impact_positions.len() as f64
3808 }
3809
3810 pub fn target_plane_cep(&self) -> Option<f64> {
3816 let mut radial_misses: Vec<f64> = self
3817 .impact_positions
3818 .iter()
3819 .filter(|position| Self::position_reached_target(position))
3820 .map(Vector3::norm)
3821 .filter(|miss| miss.is_finite())
3822 .collect();
3823 radial_misses.sort_by(f64::total_cmp);
3824 if radial_misses.is_empty() {
3825 None
3826 } else {
3827 Some(radial_misses[radial_misses.len() / 2])
3828 }
3829 }
3830
3831 pub fn hit_probability(&self, hit_radius_m: f64) -> f64 {
3840 if self.impact_positions.is_empty() {
3841 return 0.0;
3842 }
3843 let hits = self
3844 .impact_positions
3845 .iter()
3846 .filter(|position| {
3847 Self::position_reached_target(position) && position.norm() < hit_radius_m
3848 })
3849 .count();
3850 hits as f64 / self.impact_positions.len() as f64
3851 }
3852
3853 pub fn rect_hit_probability(&self, width_m: f64, height_m: f64) -> f64 {
3865 let dimensions_invalid = width_m.is_nan()
3866 || width_m <= 0.0
3867 || height_m.is_nan()
3868 || height_m <= 0.0;
3869 if self.impact_positions.is_empty() || dimensions_invalid {
3870 return 0.0;
3871 }
3872 let half_width = width_m / 2.0;
3873 let half_height = height_m / 2.0;
3874 let hits = self
3875 .impact_positions
3876 .iter()
3877 .filter(|position| {
3878 Self::position_reached_target(position)
3879 && position.z.abs() <= half_width
3880 && position.y.abs() <= half_height
3881 })
3882 .count();
3883 hits as f64 / self.impact_positions.len() as f64
3884 }
3885}
3886
3887fn wind_from_signed_speed_sample(
3888 signed_speed: f64,
3889 sampled_direction: f64,
3890 vertical_speed: f64,
3891) -> WindConditions {
3892 if signed_speed < 0.0 {
3897 WindConditions {
3898 speed: -signed_speed,
3899 direction: sampled_direction + std::f64::consts::PI,
3900 vertical_speed,
3901 }
3902 } else {
3903 WindConditions {
3904 speed: signed_speed,
3905 direction: sampled_direction,
3906 vertical_speed,
3907 }
3908 }
3909}
3910
3911struct MonteCarloWindSampler {
3912 speed: rand_distr::Normal<f64>,
3913 direction: rand_distr::Normal<f64>,
3914 vertical_speed: f64,
3916}
3917
3918impl MonteCarloWindSampler {
3919 fn new(
3920 base_wind: &WindConditions,
3921 wind_speed_std_dev: f64,
3922 wind_direction_std_dev: f64,
3923 ) -> Result<Self, BallisticsError> {
3924 use rand_distr::Normal;
3925
3926 if !wind_direction_std_dev.is_finite() || wind_direction_std_dev < 0.0 {
3927 return Err("Wind direction standard deviation must be finite and non-negative".into());
3928 }
3929
3930 let speed = Normal::new(base_wind.speed, wind_speed_std_dev)
3931 .map_err(|e| format!("Invalid wind speed distribution: {e}"))?;
3932 let direction = Normal::new(base_wind.direction, wind_direction_std_dev)
3933 .map_err(|e| format!("Invalid wind direction distribution: {e}"))?;
3934 Ok(Self { speed, direction, vertical_speed: base_wind.vertical_speed })
3935 }
3936
3937 fn sample<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> WindConditions {
3938 use rand_distr::Distribution;
3939
3940 wind_from_signed_speed_sample(
3941 self.speed.sample(rng),
3942 self.direction.sample(rng),
3943 self.vertical_speed,
3944 )
3945 }
3946}
3947
3948pub fn run_monte_carlo(
3950 base_inputs: BallisticInputs,
3951 params: MonteCarloParams,
3952) -> Result<MonteCarloResults, BallisticsError> {
3953 run_monte_carlo_with_direction_std_dev(base_inputs, params, 0.0)
3954}
3955
3956pub fn run_monte_carlo_with_direction_std_dev(
3961 base_inputs: BallisticInputs,
3962 params: MonteCarloParams,
3963 wind_direction_std_dev: f64,
3964) -> Result<MonteCarloResults, BallisticsError> {
3965 let base_wind = WindConditions {
3966 speed: params.base_wind_speed,
3967 direction: params.base_wind_direction,
3968 vertical_speed: 0.0,
3969 };
3970 run_monte_carlo_with_wind_and_direction_std_dev(
3971 base_inputs,
3972 base_wind,
3973 params,
3974 wind_direction_std_dev,
3975 )
3976}
3977
3978pub fn run_monte_carlo_with_wind(
3980 base_inputs: BallisticInputs,
3981 base_wind: WindConditions,
3982 params: MonteCarloParams,
3983) -> Result<MonteCarloResults, BallisticsError> {
3984 run_monte_carlo_with_wind_and_direction_std_dev(base_inputs, base_wind, params, 0.0)
3985}
3986
3987pub fn run_monte_carlo_with_wind_and_direction_std_dev(
3992 base_inputs: BallisticInputs,
3993 base_wind: WindConditions,
3994 params: MonteCarloParams,
3995 wind_direction_std_dev: f64,
3996) -> Result<MonteCarloResults, BallisticsError> {
3997 let mut rng = rand::rng();
3998 run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
3999 base_inputs,
4000 base_wind,
4001 params,
4002 wind_direction_std_dev,
4003 &mut rng,
4004 )
4005}
4006
4007pub fn run_monte_carlo_with_wind_and_direction_std_dev_seeded(
4014 base_inputs: BallisticInputs,
4015 base_wind: WindConditions,
4016 params: MonteCarloParams,
4017 wind_direction_std_dev: f64,
4018 seed: u64,
4019) -> Result<MonteCarloResults, BallisticsError> {
4020 use rand::{rngs::StdRng, SeedableRng};
4021 let mut rng = StdRng::seed_from_u64(seed);
4022 run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
4023 base_inputs,
4024 base_wind,
4025 params,
4026 wind_direction_std_dev,
4027 &mut rng,
4028 )
4029}
4030
4031fn run_monte_carlo_with_wind_and_direction_std_dev_using_rng<R: rand::Rng + ?Sized>(
4032 base_inputs: BallisticInputs,
4033 base_wind: WindConditions,
4034 params: MonteCarloParams,
4035 wind_direction_std_dev: f64,
4036 rng: &mut R,
4037) -> Result<MonteCarloResults, BallisticsError> {
4038 use rand_distr::{Distribution, Normal};
4039
4040 let mut ranges = Vec::new();
4041 let mut impact_velocities = Vec::new();
4042 let mut impact_positions = Vec::new();
4043
4044 let atmosphere = AtmosphericConditions {
4045 temperature: base_inputs.temperature,
4046 pressure: base_inputs.pressure,
4047 humidity: base_inputs.humidity_percent(),
4048 altitude: base_inputs.altitude,
4049 };
4050 let target_hint = params
4051 .target_distance
4052 .unwrap_or(base_inputs.target_distance);
4053 let solver_max_range = target_hint.max(1000.0) * 2.0;
4054
4055 let mut baseline_solver =
4057 TrajectorySolver::new(base_inputs.clone(), base_wind.clone(), atmosphere.clone());
4058 baseline_solver.set_max_range(solver_max_range);
4059 let baseline_result = baseline_solver.solve()?;
4060
4061 let target_distance = params.target_distance.unwrap_or(baseline_result.max_range);
4063
4064 let baseline_at_target = baseline_result
4066 .position_at_range(target_distance)
4067 .ok_or("Could not interpolate baseline at target distance")?;
4068
4069 let velocity_delta_dist = Normal::new(0.0, params.velocity_std_dev)
4074 .map_err(|e| format!("Invalid velocity distribution: {}", e))?;
4075 let angle_dist = Normal::new(base_inputs.muzzle_angle, params.angle_std_dev)
4076 .map_err(|e| format!("Invalid angle distribution: {}", e))?;
4077 let bc_dist = Normal::new(base_inputs.bc_value, params.bc_std_dev)
4078 .map_err(|e| format!("Invalid BC distribution: {}", e))?;
4079 let wind_sampler = MonteCarloWindSampler::new(
4082 &base_wind,
4083 params.wind_speed_std_dev,
4084 wind_direction_std_dev,
4085 )?;
4086 let azimuth_dist = Normal::new(base_inputs.azimuth_angle, params.azimuth_std_dev)
4087 .map_err(|e| format!("Invalid azimuth distribution: {}", e))?;
4088
4089 for _ in 0..params.num_simulations {
4090 let mut inputs = base_inputs.clone();
4092 let muzzle_velocity_delta = velocity_delta_dist.sample(&mut *rng);
4093 inputs.muzzle_angle = angle_dist.sample(&mut *rng);
4094 inputs.bc_value = bc_dist.sample(&mut *rng).max(0.01);
4095 inputs.azimuth_angle = azimuth_dist.sample(&mut *rng); let wind = wind_sampler.sample(&mut *rng);
4099
4100 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere.clone());
4102 solver.inputs.muzzle_velocity =
4103 (solver.inputs.muzzle_velocity + muzzle_velocity_delta).max(0.0);
4104 solver.set_max_range(solver_max_range);
4105 match solver.solve() {
4106 Ok(result) => {
4107 let deviation = if result.max_range < target_distance {
4113 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)
4116 } else {
4117 let pos_at_target = match result.position_at_range(target_distance) {
4118 Some(p) => p,
4119 None => continue, };
4121 Vector3::new(
4126 0.0,
4127 pos_at_target.y - baseline_at_target.y,
4128 pos_at_target.z - baseline_at_target.z,
4129 )
4130 };
4131
4132 ranges.push(result.max_range);
4133 impact_velocities.push(result.impact_velocity);
4134 impact_positions.push(deviation);
4135 }
4136 Err(_) => {
4137 continue;
4139 }
4140 }
4141 }
4142
4143 if ranges.is_empty() {
4144 return Err("No successful simulations".into());
4145 }
4146
4147 Ok(MonteCarloResults {
4148 ranges,
4149 impact_velocities,
4150 impact_positions,
4151 })
4152}
4153
4154pub fn calculate_zero_angle(
4156 inputs: BallisticInputs,
4157 target_distance: f64,
4158 target_height: f64,
4159) -> Result<f64, BallisticsError> {
4160 calculate_zero_angle_with_conditions(
4161 inputs,
4162 target_distance,
4163 target_height,
4164 WindConditions::default(),
4165 AtmosphericConditions::default(),
4166 )
4167}
4168
4169pub fn calculate_zero_angle_with_conditions(
4170 inputs: BallisticInputs,
4171 target_distance: f64,
4172 target_height: f64,
4173 wind: WindConditions,
4174 atmosphere: AtmosphericConditions,
4175) -> Result<f64, BallisticsError> {
4176 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
4177 solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
4178}
4179
4180pub fn calculate_zero_angle_with_resolved_conditions(
4186 inputs: BallisticInputs,
4187 target_distance: f64,
4188 target_height: f64,
4189 wind: WindConditions,
4190 atmosphere: AtmosphericConditions,
4191) -> Result<f64, BallisticsError> {
4192 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
4193 solver.calculate_and_set_zero_angle(target_distance, target_height, ZeroTargetFrame::SightLine)
4194}
4195
4196pub const ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M: f64 = 2000.0;
4204
4205pub fn calculate_zero_range_from_angle_with_conditions(
4220 inputs: BallisticInputs,
4221 zero_angle_rad: f64,
4222 target_height: f64,
4223 wind: WindConditions,
4224 atmosphere: AtmosphericConditions,
4225) -> Result<ZeroCrossings, BallisticsError> {
4226 let mut solver = TrajectorySolver::new(inputs, wind, atmosphere);
4227 solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
4228 solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
4229}
4230
4231pub fn calculate_zero_range_from_angle_with_resolved_conditions(
4237 inputs: BallisticInputs,
4238 zero_angle_rad: f64,
4239 target_height: f64,
4240 wind: WindConditions,
4241 atmosphere: AtmosphericConditions,
4242) -> Result<ZeroCrossings, BallisticsError> {
4243 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(inputs, wind, atmosphere);
4244 solver.set_max_range(ZERO_RANGE_FROM_ANGLE_MAX_RANGE_M);
4245 solver.find_zero_range(zero_angle_rad, target_height, ZeroTargetFrame::SightLine)
4246}
4247
4248#[derive(Debug, Clone, Copy, PartialEq, Eq)]
4250pub enum BcFitMode {
4251 Drop,
4253 Velocity,
4256}
4257
4258#[derive(Debug, Clone, Copy)]
4260pub struct BcEstimate {
4261 pub bc: f64,
4263 pub rms_error: f64,
4265 pub drag_model: DragModel,
4267 pub mode: BcFitMode,
4269 pub at_bound: bool,
4273}
4274
4275fn fit_value_at(
4283 points: &[TrajectoryPoint],
4284 target_dist: f64,
4285 mode: BcFitMode,
4286 drop_offset: f64,
4287) -> Option<f64> {
4288 let val = |p: &TrajectoryPoint| match mode {
4289 BcFitMode::Drop => drop_offset - p.position.y,
4290 BcFitMode::Velocity => p.velocity_magnitude,
4291 };
4292 for i in 0..points.len() {
4293 if points[i].position.x >= target_dist {
4294 if i == 0 {
4295 return Some(val(&points[0]));
4296 }
4297 let p1 = &points[i - 1];
4298 let p2 = &points[i];
4299 let dx = p2.position.x - p1.position.x;
4300 if dx.abs() < 1e-9 {
4301 return Some(val(p2));
4302 }
4303 let t = (target_dist - p1.position.x) / dx;
4304 return Some(val(p1) + t * (val(p2) - val(p1)));
4305 }
4306 }
4307 None
4308}
4309
4310fn fit_residual_sse(
4311 trajectory: &[TrajectoryPoint],
4312 observations: &[(f64, f64)],
4313 mode: BcFitMode,
4314 drop_offset: f64,
4315) -> Option<f64> {
4316 if observations.is_empty() {
4317 return None;
4318 }
4319 let mut total = 0.0;
4320 for (target_dist, target_val) in observations {
4321 let value = fit_value_at(trajectory, *target_dist, mode, drop_offset)?;
4324 let error = value - target_val;
4325 total += error * error;
4326 }
4327 Some(total)
4328}
4329
4330#[allow(clippy::too_many_arguments)] pub fn estimate_bc_fit(
4348 velocity: f64,
4349 mass: f64,
4350 diameter: f64,
4351 points: &[(f64, f64)],
4352 drag_model: DragModel,
4353 mode: BcFitMode,
4354 atmosphere: AtmosphericConditions,
4355 zero_range: Option<f64>,
4356 sight_height: f64,
4357) -> Result<BcEstimate, BallisticsError> {
4358 if points.is_empty() {
4359 return Err(BallisticsError::from(
4360 "No data points provided for BC estimation.".to_string(),
4361 ));
4362 }
4363 let max_dist = points.iter().map(|(d, _)| *d).fold(0.0_f64, f64::max);
4364 let drop_offset = if zero_range.is_some() { sight_height } else { 0.0 };
4367
4368 let sse = |bc_value: f64| -> Option<f64> {
4370 let mut inputs = BallisticInputs {
4371 muzzle_velocity: velocity,
4372 bc_value,
4373 bc_type: drag_model,
4374 bullet_mass: mass,
4375 bullet_diameter: diameter,
4376 sight_height,
4377 ..Default::default()
4378 };
4379 if let Some(zr) = zero_range {
4382 let za = calculate_zero_angle_with_conditions(
4388 inputs.clone(),
4389 zr,
4390 sight_height,
4391 WindConditions::default(),
4392 atmosphere.clone(),
4393 )
4394 .ok()?;
4395 inputs.muzzle_angle = za;
4396 }
4397 let mut solver =
4398 TrajectorySolver::new(inputs, WindConditions::default(), atmosphere.clone());
4399 solver.set_max_range(max_dist * 1.5);
4400 let result = solver.solve().ok()?;
4401 fit_residual_sse(&result.points, points, mode, drop_offset)
4402 };
4403
4404 let (bc_min, bc_max) = match drag_model {
4408 DragModel::G7 => (0.05, 0.70),
4409 _ => (0.10, 1.20),
4410 };
4411
4412 let mut best_bc = f64::NAN;
4414 let mut best_sse = f64::MAX;
4415 let mut bc = bc_min;
4416 while bc <= bc_max + 1e-9 {
4417 if let Some(s) = sse(bc) {
4418 if s < best_sse {
4419 best_sse = s;
4420 best_bc = bc;
4421 }
4422 }
4423 bc += 0.01;
4424 }
4425 if !best_bc.is_finite() {
4426 return Err(BallisticsError::from(
4427 "Unable to estimate BC from provided data. Check that the values and units are correct."
4428 .to_string(),
4429 ));
4430 }
4431
4432 let lo = (best_bc - 0.01).max(bc_min);
4434 let hi = (best_bc + 0.01).min(bc_max);
4435 let mut bc = lo;
4436 while bc <= hi + 1e-9 {
4437 if let Some(s) = sse(bc) {
4438 if s < best_sse {
4439 best_sse = s;
4440 best_bc = bc;
4441 }
4442 }
4443 bc += 0.001;
4444 }
4445
4446 let at_bound = best_bc <= bc_min + 0.011 || best_bc >= bc_max - 0.011;
4449 let rms_error = (best_sse / points.len() as f64).sqrt();
4452 Ok(BcEstimate {
4453 bc: best_bc,
4454 rms_error,
4455 drag_model,
4456 mode,
4457 at_bound,
4458 })
4459}
4460
4461pub fn estimate_bc_from_trajectory(
4464 velocity: f64,
4465 mass: f64,
4466 diameter: f64,
4467 points: &[(f64, f64)], ) -> Result<f64, BallisticsError> {
4469 estimate_bc_fit(
4470 velocity,
4471 mass,
4472 diameter,
4473 points,
4474 DragModel::G1,
4475 BcFitMode::Drop,
4476 AtmosphericConditions::default(),
4477 None,
4478 0.05,
4479 )
4480 .map(|e| e.bc)
4481}
4482
4483use rand;
4485use rand_distr;
4486
4487#[cfg(test)]
4488mod mba737_powder_resolution_tests {
4489 use super::*;
4490
4491 #[test]
4492 fn linear_model_cold_powder_subtracts() {
4493 let v = resolve_powder_adjusted_velocity(823.0, 11.1, true, 0.5486, 21.1, None, None);
4495 assert!((v - (823.0 + 0.5486 * (11.1 - 21.1))).abs() < 1e-12);
4496 assert!(v < 823.0);
4497 }
4498
4499 #[test]
4500 fn linear_model_hot_powder_adds() {
4501 let v = resolve_powder_adjusted_velocity(823.0, 31.1, true, 0.5486, 21.1, None, None);
4502 assert!((v - (823.0 + 0.5486 * 10.0)).abs() < 1e-12);
4503 }
4504
4505 #[test]
4506 fn disabled_flag_is_passthrough() {
4507 let v = resolve_powder_adjusted_velocity(823.0, 40.0, false, 0.5486, 21.1, None, None);
4508 assert_eq!(v, 823.0);
4509 }
4510
4511 #[test]
4512 fn curve_overrides_linear_and_interpolates_at_powder_temp() {
4513 let curve = [(4.4, 798.6), (21.1, 823.0), (37.8, 841.2)];
4514 let v = resolve_powder_adjusted_velocity(823.0, 30.0, true, 99.0, 21.1, Some(&curve), Some(4.4));
4516 assert!((v - 798.6).abs() < 1e-9);
4517 }
4518
4519 #[test]
4520 fn curve_falls_back_to_ambient_and_clamps() {
4521 let curve = [(4.4, 798.6), (37.8, 841.2)];
4522 let v = resolve_powder_adjusted_velocity(823.0, -40.0, true, 1.0, 21.1, Some(&curve), None);
4524 assert!((v - 798.6).abs() < 1e-9);
4525 let v_hot = resolve_powder_adjusted_velocity(823.0, 60.0, true, 1.0, 21.1, Some(&curve), None);
4526 assert!((v_hot - 841.2).abs() < 1e-9);
4527 }
4528
4529 #[test]
4530 fn empty_curve_suppresses_linear_fallback() {
4531 let v = resolve_powder_adjusted_velocity(823.0, 40.0, true, 0.5486, 21.1, Some(&[]), None);
4533 assert_eq!(v, 823.0);
4534 }
4535
4536 #[test]
4537 fn sweep_huge_range_errors_instead_of_overflowing() {
4538 assert!(parse_powder_sweep("0:1e20:1").is_err());
4541 assert!(parse_powder_sweep("0:1e308:1e-3").is_err());
4542 }
4543
4544 #[test]
4545 fn sweep_fractional_step_keeps_end_row() {
4546 let rows = parse_powder_sweep("0:0.3:0.1").unwrap();
4548 assert_eq!(rows.len(), 4);
4549 assert!((rows[3] - 0.3).abs() < 1e-9);
4550 }
4551
4552 #[test]
4553 fn solver_and_helper_agree_on_linear_model() {
4554 let inputs = BallisticInputs {
4556 use_powder_sensitivity: true,
4557 powder_temp_sensitivity: 0.5486,
4558 powder_temp: 21.1,
4559 temperature: 4.4,
4560 ..Default::default()
4561 };
4562 let expected = resolve_powder_adjusted_velocity(
4563 inputs.muzzle_velocity,
4564 inputs.temperature,
4565 true,
4566 0.5486,
4567 21.1,
4568 None,
4569 None,
4570 );
4571 let solver = TrajectorySolver::new(
4572 inputs,
4573 WindConditions::default(),
4574 AtmosphericConditions::default(),
4575 );
4576 assert!((solver.inputs.muzzle_velocity - expected).abs() < 1e-12);
4577 }
4578}
4579
4580#[cfg(test)]
4581mod mba1302_solver_seam_tests {
4582 use super::*;
4583 use crate::wind::WindSegment;
4584
4585 #[test]
4586 fn authoritative_station_atmosphere_preserves_explicit_standard_values_at_altitude() {
4587 let atmosphere = AtmosphericConditions {
4588 temperature: 15.0,
4589 pressure: 1013.25,
4590 humidity: 50.0,
4591 altitude: 2_000.0,
4592 };
4593 let legacy = TrajectorySolver::new(
4594 BallisticInputs::default(),
4595 WindConditions::default(),
4596 atmosphere.clone(),
4597 );
4598 let authoritative = TrajectorySolver::new_with_resolved_station_atmosphere(
4599 BallisticInputs::default(),
4600 WindConditions::default(),
4601 atmosphere,
4602 );
4603
4604 let (legacy_density, _, legacy_temp_c, legacy_pressure_hpa) = legacy.resolved_atmosphere();
4605 let (authoritative_density, _, authoritative_temp_c, authoritative_pressure_hpa) =
4606 authoritative.resolved_atmosphere();
4607 let (icao_temp_k, icao_pressure_pa) =
4608 crate::atmosphere::calculate_icao_standard_atmosphere(2_000.0);
4609 let (expected_authoritative_density, _) =
4610 crate::atmosphere::calculate_atmosphere(2_000.0, Some(15.0), Some(1013.25), 50.0);
4611
4612 assert!((legacy_temp_c - (icao_temp_k - 273.15)).abs() < 1e-12);
4613 assert!((legacy_pressure_hpa - icao_pressure_pa / 100.0).abs() < 1e-12);
4614 assert_eq!(authoritative_temp_c.to_bits(), 15.0_f64.to_bits());
4615 assert_eq!(authoritative_pressure_hpa.to_bits(), 1013.25_f64.to_bits());
4616 assert_eq!(
4617 authoritative_density.to_bits(),
4618 expected_authoritative_density.to_bits()
4619 );
4620 assert!(
4621 (authoritative_density - legacy_density).abs() > 0.1,
4622 "explicit standard values at altitude must differ from ICAO-at-altitude: explicit={authoritative_density}, ICAO={legacy_density}"
4623 );
4624 }
4625
4626 #[test]
4635 fn precomputed_absolute_resolution_via_authoritative_matches_legacy_new() {
4636 for (temperature, pressure, altitude) in [
4637 (15.0, 1013.25, 0.0), (15.0, 1013.25, 2000.0), (-5.0, 850.0, 2000.0), (22.0, 950.0, 500.0),
4641 ] {
4642 let atmosphere = AtmosphericConditions {
4643 temperature,
4644 pressure,
4645 humidity: 50.0,
4646 altitude,
4647 };
4648 let legacy = TrajectorySolver::new(
4649 BallisticInputs::default(),
4650 WindConditions::default(),
4651 atmosphere.clone(),
4652 );
4653
4654 let (resolved_temp_c, resolved_pressure_hpa) =
4655 crate::atmosphere::resolve_station_conditions_with_pressure_mode(
4656 temperature,
4657 pressure,
4658 altitude,
4659 crate::atmosphere::PressureReferenceMode::Absolute,
4660 );
4661 let precomputed_atmosphere = AtmosphericConditions {
4662 temperature: resolved_temp_c,
4663 pressure: resolved_pressure_hpa,
4664 humidity: 50.0,
4665 altitude,
4666 };
4667 let precomputed = TrajectorySolver::new_with_resolved_station_atmosphere(
4668 BallisticInputs::default(),
4669 WindConditions::default(),
4670 precomputed_atmosphere,
4671 );
4672
4673 let (legacy_density, legacy_sos, legacy_temp_c, legacy_pressure_hpa) =
4674 legacy.resolved_atmosphere();
4675 let (pre_density, pre_sos, pre_temp_c, pre_pressure_hpa) =
4676 precomputed.resolved_atmosphere();
4677
4678 assert_eq!(
4679 legacy_temp_c.to_bits(),
4680 pre_temp_c.to_bits(),
4681 "temperature=({temperature}, {pressure}, {altitude})"
4682 );
4683 assert_eq!(
4684 legacy_pressure_hpa.to_bits(),
4685 pre_pressure_hpa.to_bits(),
4686 "pressure=({temperature}, {pressure}, {altitude})"
4687 );
4688 assert_eq!(legacy_density.to_bits(), pre_density.to_bits());
4689 assert_eq!(legacy_sos.to_bits(), pre_sos.to_bits());
4690 }
4691 }
4692
4693 fn configured_euler_zero(vertical_wind_mps: f64, time_step_s: f64) -> TrajectorySolver {
4694 let inputs = BallisticInputs {
4695 muzzle_velocity: 800.0,
4696 bc_value: 0.5,
4697 bc_type: DragModel::G7,
4698 bullet_mass: 0.0109,
4699 bullet_diameter: 0.00782,
4700 bullet_length: 0.0309,
4701 sight_height: 0.05,
4702 ground_threshold: -100.0,
4703 use_rk4: false,
4704 use_adaptive_rk45: false,
4705 ..BallisticInputs::default()
4706 };
4707 let mut solver = TrajectorySolver::new_with_resolved_station_atmosphere(
4708 inputs,
4709 WindConditions::default(),
4710 AtmosphericConditions::default(),
4711 );
4712 solver.set_max_range(300.0);
4713 solver.set_time_step(time_step_s);
4714 if vertical_wind_mps != 0.0 {
4715 solver.set_wind_segments(vec![WindSegment {
4716 speed_kmh: 0.0,
4717 angle_deg: 0.0,
4718 until_m: 400.0,
4719 vertical_mps: vertical_wind_mps,
4720 }]);
4721 }
4722 solver
4723 }
4724
4725 #[test]
4726 fn inclined_shot_zeroes_like_a_level_rifle() {
4727 const ZERO_DISTANCE_M: f64 = 91.44; const SIGHT_HEIGHT_M: f64 = 0.0381; let inputs = BallisticInputs {
4736 bc_value: 0.5,
4737 bullet_mass: 150.0 * 0.06479891 / 1000.0,
4738 muzzle_velocity: 2700.0 * 0.3048,
4739 sight_height: SIGHT_HEIGHT_M,
4740 ..Default::default()
4741 };
4742
4743 let mut level = inputs.clone();
4744 level.shooting_angle = 0.0;
4745 let level_angle = TrajectorySolver::new(level, Default::default(), Default::default())
4746 .find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
4747 .expect("level zero must solve");
4748
4749 let mut inclined = inputs;
4750 inclined.shooting_angle = 5.71_f64.to_radians();
4751 let inclined_angle =
4752 TrajectorySolver::new(inclined, Default::default(), Default::default())
4753 .find_zero_angle(ZERO_DISTANCE_M, SIGHT_HEIGHT_M, ZeroTargetFrame::SightLine)
4754 .expect("MBA-1412: a 5.71 deg incline at a 100 yd zero must be solvable");
4755
4756 assert!(
4757 (inclined_angle - level_angle).abs() < 1e-9,
4758 "zeroing is level-rifle sight geometry; incline must not move the solved zero: \
4759 level={level_angle}, inclined={inclined_angle}"
4760 );
4761 }
4762
4763 #[test]
4764 fn configured_zero_keeps_segments_method_and_time_step_then_sets_base_angle() {
4765 const TARGET_DISTANCE_M: f64 = 150.0;
4766 const TARGET_HEIGHT_M: f64 = 0.05;
4767
4768 let mut segmented = configured_euler_zero(-10.0, 0.02);
4771 let coarse_height = segmented
4772 .zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
4773 .expect("coarse configured trial")
4774 .expect("coarse trial reaches target");
4775 let mut fine = segmented.clone();
4776 fine.set_time_step(0.001);
4777 let fine_height = fine
4778 .zero_trial_height_at(0.0, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
4779 .expect("fine configured trial")
4780 .expect("fine trial reaches target");
4781 assert!(
4782 (coarse_height - fine_height).abs() > 1e-5,
4783 "configured Euler step must affect zero trials: coarse={coarse_height}, fine={fine_height}"
4784 );
4785
4786 let segmented_angle = segmented
4787 .calculate_and_set_zero_angle(TARGET_DISTANCE_M, TARGET_HEIGHT_M, ZeroTargetFrame::SightLine)
4788 .expect("segmented zero");
4789 assert_eq!(
4790 segmented.inputs.muzzle_angle.to_bits(),
4791 segmented_angle.to_bits(),
4792 "successful zero must install its angle on the configured solver"
4793 );
4794 assert_eq!(segmented.time_step.to_bits(), 0.02_f64.to_bits());
4795 assert_eq!(segmented.max_range.to_bits(), 300.0_f64.to_bits());
4796 assert!(segmented.wind_sock.is_some());
4797 assert_eq!(
4798 segmented.station_atmosphere_resolution,
4799 StationAtmosphereResolution::Authoritative
4800 );
4801 let zero_height = segmented
4802 .zero_trial_height_at(segmented_angle, TARGET_DISTANCE_M, ZeroTargetFrame::SightLine)
4803 .expect("verify segmented zero")
4804 .expect("zeroed trial reaches target");
4805 assert!(
4806 (zero_height - TARGET_HEIGHT_M).abs() < 0.0001,
4807 "configured zero missed target: height={zero_height}"
4808 );
4809
4810 let mut calm = configured_euler_zero(0.0, 0.02);
4811 let calm_angle = calm
4812 .calculate_and_set_zero_angle(TARGET_DISTANCE_M, TARGET_HEIGHT_M, ZeroTargetFrame::SightLine)
4813 .expect("calm zero");
4814 assert!(
4815 (segmented_angle - calm_angle).abs() > 1e-5,
4816 "segmented vertical wind must participate in zero trials: segmented={segmented_angle}, calm={calm_angle}"
4817 );
4818 }
4819}
4820
4821#[cfg(test)]
4822mod result_sanity_tests {
4823 use super::*;
4824
4825 fn default_solver() -> TrajectorySolver {
4826 TrajectorySolver::new(
4827 BallisticInputs::default(),
4828 WindConditions::default(),
4829 AtmosphericConditions::default(),
4830 )
4831 }
4832
4833 fn minimal_result() -> TrajectoryResult {
4834 TrajectoryResult {
4835 max_range: 100.0,
4836 max_height: 1.0,
4837 time_of_flight: 0.5,
4838 impact_velocity: 700.0,
4839 impact_energy: 2450.0,
4840 projectile_mass_kg: 0.01,
4841 line_of_sight_height_m: 1.5,
4842 station_speed_of_sound_mps: 340.0,
4843 termination: TrajectoryTermination::MaxRange,
4844 points: vec![],
4845 sampled_points: None,
4846 min_pitch_damping: None,
4847 transonic_mach: None,
4848 angular_state: None,
4849 max_yaw_angle: None,
4850 max_precession_angle: None,
4851 aerodynamic_jump: None,
4852 mach_1_2_distance_m: None,
4853 mach_1_0_distance_m: None,
4854 mach_0_9_distance_m: None,
4855 }
4856 }
4857
4858 #[test]
4859 fn mba1293_negative_scalars_fail_the_result_postcondition() {
4860 let solver = default_solver();
4861 solver
4862 .validate_result_sanity(&minimal_result())
4863 .expect("a sane result must pass");
4864
4865 for (name, mutate) in [
4866 ("max_range", (|r| r.max_range = -50.588) as fn(&mut TrajectoryResult)),
4867 ("time_of_flight", |r| r.time_of_flight = -1.0),
4868 ("impact_velocity", |r| r.impact_velocity = -700.0),
4869 ("impact_energy", |r| r.impact_energy = -1.0),
4870 ] {
4871 let mut result = minimal_result();
4872 mutate(&mut result);
4873 let error = solver
4874 .validate_result_sanity(&result)
4875 .expect_err("negative scalar must fail");
4876 assert!(
4877 error.to_string().contains(name),
4878 "error for {name} did not name the field: {error}"
4879 );
4880 }
4881 }
4882
4883 #[test]
4884 fn mba1293_speed_budget_bounds_legitimate_states_and_rejects_divergence() {
4885 let solver = default_solver();
4886 let mv = solver.inputs.muzzle_velocity;
4887
4888 let position = Vector3::new(10.0, 0.0, 0.0);
4890 solver
4891 .validate_integration_state(&position, &Vector3::new(mv, 0.0, 0.0), 0.01)
4892 .expect("muzzle-speed state must pass");
4893
4894 let error = solver
4896 .validate_integration_state(&position, &Vector3::new(-13.0 * mv, 0.0, 0.0), 0.01)
4897 .expect_err("13x muzzle speed must fail the budget");
4898 assert!(error.to_string().contains("diverged"), "{error}");
4899
4900 let after_fall = mv + crate::constants::G_ACCEL_MPS2 * 60.0;
4902 solver
4903 .validate_integration_state(&position, &Vector3::new(0.0, -after_fall, 0.0), 60.0)
4904 .expect("gravity-accelerated speed within g*t must pass");
4905 }
4906}
4907
4908#[cfg(test)]
4909mod trajectory_point_budget_tests {
4910 use super::*;
4911 use crate::MAX_TRAJECTORY_SAMPLES;
4912
4913 fn solver_with_budget(
4914 use_rk4: bool,
4915 use_adaptive_rk45: bool,
4916 point_budget: usize,
4917 max_range: f64,
4918 ) -> TrajectorySolver {
4919 let inputs = BallisticInputs {
4920 use_rk4,
4921 use_adaptive_rk45,
4922 ground_threshold: f64::NEG_INFINITY,
4923 ..BallisticInputs::default()
4924 };
4925 let mut solver = TrajectorySolver::new(
4926 inputs,
4927 WindConditions::default(),
4928 AtmosphericConditions::default(),
4929 );
4930 solver.max_trajectory_points = point_budget;
4931 solver.set_max_range(max_range);
4932 solver.set_time_step(0.001);
4933 solver
4934 }
4935
4936 #[test]
4937 fn mba1283_every_solver_errors_instead_of_exceeding_point_budget() {
4938 for (mode, use_rk4, use_adaptive_rk45) in [
4939 ("Euler", false, false),
4940 ("RK4", true, false),
4941 ("RK45", true, true),
4942 ] {
4943 let error = solver_with_budget(use_rk4, use_adaptive_rk45, 3, 10.0)
4944 .solve()
4945 .expect_err("a solve requiring more than three points must fail");
4946 assert!(
4947 error.to_string().contains("point limit of 3"),
4948 "unexpected {mode} point-budget error: {error}"
4949 );
4950 }
4951 }
4952
4953 #[test]
4954 fn mba1283_interpolated_endpoint_counts_toward_point_budget() {
4955 for (mode, use_rk4, use_adaptive_rk45) in [
4956 ("Euler", false, false),
4957 ("RK4", true, false),
4958 ("RK45", true, true),
4959 ] {
4960 let result = solver_with_budget(use_rk4, use_adaptive_rk45, 2, 0.1)
4961 .solve()
4962 .expect("the initial point plus exact endpoint fit a two-point budget");
4963 assert_eq!(result.points.len(), 2, "unexpected {mode} point count");
4964
4965 let error = solver_with_budget(use_rk4, use_adaptive_rk45, 1, 0.1)
4966 .solve()
4967 .expect_err("the exact endpoint must not exceed a one-point budget");
4968 assert!(
4969 error.to_string().contains("point limit of 1"),
4970 "unexpected {mode} endpoint-budget error: {error}"
4971 );
4972 }
4973 }
4974
4975 #[test]
4976 fn mba1299_every_solver_preflights_the_sample_budget() {
4977 for (mode, use_rk4, use_adaptive_rk45) in [
4978 ("Euler", false, false),
4979 ("RK4", true, false),
4980 ("RK45", true, true),
4981 ] {
4982 let inputs = BallisticInputs {
4983 use_rk4,
4984 use_adaptive_rk45,
4985 enable_trajectory_sampling: true,
4986 sample_interval: 1.0,
4987 ground_threshold: f64::NEG_INFINITY,
4988 ..BallisticInputs::default()
4989 };
4990 let mut solver = TrajectorySolver::new(
4991 inputs,
4992 WindConditions::default(),
4993 AtmosphericConditions::default(),
4994 );
4995 solver.set_max_range(MAX_TRAJECTORY_SAMPLES as f64);
4996 solver.max_trajectory_points = 0;
4999
5000 let error = solver
5001 .solve()
5002 .expect_err("an over-limit sample grid must fail before integration");
5003 assert!(
5004 error
5005 .to_string()
5006 .contains("trajectory sample limit of 250000 exceeded"),
5007 "unexpected {mode} sample-budget error: {error}"
5008 );
5009 }
5010 }
5011
5012 #[test]
5013 fn mba1299_normal_sampling_does_not_change_solver_results() {
5014 for (mode, use_rk4, use_adaptive_rk45) in [
5015 ("Euler", false, false),
5016 ("RK4", true, false),
5017 ("RK45", true, true),
5018 ] {
5019 let solve = |enable_trajectory_sampling| {
5020 let inputs = BallisticInputs {
5021 use_rk4,
5022 use_adaptive_rk45,
5023 enable_trajectory_sampling,
5024 sample_interval: 0.5,
5025 ground_threshold: f64::NEG_INFINITY,
5026 ..BallisticInputs::default()
5027 };
5028 let mut solver = TrajectorySolver::new(
5029 inputs,
5030 WindConditions::default(),
5031 AtmosphericConditions::default(),
5032 );
5033 solver.set_max_range(2.0);
5034 solver.solve().expect("normal short-range solve")
5035 };
5036
5037 let baseline = solve(false);
5038 let sampled = solve(true);
5039 for (field, left, right) in [
5040 ("max_range", baseline.max_range, sampled.max_range),
5041 ("max_height", baseline.max_height, sampled.max_height),
5042 (
5043 "time_of_flight",
5044 baseline.time_of_flight,
5045 sampled.time_of_flight,
5046 ),
5047 (
5048 "impact_velocity",
5049 baseline.impact_velocity,
5050 sampled.impact_velocity,
5051 ),
5052 (
5053 "impact_energy",
5054 baseline.impact_energy,
5055 sampled.impact_energy,
5056 ),
5057 ] {
5058 assert_eq!(
5059 left.to_bits(),
5060 right.to_bits(),
5061 "{mode} sampling changed {field}"
5062 );
5063 }
5064 assert_eq!(baseline.points.len(), sampled.points.len());
5065 for (index, (left, right)) in baseline
5066 .points
5067 .iter()
5068 .zip(&sampled.points)
5069 .enumerate()
5070 {
5071 assert_eq!(left.time.to_bits(), right.time.to_bits(), "{mode} point {index}");
5072 assert_eq!(
5073 left.position.map(f64::to_bits),
5074 right.position.map(f64::to_bits),
5075 "{mode} point {index} position"
5076 );
5077 assert_eq!(
5078 left.velocity_magnitude.to_bits(),
5079 right.velocity_magnitude.to_bits(),
5080 "{mode} point {index} velocity"
5081 );
5082 assert_eq!(
5083 left.kinetic_energy.to_bits(),
5084 right.kinetic_energy.to_bits(),
5085 "{mode} point {index} energy"
5086 );
5087 }
5088 assert!(baseline.sampled_points.is_none());
5089 let samples = sampled
5090 .sampled_points
5091 .expect("sampling-enabled solve should return observations");
5092 assert_eq!(
5093 samples
5094 .iter()
5095 .map(|sample| sample.distance_m)
5096 .collect::<Vec<_>>(),
5097 vec![0.0, 0.5, 1.0, 1.5, 2.0],
5098 "{mode} normal sampling grid changed"
5099 );
5100 }
5101 }
5102}
5103
5104#[cfg(test)]
5105mod monte_carlo_result_tests {
5106 use super::*;
5107
5108 fn make_results(impact_positions: Vec<Vector3<f64>>) -> MonteCarloResults {
5109 let count = impact_positions.len();
5110 MonteCarloResults {
5111 ranges: vec![500.0; count],
5112 impact_velocities: vec![300.0; count],
5113 impact_positions,
5114 }
5115 }
5116
5117 #[test]
5118 fn target_plane_cep_excludes_shortfall_markers() {
5119 let mut positions: Vec<Vector3<f64>> = (1..=5)
5120 .map(|radius| Vector3::new(0.0, radius as f64, 0.0))
5121 .collect();
5122 positions.extend(
5123 (0..5).map(|_| Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0)),
5124 );
5125 let results = make_results(positions);
5126
5127 assert_eq!(results.target_arrival_count(), 5);
5128 assert_eq!(results.target_shortfall_fraction(), 0.5);
5129 assert_eq!(results.target_plane_cep(), Some(3.0));
5130
5131 let one_shortfall = make_results(vec![
5132 Vector3::new(0.0, 1.0, 0.0),
5133 Vector3::new(0.0, 2.0, 0.0),
5134 Vector3::new(0.0, 3.0, 0.0),
5135 Vector3::new(0.0, 4.0, 0.0),
5136 Vector3::new(0.0, 5.0, 0.0),
5137 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5138 ]);
5139 assert_eq!(one_shortfall.target_plane_cep(), Some(3.0));
5140 }
5141
5142 #[test]
5143 fn all_shortfalls_have_no_cep_but_still_count_as_misses() {
5144 let all_shortfalls = make_results(vec![
5145 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5146 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5147 ]);
5148 assert_eq!(all_shortfalls.target_arrival_count(), 0);
5149 assert_eq!(all_shortfalls.target_shortfall_fraction(), 1.0);
5150 assert_eq!(all_shortfalls.target_plane_cep(), None);
5151 assert_eq!(all_shortfalls.hit_probability(0.3), 0.0);
5152
5153 let one_hit_one_shortfall = make_results(vec![
5154 Vector3::new(0.0, 0.1, 0.0),
5155 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5156 ]);
5157 assert_eq!(one_hit_one_shortfall.hit_probability(0.3), 0.5);
5158 }
5159
5160 #[test]
5162 fn rect_hit_probability_checks_independent_axis_halves() {
5163 let results = make_results(vec![
5164 Vector3::new(0.0, 0.1, 0.1),
5166 Vector3::new(0.0, 0.0, 0.2),
5168 Vector3::new(0.0, 0.0, 0.201),
5170 Vector3::new(0.0, 0.301, 0.0),
5172 Vector3::new(0.0, TARGET_NOT_REACHED_SENTINEL_M, 0.0),
5174 ]);
5175 assert!((results.rect_hit_probability(0.4, 0.6) - 0.4).abs() < 1e-12);
5177 }
5178
5179 #[test]
5180 fn rect_hit_probability_matches_circular_hit_probability_for_a_centered_hit() {
5181 let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
5182 assert_eq!(results.rect_hit_probability(0.5, 0.5), 1.0);
5183 assert_eq!(results.hit_probability(0.3), 1.0);
5184 }
5185
5186 #[test]
5187 fn rect_hit_probability_is_zero_for_empty_or_nonpositive_dimensions() {
5188 let empty = make_results(vec![]);
5189 assert_eq!(empty.rect_hit_probability(1.0, 1.0), 0.0);
5190
5191 let results = make_results(vec![Vector3::new(0.0, 0.0, 0.0)]);
5192 assert_eq!(results.rect_hit_probability(0.0, 1.0), 0.0);
5193 assert_eq!(results.rect_hit_probability(1.0, 0.0), 0.0);
5194 assert_eq!(results.rect_hit_probability(-1.0, 1.0), 0.0);
5195 }
5196}
5197
5198#[cfg(test)]
5199mod monte_carlo_seeded_tests {
5200 use super::*;
5201
5202 #[test]
5203 fn seeded_runs_are_deterministic_and_match_the_using_rng_path() {
5204 let inputs = BallisticInputs {
5205 muzzle_velocity: 800.0,
5206 ..BallisticInputs::default()
5207 };
5208 let params = MonteCarloParams {
5209 num_simulations: 64,
5210 target_distance: Some(200.0),
5211 ..MonteCarloParams::default()
5212 };
5213
5214 let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5215 inputs.clone(),
5216 WindConditions::default(),
5217 params.clone(),
5218 0.01,
5219 42,
5220 )
5221 .expect("seeded run a");
5222 let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5223 inputs,
5224 WindConditions::default(),
5225 params,
5226 0.01,
5227 42,
5228 )
5229 .expect("seeded run b");
5230
5231 assert_eq!(a.ranges.len(), b.ranges.len());
5232 for (ra, rb) in a.ranges.iter().zip(b.ranges.iter()) {
5233 assert_eq!(ra.to_bits(), rb.to_bits());
5234 }
5235 for (pa, pb) in a.impact_positions.iter().zip(b.impact_positions.iter()) {
5236 assert_eq!(pa.x.to_bits(), pb.x.to_bits());
5237 assert_eq!(pa.y.to_bits(), pb.y.to_bits());
5238 assert_eq!(pa.z.to_bits(), pb.z.to_bits());
5239 }
5240 }
5241
5242 #[test]
5243 fn different_seeds_generally_produce_different_draws() {
5244 let inputs = BallisticInputs {
5245 muzzle_velocity: 800.0,
5246 ..BallisticInputs::default()
5247 };
5248 let params = MonteCarloParams {
5249 num_simulations: 32,
5250 velocity_std_dev: 5.0,
5251 target_distance: Some(200.0),
5252 ..MonteCarloParams::default()
5253 };
5254
5255 let a = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5256 inputs.clone(),
5257 WindConditions::default(),
5258 params.clone(),
5259 0.0,
5260 1,
5261 )
5262 .expect("seeded run a");
5263 let b = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
5264 inputs,
5265 WindConditions::default(),
5266 params,
5267 0.0,
5268 2,
5269 )
5270 .expect("seeded run b");
5271
5272 assert_ne!(a.impact_velocities, b.impact_velocities);
5273 }
5274}
5275
5276#[cfg(test)]
5277mod monte_carlo_powder_curve_tests {
5278 use super::*;
5279 use rand::{rngs::StdRng, SeedableRng};
5280
5281 #[test]
5282 fn powder_curve_preserves_sampled_muzzle_velocity_dispersion() {
5283 let inputs = BallisticInputs {
5284 muzzle_velocity: 700.0,
5285 powder_temp_curve: Some(vec![(15.0, 800.0)]),
5286 powder_curve_temp_c: Some(15.0),
5287 ..BallisticInputs::default()
5288 };
5289 let params = MonteCarloParams {
5290 num_simulations: 16,
5291 velocity_std_dev: 20.0,
5292 angle_std_dev: 1e-12,
5293 bc_std_dev: 1e-12,
5294 wind_speed_std_dev: 1e-12,
5295 target_distance: Some(100.0),
5296 azimuth_std_dev: 1e-12,
5297 ..MonteCarloParams::default()
5298 };
5299
5300 let mut rng = StdRng::seed_from_u64(0x5EED_1176);
5301 let results = run_monte_carlo_with_wind_and_direction_std_dev_using_rng(
5302 inputs,
5303 WindConditions::default(),
5304 params,
5305 0.0,
5306 &mut rng,
5307 )
5308 .expect("Monte Carlo solve");
5309 let min_velocity = results
5310 .impact_velocities
5311 .iter()
5312 .copied()
5313 .fold(f64::INFINITY, f64::min);
5314 let max_velocity = results
5315 .impact_velocities
5316 .iter()
5317 .copied()
5318 .fold(f64::NEG_INFINITY, f64::max);
5319
5320 assert!(
5321 max_velocity - min_velocity > 1.0,
5322 "20 m/s muzzle spread collapsed after curve resolution: impact-velocity span={} m/s",
5323 max_velocity - min_velocity
5324 );
5325 }
5326}
5327
5328#[cfg(test)]
5329mod monte_carlo_wind_sampling_tests {
5330 use super::*;
5331 use rand::{rngs::StdRng, SeedableRng};
5332
5333 #[test]
5334 fn wind_speed_sigma_does_not_change_seeded_direction_draws() {
5335 let base_wind = WindConditions {
5336 speed: 100.0,
5337 direction: 0.37,
5338 vertical_speed: 0.0,
5339 };
5340 let narrow_speed = MonteCarloWindSampler::new(&base_wind, 0.5, 0.2).unwrap();
5341 let wide_speed = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
5342 let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
5343 let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
5344 let mut speed_changed = false;
5345
5346 for _ in 0..32 {
5347 let narrow = narrow_speed.sample(&mut narrow_rng);
5348 let wide = wide_speed.sample(&mut wide_rng);
5349 assert!(narrow.speed > 0.0 && wide.speed > 0.0);
5350 assert_eq!(narrow.direction.to_bits(), wide.direction.to_bits());
5351 speed_changed |= narrow.speed.to_bits() != wide.speed.to_bits();
5352 }
5353 assert!(
5354 speed_changed,
5355 "different speed sigmas must still vary speed draws"
5356 );
5357 }
5358
5359 #[test]
5360 fn zero_direction_sigma_has_no_angular_jitter() {
5361 let base_wind = WindConditions {
5362 speed: 100.0,
5363 direction: 0.37,
5364 vertical_speed: 0.0,
5365 };
5366 let sampler = MonteCarloWindSampler::new(&base_wind, 4.0, 0.0).unwrap();
5367 let mut rng = StdRng::seed_from_u64(0x5EED_1223);
5368 let mut speed_changed = false;
5369
5370 for _ in 0..32 {
5371 let wind = sampler.sample(&mut rng);
5372 speed_changed |= wind.speed.to_bits() != base_wind.speed.to_bits();
5373 assert_eq!(wind.direction.to_bits(), base_wind.direction.to_bits());
5374 }
5375 assert!(speed_changed, "speed uncertainty should remain active");
5376 }
5377
5378 #[test]
5379 fn direction_sigma_controls_seeded_angular_spread_in_radians() {
5380 let base_wind = WindConditions {
5381 speed: 100.0,
5382 direction: 0.37,
5383 vertical_speed: 0.0,
5384 };
5385 let narrow = MonteCarloWindSampler::new(&base_wind, 4.0, 0.1).unwrap();
5386 let wide = MonteCarloWindSampler::new(&base_wind, 4.0, 0.2).unwrap();
5387 let mut narrow_rng = StdRng::seed_from_u64(0x5EED_1223);
5388 let mut wide_rng = StdRng::seed_from_u64(0x5EED_1223);
5389 let mut nonzero_direction_draw = false;
5390
5391 for _ in 0..32 {
5392 let narrow_wind = narrow.sample(&mut narrow_rng);
5393 let wide_wind = wide.sample(&mut wide_rng);
5394 assert_eq!(narrow_wind.speed.to_bits(), wide_wind.speed.to_bits());
5395
5396 let narrow_delta = narrow_wind.direction - base_wind.direction;
5397 let wide_delta = wide_wind.direction - base_wind.direction;
5398 assert!((wide_delta - 2.0 * narrow_delta).abs() < 1e-12);
5399 nonzero_direction_draw |= narrow_delta.abs() > 1e-6;
5400 }
5401 assert!(
5402 nonzero_direction_draw,
5403 "positive radians sigma must vary direction"
5404 );
5405 }
5406
5407 #[test]
5408 fn direction_sigma_rejects_negative_or_nonfinite_values() {
5409 let base_wind = WindConditions::default();
5410 for sigma in [-0.1, f64::NAN, f64::INFINITY] {
5411 assert!(MonteCarloWindSampler::new(&base_wind, 1.0, sigma).is_err());
5412 }
5413 }
5414
5415 #[test]
5416 fn base_vertical_wind_rides_into_every_mc_sample() {
5417 use rand::SeedableRng;
5421 let base_wind = WindConditions { vertical_speed: 4.2, ..Default::default() };
5422 let sampler = MonteCarloWindSampler::new(&base_wind, 1.0, 0.2).unwrap();
5423 let mut rng = rand::rngs::StdRng::seed_from_u64(7);
5424 for _ in 0..32 {
5425 let w = sampler.sample(&mut rng);
5426 assert_eq!(w.vertical_speed, 4.2);
5427 }
5428 }
5429
5430 #[test]
5431 fn negative_speed_sample_reverses_wind_direction() {
5432 let direction = 0.25;
5433 let signed_speed = -2.5;
5434 let wind = wind_from_signed_speed_sample(signed_speed, direction, 0.0);
5435 let positive_wind = wind_from_signed_speed_sample(2.5, direction, 0.0);
5436
5437 assert_eq!(wind.speed, 2.5);
5438 assert!(
5439 (wind.direction - (direction + std::f64::consts::PI)).abs() < f64::EPSILON,
5440 "negative speed must reverse direction by pi: got {}",
5441 wind.direction
5442 );
5443 assert_eq!(positive_wind.speed, 2.5);
5444 assert_eq!(positive_wind.direction, direction);
5445
5446 let normalized_x = -wind.speed * wind.direction.cos();
5447 let normalized_z = -wind.speed * wind.direction.sin();
5448 let signed_x = -signed_speed * direction.cos();
5449 let signed_z = -signed_speed * direction.sin();
5450 assert!((normalized_x - signed_x).abs() < 1e-12);
5451 assert!((normalized_z - signed_z).abs() < 1e-12);
5452 }
5453}
5454
5455#[cfg(test)]
5456mod bc_fit_objective_tests {
5457 use super::*;
5458
5459 fn velocity_point(range_m: f64, velocity_mps: f64) -> TrajectoryPoint {
5460 TrajectoryPoint {
5461 time: 0.0,
5462 position: Vector3::new(range_m, 0.0, 0.0),
5463 velocity_magnitude: velocity_mps,
5464 kinetic_energy: 0.0,
5465 drag_coefficient: None,
5466 }
5467 }
5468
5469 #[test]
5470 fn candidate_that_misses_an_observation_has_no_score() {
5471 let trajectory = vec![velocity_point(0.0, 800.0), velocity_point(100.0, 700.0)];
5472 let observations = vec![(50.0, 750.0), (150.0, 600.0)];
5473
5474 assert!(
5475 fit_residual_sse(&trajectory, &observations, BcFitMode::Velocity, 0.0).is_none(),
5476 "a candidate that reaches only one of two observations must not compete on partial SSE"
5477 );
5478
5479 let complete_observations = vec![(50.0, 740.0), (100.0, 680.0)];
5480 assert_eq!(
5481 fit_residual_sse(
5482 &trajectory,
5483 &complete_observations,
5484 BcFitMode::Velocity,
5485 0.0,
5486 ),
5487 Some(500.0)
5488 );
5489 }
5490}
5491
5492#[cfg(test)]
5493mod cluster_bc_reference_space_tests {
5494 use super::*;
5495
5496 fn acceleration_at_1100_fps(inputs: BallisticInputs) -> Vector3<f64> {
5497 let solver = TrajectorySolver::new(
5498 inputs,
5499 WindConditions::default(),
5500 AtmosphericConditions::default(),
5501 );
5502 let position = Vector3::zeros();
5503 let velocity = Vector3::new(1100.0 / 3.28084, 0.0, 0.0);
5504 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
5505 solver.calculate_acceleration(
5506 &position,
5507 &velocity,
5508 &Vector3::zeros(),
5509 (temp_c, pressure_hpa, density / 1.225),
5510 )
5511 }
5512
5513 #[test]
5514 fn solver_passes_g7_reference_model_to_cluster_classifier() {
5515 let inputs = BallisticInputs {
5516 bc_value: 0.190,
5517 bc_type: DragModel::G7,
5518 bullet_mass: 77.0 * crate::constants::GRAINS_TO_KG,
5519 bullet_diameter: 0.224 * 0.0254,
5520 use_cluster_bc: true,
5521 ..BallisticInputs::default()
5522 };
5523
5524 let solver = TrajectorySolver::new(
5525 inputs,
5526 WindConditions::default(),
5527 AtmosphericConditions::default(),
5528 );
5529 let corrected = solver.apply_cluster_bc_correction(0.190, 2800.0);
5530
5531 assert!(
5532 (corrected / 0.190 - 1.004).abs() < 1e-12,
5533 "solver selected the wrong G7 cluster multiplier: {}",
5534 corrected / 0.190
5535 );
5536 }
5537
5538 #[test]
5539 fn velocity_bc_segments_are_not_cluster_corrected_twice() {
5540 let segmented_clustered = BallisticInputs {
5541 bc_value: 0.5,
5542 bc_type: DragModel::G7,
5543 use_bc_segments: true,
5544 bc_segments_data: Some(vec![
5545 crate::BCSegmentData {
5546 velocity_min: 0.0,
5547 velocity_max: 1_600.0,
5548 bc_value: 0.4,
5549 },
5550 crate::BCSegmentData {
5551 velocity_min: 1_600.0,
5552 velocity_max: 5_000.0,
5553 bc_value: 0.45,
5554 },
5555 ]),
5556 use_cluster_bc: true,
5557 ..BallisticInputs::default()
5558 };
5559 let mut segmented_only = segmented_clustered.clone();
5560 segmented_only.use_cluster_bc = false;
5561 let mut constant_clustered = segmented_clustered.clone();
5562 constant_clustered.bc_value = 0.4;
5563 constant_clustered.bc_segments_data = None;
5564
5565 let stacked = acceleration_at_1100_fps(segmented_clustered);
5566 let segment_only = acceleration_at_1100_fps(segmented_only);
5567 let cluster_only = acceleration_at_1100_fps(constant_clustered);
5568
5569 assert!(
5570 (stacked.x - segment_only.x).abs() < 1e-12,
5571 "segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
5572 stacked.x,
5573 segment_only.x
5574 );
5575 assert!(
5576 (cluster_only.x - segment_only.x).abs() > 1e-6,
5577 "cluster correction must remain active for a constant BC"
5578 );
5579 }
5580
5581 #[test]
5582 fn mach_bc_segments_are_not_cluster_corrected_twice() {
5583 let mach_segmented_clustered = BallisticInputs {
5584 bc_value: 0.5,
5585 bc_type: DragModel::G7,
5586 use_bc_segments: false,
5587 bc_segments: Some(vec![(0.5, 0.3), (1.5, 0.5)]),
5588 use_cluster_bc: true,
5589 ..BallisticInputs::default()
5590 };
5591 let mut mach_segmented_only = mach_segmented_clustered.clone();
5592 mach_segmented_only.use_cluster_bc = false;
5593
5594 let stacked = acceleration_at_1100_fps(mach_segmented_clustered);
5595 let segment_only = acceleration_at_1100_fps(mach_segmented_only);
5596
5597 assert!(
5598 (stacked.x - segment_only.x).abs() < 1e-12,
5599 "Mach segment BC already owns the velocity shape: stacked ax={} segment-only ax={}",
5600 stacked.x,
5601 segment_only.x
5602 );
5603 }
5604}
5605
5606#[cfg(test)]
5607mod velocity_bc_flag_tests {
5608 use super::*;
5609
5610 fn acceleration_at_600_mps(inputs: BallisticInputs) -> Vector3<f64> {
5611 let solver = TrajectorySolver::new(
5612 inputs,
5613 WindConditions::default(),
5614 AtmosphericConditions::default(),
5615 );
5616 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
5617 solver.calculate_acceleration(
5618 &Vector3::zeros(),
5619 &Vector3::new(600.0, 0.0, 0.0),
5620 &Vector3::zeros(),
5621 (temp_c, pressure_hpa, density / 1.225),
5622 )
5623 }
5624
5625 #[test]
5626 fn velocity_bc_data_requires_opt_in_in_trajectory_solver() {
5627 let scalar_inputs = BallisticInputs {
5628 bc_value: 0.5,
5629 bc_type: DragModel::G7,
5630 ..BallisticInputs::default()
5631 };
5632 let mut disabled_inputs = scalar_inputs.clone();
5633 disabled_inputs.bc_segments_data = Some(vec![crate::BCSegmentData {
5634 velocity_min: 0.0,
5635 velocity_max: 4_000.0,
5636 bc_value: 0.46,
5637 }]);
5638 disabled_inputs.use_bc_segments = false;
5639 let mut enabled_inputs = disabled_inputs.clone();
5640 enabled_inputs.use_bc_segments = true;
5641 let mut mach_only_inputs = scalar_inputs.clone();
5642 mach_only_inputs.bc_segments = Some(vec![(0.0, 0.4), (3.0, 0.4)]);
5643 let mut disabled_with_both = mach_only_inputs.clone();
5644 disabled_with_both.bc_segments_data = disabled_inputs.bc_segments_data.clone();
5645
5646 let scalar = acceleration_at_600_mps(scalar_inputs);
5647 let disabled = acceleration_at_600_mps(disabled_inputs);
5648 let enabled = acceleration_at_600_mps(enabled_inputs);
5649 let mach_only = acceleration_at_600_mps(mach_only_inputs);
5650 let disabled_with_both = acceleration_at_600_mps(disabled_with_both);
5651
5652 assert_eq!(
5653 disabled.x.to_bits(),
5654 scalar.x.to_bits(),
5655 "a populated velocity table must not change drag while use_bc_segments is false"
5656 );
5657 assert!(
5658 enabled.x < disabled.x - 1.0,
5659 "enabling the lower BC table must increase drag: disabled ax={} enabled ax={}",
5660 disabled.x,
5661 enabled.x
5662 );
5663 assert_eq!(
5664 disabled_with_both.x.to_bits(),
5665 mach_only.x.to_bits(),
5666 "disabling velocity data must fall through to an explicit Mach table"
5667 );
5668 }
5669}
5670
5671#[cfg(test)]
5672mod mach_bc_segment_tests {
5673 use super::*;
5674
5675 #[test]
5676 fn trajectory_solver_interpolates_explicit_mach_bc_segments() {
5677 let segmented_inputs = BallisticInputs {
5678 bc_value: 0.8,
5679 use_bc_segments: false,
5680 bc_segments: Some(vec![(1.0, 0.2), (2.0, 0.4)]),
5681 bc_segments_data: None,
5682 ..BallisticInputs::default()
5683 };
5684
5685 let mut expected_inputs = segmented_inputs.clone();
5686 expected_inputs.bc_value = 0.3;
5687 expected_inputs.bc_segments = None;
5688
5689 let atmosphere = AtmosphericConditions::default();
5690 let segmented_solver = TrajectorySolver::new(
5691 segmented_inputs,
5692 WindConditions::default(),
5693 atmosphere.clone(),
5694 );
5695 let expected_solver = TrajectorySolver::new(
5696 expected_inputs,
5697 WindConditions::default(),
5698 atmosphere,
5699 );
5700 let position = Vector3::zeros();
5701 let (density, _, temp_c, pressure_hpa) = segmented_solver.resolved_atmosphere();
5702 let (_, local_speed_of_sound) = crate::atmosphere::get_local_atmosphere_humid(
5703 segmented_solver.atmosphere.altitude,
5704 segmented_solver.atmosphere.altitude,
5705 temp_c,
5706 pressure_hpa,
5707 density / 1.225,
5708 segmented_solver.atmosphere.humidity,
5709 );
5710 let velocity = Vector3::new(1.5 * local_speed_of_sound, 0.0, 0.0);
5711 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
5712
5713 let segmented_acceleration = segmented_solver.calculate_acceleration(
5714 &position,
5715 &velocity,
5716 &Vector3::zeros(),
5717 resolved_atmo,
5718 );
5719 let expected_acceleration = expected_solver.calculate_acceleration(
5720 &position,
5721 &velocity,
5722 &Vector3::zeros(),
5723 resolved_atmo,
5724 );
5725
5726 assert!(
5727 (segmented_acceleration.x - expected_acceleration.x).abs() < 1e-12,
5728 "Mach 1.5 must interpolate BC 0.3: segmented ax={} expected ax={}",
5729 segmented_acceleration.x,
5730 expected_acceleration.x
5731 );
5732 }
5733}
5734
5735#[cfg(test)]
5736mod custom_drag_table_validation_tests {
5737 use super::*;
5738
5739 #[test]
5740 fn solve_accepts_zero_bc_when_custom_table_present() {
5741 let inputs = BallisticInputs {
5742 bc_value: 0.0, bullet_mass: 0.0106,
5744 bullet_diameter: 0.00782,
5745 muzzle_velocity: 850.0,
5746 custom_drag_table: Some(crate::drag::DragTable::new(
5747 vec![0.5, 1.0, 2.0, 3.0],
5748 vec![0.23, 0.40, 0.30, 0.26],
5749 )),
5750 ..BallisticInputs::default()
5751 };
5752 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
5753 assert!(solver.solve().is_ok());
5755 }
5756
5757 #[test]
5758 fn solve_still_requires_bc_without_table() {
5759 let inputs = BallisticInputs {
5760 bc_value: 0.0,
5761 bullet_mass: 0.0106,
5762 bullet_diameter: 0.00782,
5763 muzzle_velocity: 850.0,
5764 ..BallisticInputs::default()
5765 };
5766 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
5767 assert!(solver.solve().is_err());
5768 }
5769}
5770
5771#[cfg(test)]
5773mod cd_scale_tests {
5774 use super::*;
5775
5776 fn deck() -> crate::drag::DragTable {
5777 crate::drag::DragTable::new(vec![0.5, 1.0, 2.0, 3.0], vec![0.23, 0.40, 0.30, 0.26])
5778 }
5779
5780 fn deck_inputs(cd_scale: f64) -> BallisticInputs {
5781 BallisticInputs {
5782 bullet_mass: 0.0106,
5783 bullet_diameter: 0.00782,
5784 muzzle_velocity: 850.0,
5785 custom_drag_table: Some(deck()),
5786 cd_scale,
5787 ..BallisticInputs::default()
5788 }
5789 }
5790
5791 #[test]
5792 fn default_cd_scale_is_one() {
5793 assert_eq!(BallisticInputs::default().cd_scale, 1.0);
5794 }
5795
5796 #[test]
5799 fn cd_scale_absent_is_byte_identical_to_explicit_one() {
5800 let omitted = BallisticInputs {
5801 bullet_mass: 0.0106,
5802 bullet_diameter: 0.00782,
5803 muzzle_velocity: 850.0,
5804 custom_drag_table: Some(deck()),
5805 ..BallisticInputs::default()
5806 };
5807 let explicit = BallisticInputs {
5808 cd_scale: 1.0,
5809 ..omitted.clone()
5810 };
5811
5812 let solver_omitted =
5813 TrajectorySolver::new(omitted, WindConditions::default(), AtmosphericConditions::default());
5814 let solver_explicit =
5815 TrajectorySolver::new(explicit, WindConditions::default(), AtmosphericConditions::default());
5816
5817 let cd_omitted = solver_omitted.calculate_drag_coefficient(700.0, 340.0);
5818 let cd_explicit = solver_explicit.calculate_drag_coefficient(700.0, 340.0);
5819 assert_eq!(
5820 cd_omitted.to_bits(),
5821 cd_explicit.to_bits(),
5822 "default cd_scale must be bit-identical to an explicit 1.0"
5823 );
5824
5825 let result = solver_omitted.solve();
5828 assert!(result.is_ok(), "existing custom-deck solves must pass unchanged");
5829 }
5830
5831 #[test]
5833 fn cd_scale_multiplies_the_interpolated_cd_exactly() {
5834 let velocity = 700.0;
5835 let speed_of_sound = 340.0;
5836 let mach = velocity / speed_of_sound;
5837 let expected_unscaled = deck().interpolate(mach);
5838
5839 for &scale in &[0.90, 1.0, 1.10, 1.5] {
5840 let solver = TrajectorySolver::new(
5841 deck_inputs(scale),
5842 WindConditions::default(),
5843 AtmosphericConditions::default(),
5844 );
5845 let cd = solver.calculate_drag_coefficient(velocity, speed_of_sound);
5846 assert!(
5847 (cd - expected_unscaled * scale).abs() < 1e-12,
5848 "scale={scale}: cd={cd} expected={}",
5849 expected_unscaled * scale
5850 );
5851 }
5852 }
5853
5854 #[test]
5858 fn cd_scale_direction_on_cli_api_solver() {
5859 let solve = |scale: f64| {
5860 TrajectorySolver::new(
5861 deck_inputs(scale),
5862 WindConditions::default(),
5863 AtmosphericConditions::default(),
5864 )
5865 .solve()
5866 .expect("custom-deck solve should succeed")
5867 };
5868
5869 let baseline = solve(1.0);
5870 let scaled_up = solve(1.10);
5871 let scaled_down = solve(0.90);
5872
5873 assert!(
5874 scaled_up.impact_velocity < baseline.impact_velocity,
5875 "cd_scale=1.10 must increase drag -> lower impact velocity: base={} up={}",
5876 baseline.impact_velocity,
5877 scaled_up.impact_velocity
5878 );
5879 assert!(
5880 scaled_down.impact_velocity > baseline.impact_velocity,
5881 "cd_scale=0.90 must decrease drag -> higher impact velocity: base={} down={}",
5882 baseline.impact_velocity,
5883 scaled_down.impact_velocity
5884 );
5885 }
5886
5887 #[test]
5889 fn validate_for_solve_rejects_invalid_cd_scale() {
5890 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
5891 let solver = TrajectorySolver::new(
5892 deck_inputs(bad),
5893 WindConditions::default(),
5894 AtmosphericConditions::default(),
5895 );
5896 assert!(
5897 solver.solve().is_err(),
5898 "cd_scale={bad} must be rejected by validate_for_solve"
5899 );
5900 }
5901 }
5902
5903 #[test]
5910 fn validate_for_solve_rejects_invalid_cd_scale_without_a_custom_drag_table() {
5911 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
5912 let inputs = BallisticInputs {
5913 bc_value: 0.5,
5914 bc_type: crate::DragModel::G1,
5915 bullet_mass: 0.0106,
5916 bullet_diameter: 0.00782,
5917 muzzle_velocity: 850.0,
5918 cd_scale: bad,
5919 ..BallisticInputs::default()
5920 };
5921 assert!(inputs.custom_drag_table.is_none(), "precondition: no custom deck");
5922 let solver = TrajectorySolver::new(
5923 inputs,
5924 WindConditions::default(),
5925 AtmosphericConditions::default(),
5926 );
5927 assert!(
5928 solver.solve().is_err(),
5929 "cd_scale={bad} must be rejected by validate_for_solve even without a custom \
5930 drag table"
5931 );
5932 }
5933 }
5934
5935 #[test]
5939 fn cd_scale_is_inert_without_a_custom_drag_table() {
5940 let make = |cd_scale: f64| BallisticInputs {
5941 bc_value: 0.5,
5942 bc_type: crate::DragModel::G1,
5943 bullet_mass: 0.0106,
5944 bullet_diameter: 0.00782,
5945 muzzle_velocity: 850.0,
5946 cd_scale,
5947 ..BallisticInputs::default()
5948 };
5949 let solver_neutral = TrajectorySolver::new(
5950 make(1.0),
5951 WindConditions::default(),
5952 AtmosphericConditions::default(),
5953 );
5954 let solver_far = TrajectorySolver::new(
5955 make(1.5),
5956 WindConditions::default(),
5957 AtmosphericConditions::default(),
5958 );
5959 let cd_neutral = solver_neutral.calculate_drag_coefficient(700.0, 340.0);
5960 let cd_far = solver_far.calculate_drag_coefficient(700.0, 340.0);
5961 assert_eq!(
5962 cd_neutral.to_bits(),
5963 cd_far.to_bits(),
5964 "cd_scale must not affect the G-model/BC drag path"
5965 );
5966 }
5967
5968 #[test]
5971 fn cd_scale_shifts_all_three_solver_paths_in_the_same_direction() {
5972 let cli_solve = |scale: f64| {
5974 TrajectorySolver::new(
5975 deck_inputs(scale),
5976 WindConditions::default(),
5977 AtmosphericConditions::default(),
5978 )
5979 .solve()
5980 .expect("cli_api custom-deck solve should succeed")
5981 };
5982 let cli_baseline = cli_solve(1.0);
5983 let cli_scaled = cli_solve(1.10);
5984 assert!(
5985 cli_scaled.impact_velocity < cli_baseline.impact_velocity,
5986 "cli_api: cd_scale=1.10 must lower impact velocity"
5987 );
5988
5989 let derivatives_accel_x = |scale: f64| {
5991 let inputs = deck_inputs(scale);
5992 crate::derivatives::compute_derivatives(
5993 nalgebra::Vector3::zeros(),
5994 nalgebra::Vector3::new(700.0, 0.0, 0.0),
5995 &inputs,
5996 nalgebra::Vector3::zeros(),
5997 (1.225, 340.0, 0.0, 0.0),
5998 inputs.bc_value,
5999 None,
6000 0.0,
6001 None,
6002 )[3]
6003 };
6004 let deriv_baseline = derivatives_accel_x(1.0);
6005 let deriv_scaled = derivatives_accel_x(1.10);
6006 assert!(
6007 deriv_scaled < deriv_baseline,
6008 "derivatives: cd_scale=1.10 must make x-acceleration more negative (more drag): \
6009 base={deriv_baseline} scaled={deriv_scaled}"
6010 );
6011
6012 let fast_final_speed = |scale: f64| {
6014 let inputs = deck_inputs(scale);
6015 let wind_sock = crate::wind::WindSock::new(vec![]);
6016 let params = crate::fast_trajectory::FastIntegrationParams {
6017 horiz: 500.0,
6018 vert: 0.0,
6019 initial_state: [0.0, 0.0, 0.0, 850.0, 0.0, 0.0],
6020 t_span: (0.0, 5.0),
6021 atmo_params: (0.0, 15.0, 1013.25, 1.0),
6022 atmo_sock: None,
6023 };
6024 let solution = crate::fast_trajectory::fast_integrate(&inputs, &wind_sock, params);
6025 assert!(solution.success, "fast_integrate must succeed for scale={scale}");
6026 let last = solution.t.len() - 1;
6027 let (vx, vy, vz) = (
6028 solution.y[3][last],
6029 solution.y[4][last],
6030 solution.y[5][last],
6031 );
6032 (vx * vx + vy * vy + vz * vz).sqrt()
6033 };
6034 let fast_baseline = fast_final_speed(1.0);
6035 let fast_scaled = fast_final_speed(1.10);
6036 assert!(
6037 fast_scaled < fast_baseline,
6038 "fast_trajectory: cd_scale=1.10 must lower final speed: base={fast_baseline} scaled={fast_scaled}"
6039 );
6040 }
6041}
6042
6043#[cfg(test)]
6044mod humid_local_mach_tests {
6045 use super::*;
6046
6047 fn solver_with_station_humidity(humidity_percent: f64) -> TrajectorySolver {
6048 let inputs = BallisticInputs {
6049 custom_drag_table: Some(crate::drag::DragTable::new(vec![0.5, 1.5], vec![0.1, 1.1])),
6050 ..BallisticInputs::default()
6051 };
6052 TrajectorySolver::new(
6053 inputs,
6054 WindConditions::default(),
6055 AtmosphericConditions {
6056 temperature: 30.0,
6057 pressure: 1013.25,
6058 humidity: humidity_percent,
6059 altitude: 0.0,
6060 },
6061 )
6062 }
6063
6064 fn acceleration(solver: &TrajectorySolver, base_ratio: f64) -> Vector3<f64> {
6065 solver.calculate_acceleration(
6066 &Vector3::zeros(),
6067 &Vector3::new(350.0, 0.0, 0.0),
6068 &Vector3::zeros(),
6069 (30.0, 1013.25, base_ratio),
6070 )
6071 }
6072
6073 #[test]
6074 fn local_mach_uses_station_humidity_when_density_is_held_constant() {
6075 let dry = acceleration(&solver_with_station_humidity(0.0), 1.0);
6076 let humid = acceleration(&solver_with_station_humidity(100.0), 1.0);
6077
6078 assert!(
6079 humid.x > dry.x,
6080 "humid sound speed should lower Mach and drag on the rising test curve: dry ax={} humid ax={}",
6081 dry.x,
6082 humid.x
6083 );
6084 }
6085
6086 #[test]
6087 fn active_atmosphere_zone_uses_zone_humidity_instead_of_station_humidity() {
6088 let zone_humidity = 80.0;
6089 let zone_ratio =
6090 crate::atmosphere::calculate_air_density_cimp(30.0, 1013.25, zone_humidity) / 1.225;
6091 let station_solver = solver_with_station_humidity(zone_humidity);
6092 let mut zoned_solver = solver_with_station_humidity(0.0);
6093 zoned_solver.set_atmo_segments(vec![(30.0, 1013.25, zone_humidity, 1_000.0)]);
6094
6095 let station = acceleration(&station_solver, zone_ratio);
6096 let zoned = acceleration(&zoned_solver, zone_ratio);
6097
6098 assert!(
6099 (zoned - station).norm() < 1e-12,
6100 "active zone T/P/RH should override the station atmosphere: station={station:?} zoned={zoned:?}"
6101 );
6102 }
6103}
6104
6105#[cfg(test)]
6106mod inclined_atmosphere_frame_tests {
6107 use super::*;
6108
6109 fn expected_shot_frame_vector(level: Vector3<f64>, angle: f64) -> Vector3<f64> {
6110 let (sin_angle, cos_angle) = angle.sin_cos();
6111 Vector3::new(
6112 level.x * cos_angle + level.y * sin_angle,
6113 -level.x * sin_angle + level.y * cos_angle,
6114 level.z,
6115 )
6116 }
6117
6118 #[test]
6119 fn inclined_positions_at_same_world_altitude_have_same_solver_acceleration() {
6120 let angle = std::f64::consts::FRAC_PI_6;
6121 let inputs = BallisticInputs {
6122 shooting_angle: angle,
6123 ..BallisticInputs::default()
6124 };
6125 let atmosphere = AtmosphericConditions {
6126 altitude: 100.0,
6127 ..AtmosphericConditions::default()
6128 };
6129 let solver = TrajectorySolver::new(inputs, WindConditions::default(), atmosphere);
6130 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6131 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
6132 let velocity = Vector3::new(600.0, 0.0, 0.0);
6133 let along_slant = Vector3::new(1_000.0, 0.0, 0.0);
6134 let across_slant = Vector3::new(0.0, 500.0 / angle.cos(), 0.0);
6135
6136 let a = solver.calculate_acceleration(
6137 &along_slant,
6138 &velocity,
6139 &Vector3::zeros(),
6140 resolved_atmo,
6141 );
6142 let b = solver.calculate_acceleration(
6143 &across_slant,
6144 &velocity,
6145 &Vector3::zeros(),
6146 resolved_atmo,
6147 );
6148
6149 assert!(
6150 (a - b).norm() < 1e-10,
6151 "solver acceleration differs at equal world altitude: {a:?} vs {b:?}"
6152 );
6153 }
6154
6155 #[test]
6156 fn inclined_headwind_is_rotated_into_solver_frame() {
6157 let angle = std::f64::consts::FRAC_PI_6;
6158 let inputs = BallisticInputs {
6159 shooting_angle: angle,
6160 ..BallisticInputs::default()
6161 };
6162 let solver = TrajectorySolver::new(
6163 inputs,
6164 WindConditions::default(),
6165 AtmosphericConditions::default(),
6166 );
6167 let level_headwind = Vector3::new(-100.0, 0.0, 0.0);
6168 let velocity = expected_shot_frame_vector(level_headwind, angle);
6169 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6170 let actual = solver.calculate_acceleration(
6171 &Vector3::zeros(),
6172 &velocity,
6173 &level_headwind,
6174 (temp_c, pressure_hpa, density / 1.225),
6175 );
6176
6177 assert!(
6178 (actual - solver.gravity_acceleration()).norm() < 1e-12,
6179 "co-moving horizontal wind must leave only shot-frame gravity: {actual:?}"
6180 );
6181 }
6182
6183 #[test]
6184 fn inclined_coriolis_is_rotated_into_solver_frame() {
6185 let angle = std::f64::consts::FRAC_PI_6;
6186 let latitude_deg = 45.0_f64;
6187 let shot_azimuth = 0.4_f64;
6188 let velocity = Vector3::new(600.0, 20.0, 5.0);
6189 let base_inputs = BallisticInputs {
6190 shooting_angle: angle,
6191 latitude: Some(latitude_deg),
6192 shot_azimuth,
6193 ..BallisticInputs::default()
6194 };
6195 let acceleration = |enable_coriolis| {
6196 let mut inputs = base_inputs.clone();
6197 inputs.enable_coriolis = enable_coriolis;
6198 let solver = TrajectorySolver::new(
6199 inputs,
6200 WindConditions::default(),
6201 AtmosphericConditions::default(),
6202 );
6203 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6204 solver.calculate_acceleration(
6205 &Vector3::zeros(),
6206 &velocity,
6207 &Vector3::zeros(),
6208 (temp_c, pressure_hpa, density / 1.225),
6209 )
6210 };
6211
6212 let omega_earth = 7.2921159e-5_f64;
6213 let latitude = latitude_deg.to_radians();
6214 let level_omega = Vector3::new(
6215 omega_earth * latitude.cos() * shot_azimuth.cos(),
6216 omega_earth * latitude.sin(),
6217 -omega_earth * latitude.cos() * shot_azimuth.sin(),
6218 );
6219 let expected = -2.0 * expected_shot_frame_vector(level_omega, angle).cross(&velocity);
6220 let actual = acceleration(true) - acceleration(false);
6221
6222 assert!(
6223 (actual - expected).norm() < 1e-12,
6224 "inclined Coriolis mismatch: actual={actual:?}, expected={expected:?}"
6225 );
6226 }
6227}
6228
6229#[cfg(test)]
6230mod terminal_range_interpolation_tests {
6231 use super::*;
6232
6233 #[test]
6234 fn terminal_finalizer_selects_the_earliest_crossed_boundary() {
6235 let inputs = BallisticInputs {
6236 ground_threshold: 0.0,
6237 ..BallisticInputs::default()
6238 };
6239 let mut solver = TrajectorySolver::new(
6240 inputs,
6241 WindConditions::default(),
6242 AtmosphericConditions::default(),
6243 );
6244 solver.set_max_range(120.0);
6245
6246 let previous_speed = 700.0;
6247 let mut points = vec![TrajectoryPoint {
6248 time: 99.0,
6249 position: Vector3::new(90.0, 1.0, -1.0),
6250 velocity_magnitude: previous_speed,
6251 kinetic_energy: 0.5 * solver.inputs.bullet_mass * previous_speed.powi(2),
6252 drag_coefficient: None,
6253 }];
6254 let mut max_height = 1.0;
6255 let termination = solver
6256 .append_terminal_endpoint(
6257 &mut points,
6258 Vector3::new(130.0, -3.0, 3.0),
6259 Vector3::new(600.0, 0.0, 0.0),
6260 101.0,
6261 &mut max_height,
6262 )
6263 .expect("the final step brackets supported boundaries");
6264
6265 assert_eq!(termination, TrajectoryTermination::GroundThreshold);
6266 assert_eq!(points.len(), 2);
6267 let terminal = points.last().expect("terminal point");
6268 assert_eq!(terminal.time, 99.5);
6269 assert_eq!(terminal.position, Vector3::new(100.0, 0.0, 0.0));
6270 assert_eq!(terminal.velocity_magnitude, 675.0);
6271 assert_eq!(
6272 terminal.kinetic_energy,
6273 0.5 * solver.inputs.bullet_mass * 675.0_f64.powi(2)
6274 );
6275
6276 solver.set_max_range(100.0);
6278 let mut tied_points = vec![points[0].clone()];
6279 assert_eq!(
6280 solver
6281 .append_terminal_endpoint(
6282 &mut tied_points,
6283 Vector3::new(130.0, -3.0, 3.0),
6284 Vector3::new(600.0, 0.0, 0.0),
6285 101.0,
6286 &mut max_height,
6287 )
6288 .expect("tied boundaries remain a valid terminal"),
6289 TrajectoryTermination::GroundThreshold
6290 );
6291 }
6292
6293 #[test]
6294 fn sub_ulp_terminal_crossing_replaces_instead_of_duplicating_range() {
6295 let ground_threshold = f64::from_bits(1.0_f64.to_bits() - 1);
6296 let inputs = BallisticInputs {
6297 ground_threshold,
6298 ..BallisticInputs::default()
6299 };
6300 let mut solver = TrajectorySolver::new(
6301 inputs,
6302 WindConditions::default(),
6303 AtmosphericConditions::default(),
6304 );
6305 solver.set_max_range(1_000.0);
6306
6307 let speed = 700.0;
6308 let mut points = vec![TrajectoryPoint {
6309 time: 0.0,
6310 position: Vector3::new(100.0, 1.0, 0.0),
6311 velocity_magnitude: speed,
6312 kinetic_energy: 0.5 * solver.inputs.bullet_mass * speed.powi(2),
6313 drag_coefficient: None,
6314 }];
6315 let mut max_height = 1.0;
6316 let termination = solver
6317 .append_terminal_endpoint(
6318 &mut points,
6319 Vector3::new(101.0, 0.0, 0.0),
6320 Vector3::new(699.0, 0.0, 0.0),
6321 1.0,
6322 &mut max_height,
6323 )
6324 .expect("sub-ULP ground crossing remains representable as one terminal state");
6325
6326 assert_eq!(termination, TrajectoryTermination::GroundThreshold);
6327 assert_eq!(points.len(), 1);
6328 assert_eq!(points[0].position.x, 100.0);
6329 assert_eq!(points[0].position.y.to_bits(), ground_threshold.to_bits());
6330 assert!(points[0].time > 0.0);
6331 }
6332
6333 #[test]
6334 fn every_solver_appends_an_exact_max_range_endpoint() {
6335 let target_range = 0.1;
6336 let modes = [
6337 ("Euler", false, false),
6338 ("RK4", true, false),
6339 ("RK45", true, true),
6340 ];
6341
6342 for (name, use_rk4, use_adaptive_rk45) in modes {
6343 let inputs = BallisticInputs {
6344 use_rk4,
6345 use_adaptive_rk45,
6346 ground_threshold: f64::NEG_INFINITY,
6347 enable_trajectory_sampling: true,
6348 sample_interval: target_range,
6349 ..BallisticInputs::default()
6350 };
6351 let mut solver = TrajectorySolver::new(
6352 inputs,
6353 WindConditions::default(),
6354 AtmosphericConditions::default(),
6355 );
6356 solver.set_max_range(target_range);
6357
6358 let result = solver.solve().expect("short-range solve should succeed");
6359 let terminal = result.points.last().expect("terminal point is missing");
6360 let muzzle = result.points.first().expect("muzzle point is missing");
6361
6362 assert_eq!(result.termination, TrajectoryTermination::MaxRange);
6363 assert_eq!(
6364 terminal.position.x.to_bits(),
6365 target_range.to_bits(),
6366 "{name} did not terminate exactly at max_range"
6367 );
6368 assert_eq!(result.max_range.to_bits(), target_range.to_bits());
6369 assert!(
6370 result.time_of_flight > 0.0 && result.time_of_flight < solver.time_step,
6371 "{name} terminal time was not interpolated within the crossing step: {}",
6372 result.time_of_flight
6373 );
6374 assert_eq!(result.time_of_flight.to_bits(), terminal.time.to_bits());
6375 assert_eq!(
6376 result.impact_velocity.to_bits(),
6377 terminal.velocity_magnitude.to_bits()
6378 );
6379 assert_eq!(
6380 result.impact_energy.to_bits(),
6381 terminal.kinetic_energy.to_bits()
6382 );
6383 let expected_energy = 0.5 * solver.inputs.bullet_mass * result.impact_velocity.powi(2);
6384 assert!((result.impact_energy - expected_energy).abs() < 1e-12);
6385 assert!(terminal.velocity_magnitude < muzzle.velocity_magnitude);
6386 assert!(terminal.kinetic_energy < muzzle.kinetic_energy);
6387
6388 let terminal_sample = result
6389 .sampled_points
6390 .as_ref()
6391 .and_then(|samples| samples.last())
6392 .expect("terminal trajectory sample is missing");
6393 assert_eq!(
6394 terminal_sample.distance_m.to_bits(),
6395 target_range.to_bits(),
6396 "{name} sampling did not include max_range"
6397 );
6398 assert_eq!(
6399 terminal_sample.time_s.to_bits(),
6400 result.time_of_flight.to_bits()
6401 );
6402 assert_eq!(
6403 terminal_sample.velocity_mps.to_bits(),
6404 result.impact_velocity.to_bits()
6405 );
6406 assert!((terminal_sample.energy_j - result.impact_energy).abs() < 1e-12);
6407 }
6408 }
6409}
6410
6411#[cfg(test)]
6412mod precession_inertia_wiring_tests {
6413 use super::*;
6414
6415 #[test]
6416 fn solver_uses_projectile_specific_moments_of_inertia() {
6417 let mass_kg = 55.0 * crate::constants::GRAINS_TO_KG;
6418 let caliber_m = 0.224 * 0.0254;
6419 let length_m = 0.75 * 0.0254;
6420 let inputs = BallisticInputs {
6421 bullet_mass: mass_kg,
6422 bullet_diameter: caliber_m,
6423 bullet_length: length_m,
6424 muzzle_velocity: 800.0,
6425 twist_rate: 7.0,
6426 enable_precession_nutation: true,
6427 use_rk4: false,
6428 use_adaptive_rk45: false,
6429 ..BallisticInputs::default()
6430 };
6431 let mut solver = TrajectorySolver::new(
6432 inputs,
6433 WindConditions::default(),
6434 AtmosphericConditions::default(),
6435 );
6436 solver.set_max_range(0.1);
6437
6438 let (air_density, speed_of_sound, _, _) = solver.resolved_atmosphere();
6439 let velocity_mps = solver.inputs.muzzle_velocity;
6440 let velocity_fps = velocity_mps * 3.28084;
6441 let twist_rate_ft = solver.inputs.twist_rate / 12.0;
6442 let spin_rate_rad_s = (velocity_fps / twist_rate_ft) * 2.0 * std::f64::consts::PI;
6443 let initial_state = AngularState {
6444 pitch_angle: 0.001,
6445 yaw_angle: 0.001,
6446 pitch_rate: 0.0,
6447 yaw_rate: 0.0,
6448 precession_angle: 0.0,
6449 nutation_phase: 0.0,
6450 };
6451 let params = PrecessionNutationParams {
6452 mass_kg,
6453 caliber_m,
6454 length_m,
6455 spin_rate_rad_s,
6456 spin_inertia: crate::spin_decay::calculate_moment_of_inertia(
6457 mass_kg, caliber_m, length_m, "ogive",
6458 ),
6459 transverse_inertia: crate::pitch_damping::calculate_transverse_moment_of_inertia(
6460 mass_kg, caliber_m, length_m, "ogive",
6461 ),
6462 velocity_mps,
6463 air_density_kg_m3: air_density,
6464 mach: velocity_mps / speed_of_sound,
6465 pitch_damping_coeff: PitchDampingCoefficients::default().subsonic,
6466 nutation_damping_factor: 0.05,
6467 };
6468 let expected = calculate_combined_angular_motion(
6469 ¶ms,
6470 &initial_state,
6471 0.0,
6472 solver.time_step,
6473 0.001,
6474 );
6475 let actual = solver
6476 .solve()
6477 .expect("one-step solve should succeed")
6478 .angular_state
6479 .expect("precession/nutation was enabled");
6480
6481 assert!(
6482 (actual.precession_angle - expected.precession_angle).abs() < 1e-15,
6483 "precession phase used the wrong inertia: actual={}, expected={}",
6484 actual.precession_angle,
6485 expected.precession_angle
6486 );
6487 assert!(
6488 (actual.nutation_phase - expected.nutation_phase).abs() < 1e-15,
6489 "nutation phase used the wrong inertia: actual={}, expected={}",
6490 actual.nutation_phase,
6491 expected.nutation_phase
6492 );
6493 }
6494}
6495
6496#[cfg(test)]
6497mod form_factor_drag_tests {
6498 use super::*;
6499
6500 fn acceleration_with_form_factor_flag(enabled: bool) -> Vector3<f64> {
6501 let inputs = BallisticInputs {
6502 bc_value: 0.462,
6503 bc_type: DragModel::G1,
6504 bullet_model: Some("168gr SMK Match".to_string()),
6505 use_form_factor: enabled,
6506 ..BallisticInputs::default()
6507 };
6508 let solver = TrajectorySolver::new(
6509 inputs,
6510 WindConditions::default(),
6511 AtmosphericConditions::default(),
6512 );
6513 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6514 solver.calculate_acceleration(
6515 &Vector3::zeros(),
6516 &Vector3::new(600.0, 0.0, 0.0),
6517 &Vector3::zeros(),
6518 (temp_c, pressure_hpa, density / 1.225),
6519 )
6520 }
6521
6522 #[test]
6523 fn measured_bc_drag_does_not_apply_name_based_form_factor_again() {
6524 let baseline = acceleration_with_form_factor_flag(false);
6525 let flagged = acceleration_with_form_factor_flag(true);
6526
6527 assert!(
6528 (flagged - baseline).norm() < 1e-12,
6529 "published BC already encodes form factor: baseline={baseline:?} flagged={flagged:?}"
6530 );
6531 }
6532}
6533
6534#[cfg(test)]
6535mod rk45_adaptivity_tests {
6536 use super::*;
6537
6538 #[test]
6539 fn cli_rk45_error_norm_scales_components_independently() {
6540 let position = Vector3::new(1.0e9, 0.0, 0.0);
6541 let velocity = Vector3::new(800.0, 0.0, 0.0);
6542 let fifth_position = position;
6543 let fifth_velocity = velocity;
6544 let fourth_position = position;
6545 let fourth_velocity = Vector3::new(800.0, 1.0e-3, 0.0);
6546
6547 let error = cli_rk45_error_norm(
6548 &position,
6549 &velocity,
6550 &fifth_position,
6551 &fifth_velocity,
6552 &fourth_position,
6553 &fourth_velocity,
6554 );
6555 let expected = 1.0e-3 / 6.0_f64.sqrt();
6556
6557 assert!(
6558 (error - expected).abs() <= 1e-15,
6559 "large downrange position masked a velocity-component error: {error}"
6560 );
6561 }
6562
6563 fn discontinuous_wind_solver() -> TrajectorySolver {
6564 let inputs = BallisticInputs::default();
6565 let mut solver = TrajectorySolver::new(
6566 inputs,
6567 WindConditions::default(),
6568 AtmosphericConditions::default(),
6569 );
6570 solver.set_wind_segments(vec![
6571 crate::wind::WindSegment::new(0.0, 90.0, 4.0),
6572 crate::wind::WindSegment::new(1_000.0, 90.0, 10_000.0),
6573 ]);
6574 solver
6575 }
6576
6577 #[test]
6578 fn rk45_retries_discontinuous_trial_before_advancing() {
6579 let solver = discontinuous_wind_solver();
6580 let position = Vector3::new(0.0, solver.inputs.muzzle_height, 0.0);
6581 let velocity = Vector3::new(solver.inputs.muzzle_velocity, 0.0, 0.0);
6582 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6583 let resolved_atmo = (temp_c, pressure_hpa, density / 1.225);
6584 let dt = 0.01;
6585
6586 let rejected_trial = solver.rk45_step(
6587 &position,
6588 &velocity,
6589 dt,
6590 &Vector3::zeros(),
6591 RK45_TOLERANCE,
6592 resolved_atmo,
6593 );
6594 assert!(
6595 rejected_trial.error > RK45_TOLERANCE,
6596 "discontinuous full step must exceed tolerance, got {}",
6597 rejected_trial.error
6598 );
6599
6600 let accepted = solver.adaptive_rk45_step(
6601 &position,
6602 &velocity,
6603 dt,
6604 &Vector3::zeros(),
6605 resolved_atmo,
6606 );
6607 assert!(accepted.used_dt < dt, "oversized trial was not retried");
6608 assert!(
6609 accepted.error <= RK45_TOLERANCE || accepted.used_dt <= RK45_MIN_DT,
6610 "accepted error {} exceeds tolerance at dt {}",
6611 accepted.error,
6612 accepted.used_dt
6613 );
6614
6615 let accepted_trial = solver.rk45_step(
6616 &position,
6617 &velocity,
6618 accepted.used_dt,
6619 &Vector3::zeros(),
6620 RK45_TOLERANCE,
6621 resolved_atmo,
6622 );
6623 assert_eq!(accepted.position, accepted_trial.position);
6624 assert_eq!(accepted.velocity, accepted_trial.velocity);
6625 assert!((RK45_MIN_DT..=RK45_MAX_DT).contains(&accepted.next_dt));
6626 }
6627}
6628
6629#[cfg(test)]
6630mod ground_termination_tests {
6631 use super::*;
6632 use crate::trajectory_observation::TrajectoryObservationFlag;
6633
6634 #[test]
6635 fn every_solver_reports_one_exact_early_ground_endpoint() {
6636 for (name, use_rk4, use_adaptive_rk45) in [
6637 ("Euler", false, false),
6638 ("RK4", true, false),
6639 ("RK45", true, true),
6640 ] {
6641 let inputs = BallisticInputs {
6642 muzzle_height: 1.0,
6643 muzzle_angle: -0.2,
6644 ground_threshold: 0.0,
6645 use_rk4,
6646 use_adaptive_rk45,
6647 ..BallisticInputs::default()
6648 };
6649 let mut solver = TrajectorySolver::new(
6650 inputs,
6651 WindConditions::default(),
6652 AtmosphericConditions::default(),
6653 );
6654 solver.set_max_range(1_000.0);
6655
6656 let result = solver.solve().expect("early-ground solve should succeed");
6657 let terminal = result.points.last().expect("terminal point is missing");
6658
6659 assert_eq!(result.termination, TrajectoryTermination::GroundThreshold);
6660 assert_eq!(terminal.position.y.to_bits(), 0.0_f64.to_bits());
6661 assert!(
6662 terminal.position.x < 1_000.0,
6663 "{name} incorrectly reached max range"
6664 );
6665 assert_eq!(result.max_range.to_bits(), terminal.position.x.to_bits());
6666 assert_eq!(
6667 result
6668 .points
6669 .iter()
6670 .filter(|point| point.position.y == 0.0)
6671 .count(),
6672 1,
6673 "{name} did not retain exactly one ground endpoint"
6674 );
6675
6676 let observations = result
6677 .sample_observations(1.0, 100)
6678 .expect("checked early-ground sampling should succeed");
6679 assert!(observations[..observations.len() - 1]
6680 .iter()
6681 .all(|observation| observation.distance_m < terminal.position.x));
6682 let terminal_observation = observations.last().expect("terminal observation");
6683 assert_eq!(
6684 terminal_observation.distance_m.to_bits(),
6685 terminal.position.x.to_bits()
6686 );
6687 assert!(terminal_observation
6688 .flags
6689 .contains(&TrajectoryObservationFlag::Terminal));
6690 assert!(terminal_observation
6691 .flags
6692 .contains(&TrajectoryObservationFlag::GroundThreshold));
6693 assert_eq!(
6694 observations
6695 .iter()
6696 .filter(|observation| observation
6697 .flags
6698 .contains(&TrajectoryObservationFlag::Terminal))
6699 .count(),
6700 1,
6701 "{name} repeated the terminal observation"
6702 );
6703 }
6704 }
6705
6706 #[test]
6711 fn rk4_and_rk45_descend_to_ground_threshold() {
6712 for adaptive in [false, true] {
6713 let inputs = BallisticInputs {
6714 muzzle_angle: 0.1, use_rk4: true,
6716 use_adaptive_rk45: adaptive,
6717 ..BallisticInputs::default()
6718 };
6719 assert_eq!(
6720 inputs.ground_threshold, -100.0,
6721 "default ground_threshold is -100 m"
6722 );
6723
6724 let mut solver = TrajectorySolver::new(
6725 inputs,
6726 WindConditions::default(),
6727 AtmosphericConditions::default(),
6728 );
6729 solver.set_max_range(1.0e7);
6731
6732 let result = solver.solve().expect("solve should succeed");
6733 let final_y = result
6734 .points
6735 .last()
6736 .expect("trajectory has points")
6737 .position
6738 .y;
6739 assert!(
6740 final_y < -1.0,
6741 "adaptive_rk45={adaptive}: final y = {final_y} m; a lofted shot should descend \
6742 past launch level toward the ground_threshold floor, not stop at y = 0"
6743 );
6744 }
6745 }
6746}
6747
6748#[cfg(test)]
6749mod magnus_stability_tests {
6750 use super::*;
6751
6752 #[test]
6753 fn yaw_of_repose_magnus_force_is_vertical_and_twist_signed() {
6754 let acceleration = |enable_magnus, is_twist_right| {
6755 let inputs = BallisticInputs {
6756 muzzle_velocity: 822.96,
6757 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
6758 bullet_diameter: 0.308 * 0.0254,
6759 bullet_length: 1.215 * 0.0254,
6760 twist_rate: 10.0,
6761 is_twist_right,
6762 enable_magnus,
6763 ..BallisticInputs::default()
6764 };
6765 let solver = TrajectorySolver::new(
6766 inputs,
6767 WindConditions::default(),
6768 AtmosphericConditions::default(),
6769 );
6770 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6771 solver.calculate_acceleration(
6772 &Vector3::zeros(),
6773 &Vector3::new(822.96, 0.0, 0.0),
6774 &Vector3::zeros(),
6775 (temp_c, pressure_hpa, density / 1.225),
6776 )
6777 };
6778
6779 let baseline = acceleration(false, true);
6780 let right_twist = acceleration(true, true) - baseline;
6781 let left_twist = acceleration(true, false) - baseline;
6782
6783 assert!(
6784 right_twist.y < 0.0,
6785 "right-hand Magnus must point down, got {right_twist:?}"
6786 );
6787 assert!(
6788 left_twist.y > 0.0,
6789 "left-hand Magnus must point up, got {left_twist:?}"
6790 );
6791 assert!((right_twist.y + left_twist.y).abs() < 1e-12);
6792 assert!(right_twist.x.abs() < 1e-12 && right_twist.z.abs() < 1e-12);
6793 assert!(left_twist.x.abs() < 1e-12 && left_twist.z.abs() < 1e-12);
6794 }
6795
6796 #[test]
6797 fn magnus_uses_velocity_corrected_muzzle_stability_gate() {
6798 let muzzle_velocity = 1_400.0 / 3.28084;
6799 let inputs = BallisticInputs {
6800 muzzle_velocity,
6801 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
6802 bullet_diameter: 0.308 * 0.0254,
6803 bullet_length: 1.215 * 0.0254,
6804 twist_rate: 15.0,
6805 enable_magnus: true,
6806 ..BallisticInputs::default()
6807 };
6808 let solver = TrajectorySolver::new(
6809 inputs.clone(),
6810 WindConditions::default(),
6811 AtmosphericConditions::default(),
6812 );
6813
6814 let bare_sg = crate::spin_drift::miller_stability(0.308, 168.0, 15.0, 1.215);
6815 let canonical_sg = solver.effective_spin_drift_sg();
6816 assert!(bare_sg > 1.0, "test requires bare Sg above the Magnus gate");
6817 assert!(
6818 canonical_sg < 1.0,
6819 "velocity-corrected Sg must be below the gate, got {canonical_sg}"
6820 );
6821
6822 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6823 let acceleration = solver.calculate_acceleration(
6824 &Vector3::zeros(),
6825 &Vector3::new(muzzle_velocity, 0.0, 0.0),
6826 &Vector3::zeros(),
6827 (temp_c, pressure_hpa, density / 1.225),
6828 );
6829 let mut baseline_inputs = inputs;
6830 baseline_inputs.enable_magnus = false;
6831 let baseline_solver = TrajectorySolver::new(
6832 baseline_inputs,
6833 WindConditions::default(),
6834 AtmosphericConditions::default(),
6835 );
6836 let baseline = baseline_solver.calculate_acceleration(
6837 &Vector3::zeros(),
6838 &Vector3::new(muzzle_velocity, 0.0, 0.0),
6839 &Vector3::zeros(),
6840 (temp_c, pressure_hpa, density / 1.225),
6841 );
6842
6843 assert_eq!(
6844 acceleration, baseline,
6845 "canonical Sg below 1 must suppress every Magnus acceleration component"
6846 );
6847 }
6848
6849 #[test]
6850 fn magnus_force_grows_as_fixed_spin_projectile_slows() {
6851 let inputs = BallisticInputs {
6852 muzzle_velocity: 800.0,
6853 bullet_mass: 168.0 * crate::constants::GRAINS_TO_KG,
6854 bullet_diameter: 0.308 * 0.0254,
6855 bullet_length: 1.215 * 0.0254,
6856 twist_rate: 12.0,
6857 enable_magnus: true,
6858 ..BallisticInputs::default()
6859 };
6860
6861 let magnus_acceleration = |speed_mps| {
6862 let evaluate = |enable_magnus| {
6863 let mut run_inputs = inputs.clone();
6864 run_inputs.enable_magnus = enable_magnus;
6865 let solver = TrajectorySolver::new(
6866 run_inputs,
6867 WindConditions::default(),
6868 AtmosphericConditions::default(),
6869 );
6870 let (density, _, temp_c, pressure_hpa) = solver.resolved_atmosphere();
6871 solver
6872 .calculate_acceleration(
6873 &Vector3::zeros(),
6874 &Vector3::new(speed_mps, 0.0, 0.0),
6875 &Vector3::zeros(),
6876 (temp_c, pressure_hpa, density / 1.225),
6877 )
6878 .y
6879 };
6880 (evaluate(true) - evaluate(false)).abs()
6881 };
6882
6883 let fast = magnus_acceleration(200.0);
6884 let slow = magnus_acceleration(100.0);
6885 let ratio = slow / fast;
6886 let expected_ratio = 2.0_f64.powf(5.0 / 3.0);
6887
6888 assert!(fast > 0.0 && slow > 0.0, "fast={fast}, slow={slow}");
6889 assert!(
6890 (ratio - expected_ratio).abs() < 1e-3,
6891 "fixed-spin Magnus acceleration must grow downrange; slow/fast={ratio}, \
6892 expected={expected_ratio}"
6893 );
6894 }
6895}
6896
6897#[cfg(test)]
6898mod coriolis_direction_tests {
6899 use super::*;
6900 use std::f64::consts::FRAC_PI_2;
6901
6902 #[test]
6903 fn supersonic_crossing_flags_a_positive_range_sample() {
6904 use crate::trajectory_sampling::TrajectoryFlag;
6908
6909 for (solver_name, use_rk4, use_adaptive_rk45) in [
6910 ("Euler", false, false),
6911 ("RK4", true, false),
6912 ("RK45", true, true),
6913 ] {
6914 let inputs = BallisticInputs {
6915 muzzle_velocity: 850.0,
6916 bc_value: 0.2,
6917 bc_type: DragModel::G7,
6918 muzzle_angle: 0.03,
6919 enable_trajectory_sampling: true,
6920 sample_interval: 50.0,
6921 use_rk4,
6922 use_adaptive_rk45,
6923 ..BallisticInputs::default()
6924 };
6925 let mut solver = TrajectorySolver::new(
6926 inputs,
6927 WindConditions::default(),
6928 AtmosphericConditions::default(),
6929 );
6930 solver.set_max_range(2000.0);
6931 let samples = solver
6932 .solve()
6933 .expect("supersonic solve should succeed")
6934 .sampled_points
6935 .expect("sampling was enabled");
6936 let flagged_distances: Vec<_> = samples
6937 .iter()
6938 .filter(|sample| sample.flags.contains(&TrajectoryFlag::MachTransition))
6939 .map(|sample| sample.distance_m)
6940 .collect();
6941
6942 assert!(
6943 !flagged_distances.is_empty()
6944 && flagged_distances.iter().all(|distance| *distance > 0.0),
6945 "{solver_name} must flag genuine crossings only at positive range: {flagged_distances:?}"
6946 );
6947 }
6948 }
6949
6950 #[test]
6951 fn subsonic_launch_does_not_flag_a_muzzle_transition() {
6952 use crate::trajectory_sampling::TrajectoryFlag;
6953
6954 for (solver_name, use_rk4, use_adaptive_rk45) in [
6955 ("Euler", false, false),
6956 ("RK4", true, false),
6957 ("RK45", true, true),
6958 ] {
6959 let inputs = BallisticInputs {
6960 muzzle_velocity: 250.0,
6961 muzzle_angle: 0.02,
6962 enable_trajectory_sampling: true,
6963 sample_interval: 25.0,
6964 use_rk4,
6965 use_adaptive_rk45,
6966 ..BallisticInputs::default()
6967 };
6968 let mut solver = TrajectorySolver::new(
6969 inputs,
6970 WindConditions::default(),
6971 AtmosphericConditions::default(),
6972 );
6973 solver.set_max_range(300.0);
6974 let samples = solver
6975 .solve()
6976 .expect("subsonic solve should succeed")
6977 .sampled_points
6978 .expect("sampling was enabled");
6979
6980 assert!(
6981 samples
6982 .iter()
6983 .all(|sample| !sample.flags.contains(&TrajectoryFlag::MachTransition)),
6984 "{solver_name} marked a Mach transition for a launch already below Mach 1"
6985 );
6986 }
6987 }
6988
6989 #[test]
6990 fn mach_transition_tracker_requires_a_downward_crossing() {
6991 fn record(mach_values: &[f64]) -> Vec<f64> {
6992 let mut tracker = MachTransitionTracker::default();
6993 let mut distances = Vec::new();
6994 for (index, mach) in mach_values.iter().copied().enumerate() {
6995 tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
6996 }
6997 distances
6998 }
6999
7000 assert!(record(&[0.9, 0.8, 0.7]).is_empty());
7001 assert_eq!(record(&[1.1, 1.05, 0.99]), vec![20.0]);
7002 assert_eq!(record(&[1.2, 1.19, 1.0, 0.99]), vec![10.0, 30.0]);
7003 assert_eq!(record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]), vec![20.0, 30.0]);
7004 assert!(record(&[1.3, f64::NAN, 1.1]).is_empty());
7005 }
7006
7007 #[test]
7008 fn mach_transition_tracker_labels_0_9_without_touching_the_flat_vec() {
7009 fn record(mach_values: &[f64]) -> (Vec<f64>, MachTransitionTracker) {
7015 let mut tracker = MachTransitionTracker::default();
7016 let mut distances = Vec::new();
7017 for (index, mach) in mach_values.iter().copied().enumerate() {
7018 tracker.record_downward_crossings(mach, index as f64 * 10.0, &mut distances);
7019 }
7020 (distances, tracker)
7021 }
7022
7023 let (distances, tracker) = record(&[0.9, 0.8, 0.7]);
7026 assert!(distances.is_empty()); assert_eq!(tracker.mach_1_2_distance_m, None);
7028 assert_eq!(tracker.mach_1_0_distance_m, None);
7029 assert_eq!(tracker.mach_0_9_distance_m, Some(10.0));
7030
7031 let (distances, tracker) = record(&[1.1, 1.05, 0.99]);
7033 assert_eq!(distances, vec![20.0]);
7034 assert_eq!(tracker.mach_1_2_distance_m, None);
7035 assert_eq!(tracker.mach_1_0_distance_m, Some(20.0));
7036 assert_eq!(tracker.mach_0_9_distance_m, None);
7037
7038 let (distances, tracker) = record(&[1.2, 1.19, 1.0, 0.99]);
7040 assert_eq!(distances, vec![10.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(10.0));
7042 assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
7043 assert_eq!(tracker.mach_0_9_distance_m, None);
7044
7045 let (distances, tracker) = record(&[0.9, 1.3, 1.1, 0.9, 1.3, 0.8]);
7048 assert_eq!(distances, vec![20.0, 30.0]); assert_eq!(tracker.mach_1_2_distance_m, Some(20.0));
7050 assert_eq!(tracker.mach_1_0_distance_m, Some(30.0));
7051 assert_eq!(tracker.mach_0_9_distance_m, Some(50.0));
7052 assert!(
7053 tracker.mach_1_2_distance_m < tracker.mach_1_0_distance_m
7054 && tracker.mach_1_0_distance_m < tracker.mach_0_9_distance_m,
7055 "labeled crossings must be strictly increasing downrange"
7056 );
7057
7058 let (distances, tracker) = record(&[1.3, f64::NAN, 1.1]);
7060 assert!(distances.is_empty());
7061 assert_eq!(tracker.mach_1_2_distance_m, None);
7062 assert_eq!(tracker.mach_1_0_distance_m, None);
7063 assert_eq!(tracker.mach_0_9_distance_m, None);
7064 }
7065
7066 #[test]
7067 fn humidity_percent_converts_and_clamps() {
7068 let mut i = BallisticInputs {
7070 humidity: 0.5,
7071 ..BallisticInputs::default()
7072 };
7073 assert!((i.humidity_percent() - 50.0).abs() < 1e-9, "0.5 -> 50%");
7074 i.humidity = 0.0;
7075 assert_eq!(i.humidity_percent(), 0.0);
7076 i.humidity = 1.0;
7077 assert_eq!(i.humidity_percent(), 100.0);
7078 i.humidity = 1.5; assert_eq!(i.humidity_percent(), 100.0);
7080 }
7081
7082 fn vertical_at(shot_azimuth: f64, range_m: f64) -> f64 {
7085 let inputs = BallisticInputs {
7086 muzzle_velocity: 800.0,
7087 bc_value: 0.5,
7088 bc_type: DragModel::G7,
7089 muzzle_angle: 0.02, enable_coriolis: true,
7091 latitude: Some(45.0),
7092 shot_azimuth,
7093 ground_threshold: f64::NEG_INFINITY, ..BallisticInputs::default()
7095 };
7096 let mut solver = TrajectorySolver::new(
7097 inputs,
7098 WindConditions::default(),
7099 AtmosphericConditions::default(),
7100 );
7101 solver.set_max_range(range_m + 50.0);
7102 let r = solver.solve().expect("solve");
7103 let pts = &r.points;
7104 for i in 1..pts.len() {
7105 if pts[i].position.x >= range_m {
7106 let p1 = &pts[i - 1];
7107 let p2 = &pts[i];
7108 let t = (range_m - p1.position.x) / (p2.position.x - p1.position.x);
7109 return p1.position.y + t * (p2.position.y - p1.position.y);
7110 }
7111 }
7112 panic!("range {range_m} not reached");
7113 }
7114
7115 #[test]
7120 fn eotvos_east_higher_than_west() {
7121 let range = 600.0;
7122 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!(
7126 east > west,
7127 "east ({east:.5}) must be higher than west ({west:.5}) at {range} m (Eötvös)"
7128 );
7129 assert!(
7130 east > north && north > west,
7131 "north ({north:.5}) must lie between east ({east:.5}) and west ({west:.5})"
7132 );
7133 assert!(
7134 (east - west) > 1e-3,
7135 "E-W vertical separation ({:.6} m) should be physically meaningful, not FP noise",
7136 east - west
7137 );
7138 }
7139
7140 #[test]
7148 fn labeled_mach_crossings_match_pinned_pre_change_flat_vec_across_solvers() {
7149 let cases = [
7151 ("Euler", false, false, 670.9878683238721_f64, 805.5274119916264_f64),
7152 ("RK4", true, false, 671.7257336844475_f64, 805.933409072171_f64),
7153 ("RK45", true, true, 672.4905711917901_f64, 806.5709746782849_f64),
7154 ];
7155
7156 for (solver_name, use_rk4, use_adaptive_rk45, expected_1_2, expected_1_0) in cases {
7157 let inputs = BallisticInputs {
7158 muzzle_velocity: 850.0,
7159 bc_value: 0.2,
7160 bc_type: DragModel::G7,
7161 muzzle_angle: 0.03,
7162 use_rk4,
7163 use_adaptive_rk45,
7164 ..BallisticInputs::default()
7165 };
7166 let mut solver = TrajectorySolver::new(
7167 inputs,
7168 WindConditions::default(),
7169 AtmosphericConditions::default(),
7170 );
7171 solver.set_max_range(2000.0);
7172 let result = solver.solve().expect("solve should succeed");
7173
7174 assert_eq!(
7175 result.mach_1_2_distance_m,
7176 Some(expected_1_2),
7177 "{solver_name}: mach_1_2_distance_m must match the pinned pre-change flat-Vec value"
7178 );
7179 assert_eq!(
7180 result.mach_1_0_distance_m,
7181 Some(expected_1_0),
7182 "{solver_name}: mach_1_0_distance_m must match the pinned pre-change flat-Vec value"
7183 );
7184
7185 let mach_1_2 = result.mach_1_2_distance_m.expect("crosses 1.2");
7186 let mach_1_0 = result.mach_1_0_distance_m.expect("crosses 1.0");
7187 let mach_0_9 = result
7188 .mach_0_9_distance_m
7189 .expect("this trajectory also goes past 0.9 within 2000 m");
7190 assert!(
7191 mach_1_2 < mach_1_0 && mach_1_0 < mach_0_9,
7192 "{solver_name}: labeled crossings must be strictly increasing downrange \
7193 (1.2={mach_1_2}, 1.0={mach_1_0}, 0.9={mach_0_9})"
7194 );
7195 }
7196 }
7197
7198 #[test]
7201 fn labeled_mach_crossings_are_none_for_a_fully_supersonic_trajectory() {
7202 for (solver_name, use_rk4, use_adaptive_rk45) in [
7203 ("Euler", false, false),
7204 ("RK4", true, false),
7205 ("RK45", true, true),
7206 ] {
7207 let inputs = BallisticInputs {
7208 muzzle_velocity: 850.0,
7209 bc_value: 0.2,
7210 bc_type: DragModel::G7,
7211 muzzle_angle: 0.03,
7212 use_rk4,
7213 use_adaptive_rk45,
7214 ..BallisticInputs::default()
7215 };
7216 let mut solver = TrajectorySolver::new(
7217 inputs,
7218 WindConditions::default(),
7219 AtmosphericConditions::default(),
7220 );
7221 solver.set_max_range(200.0);
7223 let result = solver.solve().expect("solve should succeed");
7224
7225 assert_eq!(
7226 result.mach_1_2_distance_m, None,
7227 "{solver_name}: must not report a 1.2 crossing that never happens"
7228 );
7229 assert_eq!(
7230 result.mach_1_0_distance_m, None,
7231 "{solver_name}: must not report a 1.0 crossing that never happens"
7232 );
7233 assert_eq!(
7234 result.mach_0_9_distance_m, None,
7235 "{solver_name}: must not report a 0.9 crossing that never happens"
7236 );
7237 }
7238 }
7239}
7240
7241#[cfg(test)]
7242mod cant_tests {
7243 use super::*;
7244
7245 fn base_inputs() -> BallisticInputs {
7246 BallisticInputs {
7247 muzzle_velocity: 800.0,
7248 bc_value: 0.5,
7249 bc_type: DragModel::G7,
7250 bullet_mass: 0.0109,
7251 bullet_diameter: 0.00782,
7252 bullet_length: 0.0309,
7253 sight_height: 0.05,
7254 twist_rate: 10.0,
7255 use_rk4: true,
7256 ..BallisticInputs::default()
7257 }
7258 }
7259
7260 fn solve_with(inputs: BallisticInputs, max_range: f64) -> TrajectoryResult {
7261 let mut s = TrajectorySolver::new(
7262 inputs,
7263 WindConditions::default(),
7264 AtmosphericConditions::default(),
7265 );
7266 s.set_max_range(max_range);
7267 s.solve().expect("solve")
7268 }
7269
7270 fn yz_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
7272 let pts = &result.points;
7273 for i in 1..pts.len() {
7274 if pts[i].position.x >= x {
7275 let (p1, p2) = (&pts[i - 1], &pts[i]);
7276 let dx = p2.position.x - p1.position.x;
7277 let t = if dx.abs() < 1e-12 { 0.0 } else { (x - p1.position.x) / dx };
7278 return (
7279 p1.position.y + t * (p2.position.y - p1.position.y),
7280 p1.position.z + t * (p2.position.z - p1.position.z),
7281 );
7282 }
7283 }
7284 panic!("trajectory never reached {x} m");
7285 }
7286
7287 #[test]
7288 fn cant_sign_clockwise_up_offset_goes_right_and_low() {
7289 let mut level = base_inputs();
7291 level.muzzle_angle = 0.003; let mut canted = level.clone();
7293 canted.cant_angle = 10f64.to_radians();
7294
7295 let (y0, z0) = yz_at(&solve_with(level, 400.0), 300.0);
7296 let (y1, z1) = yz_at(&solve_with(canted, 400.0), 300.0);
7297 assert!(z1 > z0 + 0.01, "clockwise cant must move POI right: z0={z0} z1={z1}");
7298 assert!(y1 < y0 - 0.001, "clockwise cant must move POI low: y0={y0} y1={y1}");
7299 }
7300
7301 #[test]
7302 fn pure_cant_shows_bore_offset_near_range() {
7303 let mut i = base_inputs();
7306 i.muzzle_angle = 0.0;
7307 i.cant_angle = 10f64.to_radians();
7308 let sh = i.sight_height;
7309 let r = solve_with(i, 60.0);
7310 let first = &r.points[1]; let expected = -sh * 10f64.to_radians().sin();
7312 assert!(
7313 (first.position.z - expected).abs() < 0.005,
7314 "near-muzzle lateral {} should be ~bore offset {expected}",
7315 first.position.z
7316 );
7317 }
7318
7319 #[test]
7320 fn zero_angle_is_independent_of_cant() {
7321 let a = base_inputs();
7322 let mut b = base_inputs();
7323 b.cant_angle = 15f64.to_radians();
7324 let za = calculate_zero_angle(a.clone(), 100.0, 0.0).expect("zero a");
7325 let zb = calculate_zero_angle(b.clone(), 100.0, 0.0).expect("zero b");
7326 assert_eq!(za.to_bits(), zb.to_bits(), "zeroing must ignore cant: {za} vs {zb}");
7327 let _ = (a.cant_angle, b.cant_angle);
7329 }
7330
7331 #[test]
7332 fn nonfinite_cant_is_rejected() {
7333 let mut i = base_inputs();
7334 i.cant_angle = f64::NAN;
7335 let s = TrajectorySolver::new(i, WindConditions::default(), AtmosphericConditions::default());
7336 assert!(s.solve().is_err());
7337 }
7338
7339 #[test]
7340 fn incline_and_cant_compose_without_breaking() {
7341 let mut flat = base_inputs();
7343 flat.muzzle_angle = 0.003;
7344 flat.shooting_angle = 15f64.to_radians();
7345 let mut canted = flat.clone();
7346 canted.cant_angle = 10f64.to_radians();
7347 let (_, z_flat) = yz_at(&solve_with(flat, 400.0), 300.0);
7348 let (_, z_cant) = yz_at(&solve_with(canted, 400.0), 300.0);
7349 assert!(z_cant > z_flat, "cant must still deflect right on an incline");
7350 }
7351}
7352
7353#[cfg(test)]
7354mod vertical_wind_tests {
7355 use super::*;
7356
7357 fn base_inputs() -> BallisticInputs {
7358 BallisticInputs {
7359 muzzle_velocity: 800.0,
7360 bc_value: 0.5,
7361 bc_type: DragModel::G7,
7362 bullet_mass: 0.0109,
7363 bullet_diameter: 0.00782,
7364 bullet_length: 0.0309,
7365 sight_height: 0.05,
7366 twist_rate: 10.0,
7367 use_rk4: true,
7368 ..BallisticInputs::default()
7369 }
7370 }
7371
7372 fn y_at(result: &TrajectoryResult, x: f64) -> f64 {
7374 let pts = &result.points;
7375 for i in 1..pts.len() {
7376 if pts[i].position.x >= x {
7377 let (p1, p2) = (&pts[i - 1], &pts[i]);
7378 let dx = p2.position.x - p1.position.x;
7379 let t = if dx.abs() < 1e-12 { 0.0 } else { (x - p1.position.x) / dx };
7380 return p1.position.y + t * (p2.position.y - p1.position.y);
7381 }
7382 }
7383 panic!("trajectory never reached {x} m");
7384 }
7385
7386 fn solve_with(inputs: BallisticInputs, wind: WindConditions, max_range: f64) -> TrajectoryResult {
7387 let mut s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
7388 s.set_max_range(max_range);
7389 s.solve().expect("solve")
7390 }
7391
7392 #[test]
7393 fn updraft_raises_poi_downrange() {
7394 let calm_inputs = base_inputs();
7397 let calm_wind = WindConditions::default();
7398 let updraft = WindConditions {
7399 vertical_speed: 5.0,
7400 ..Default::default()
7401 };
7402
7403 let calm = solve_with(calm_inputs.clone(), calm_wind, 500.0);
7404 let updraft_result = solve_with(calm_inputs, updraft, 500.0);
7405
7406 let y_calm = y_at(&calm, 400.0);
7407 let y_updraft = y_at(&updraft_result, 400.0);
7408 assert!(
7409 y_updraft > y_calm,
7410 "5 m/s updraft must raise POI at 400m: calm={y_calm}, updraft={y_updraft}"
7411 );
7412 }
7413
7414 #[test]
7415 fn zero_vertical_is_default_and_finite_required() {
7416 assert_eq!(WindConditions::default().vertical_speed, 0.0);
7417
7418 let inputs = base_inputs();
7419 let wind = WindConditions {
7420 vertical_speed: f64::NAN,
7421 ..Default::default()
7422 };
7423 let s = TrajectorySolver::new(inputs, wind, AtmosphericConditions::default());
7424 assert!(
7425 s.solve().is_err(),
7426 "NaN wind.vertical_speed must be rejected by validate_for_solve"
7427 );
7428 }
7429}
7430
7431#[cfg(test)]
7433mod bc_reference_standard_tests {
7434 use super::*;
7435
7436 fn base_inputs() -> BallisticInputs {
7437 BallisticInputs {
7438 muzzle_velocity: 800.0,
7439 bc_value: 0.5,
7440 bc_type: DragModel::G7,
7441 bullet_mass: 0.0109,
7442 bullet_diameter: 0.00782,
7443 bullet_length: 0.0309,
7444 sight_height: 0.05,
7445 twist_rate: 10.0,
7446 use_rk4: true,
7447 ..BallisticInputs::default()
7448 }
7449 }
7450
7451 fn y_and_speed_at(result: &TrajectoryResult, x: f64) -> (f64, f64) {
7453 let pts = &result.points;
7454 for i in 1..pts.len() {
7455 if pts[i].position.x >= x {
7456 let (p1, p2) = (&pts[i - 1], &pts[i]);
7457 let dx = p2.position.x - p1.position.x;
7458 let t = if dx.abs() < 1e-12 {
7459 0.0
7460 } else {
7461 (x - p1.position.x) / dx
7462 };
7463 return (
7464 p1.position.y + t * (p2.position.y - p1.position.y),
7465 p1.velocity_magnitude + t * (p2.velocity_magnitude - p1.velocity_magnitude),
7466 );
7467 }
7468 }
7469 panic!("trajectory never reached {x} m");
7470 }
7471
7472 #[test]
7475 fn asm_to_icao_ratio_matches_documented_value() {
7476 assert!(
7477 (crate::constants::ASM_TO_ICAO_BC - 0.98237).abs() < 1e-5,
7478 "ASM_TO_ICAO_BC = {} must equal 0.98237 to 5 decimal places",
7479 crate::constants::ASM_TO_ICAO_BC
7480 );
7481 assert_eq!(
7486 crate::constants::ASM_TO_ICAO_BC,
7487 crate::constants::ASM_DENSITY_LB_FT3 / crate::constants::ICAO_DENSITY_LB_FT3
7488 );
7489 }
7490
7491 #[test]
7494 fn default_bc_reference_standard_is_icao() {
7495 assert_eq!(
7496 BallisticInputs::default().bc_reference_standard,
7497 BcReferenceStandard::Icao
7498 );
7499 }
7500
7501 #[test]
7506 fn icao_reference_leaves_bc_value_bit_identical() {
7507 let raw_bc: f64 = 0.4372911; let inputs = BallisticInputs {
7509 bc_value: raw_bc,
7510 bc_reference_standard: BcReferenceStandard::Icao,
7511 ..base_inputs()
7512 };
7513 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
7514 assert_eq!(solver.inputs.bc_value.to_bits(), raw_bc.to_bits());
7515 }
7516
7517 #[test]
7518 fn default_inputs_solve_is_unaffected_by_the_new_field_existing() {
7519 let a = TrajectorySolver::new(base_inputs(), WindConditions::default(), AtmosphericConditions::default())
7523 .solve()
7524 .expect("solve a");
7525 let b = TrajectorySolver::new(
7526 BallisticInputs { ..base_inputs() },
7527 WindConditions::default(),
7528 AtmosphericConditions::default(),
7529 )
7530 .solve()
7531 .expect("solve b");
7532 assert_eq!(a.impact_velocity.to_bits(), b.impact_velocity.to_bits());
7533 assert_eq!(a.max_range.to_bits(), b.max_range.to_bits());
7534 }
7535
7536 #[test]
7539 fn army_standard_metro_scales_bc_value_by_exactly_the_derived_ratio() {
7540 let raw_bc = 0.5;
7541 let inputs = BallisticInputs {
7542 bc_value: raw_bc,
7543 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
7544 ..base_inputs()
7545 };
7546 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
7547 assert_eq!(
7548 solver.inputs.bc_value,
7549 raw_bc * crate::constants::ASM_TO_ICAO_BC
7550 );
7551 }
7552
7553 #[test]
7554 fn army_standard_metro_scales_mach_keyed_bc_segments() {
7555 let inputs = BallisticInputs {
7556 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
7557 bc_segments: Some(vec![(0.5, 0.40), (1.5, 0.30)]),
7558 ..base_inputs()
7559 };
7560 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
7561 let segments = solver.inputs.bc_segments.as_ref().expect("segments");
7562 assert_eq!(segments[0], (0.5, 0.40 * crate::constants::ASM_TO_ICAO_BC));
7563 assert_eq!(segments[1], (1.5, 0.30 * crate::constants::ASM_TO_ICAO_BC));
7564 }
7565
7566 #[test]
7567 fn army_standard_metro_scales_velocity_keyed_bc_segments_data() {
7568 let inputs = BallisticInputs {
7569 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
7570 bc_segments_data: Some(vec![
7571 crate::BCSegmentData {
7572 velocity_min: 0.0,
7573 velocity_max: 500.0,
7574 bc_value: 0.40,
7575 },
7576 crate::BCSegmentData {
7577 velocity_min: 500.0,
7578 velocity_max: 900.0,
7579 bc_value: 0.45,
7580 },
7581 ]),
7582 ..base_inputs()
7583 };
7584 let solver = TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
7585 let segments = solver.inputs.bc_segments_data.as_ref().expect("segments");
7586 assert_eq!(segments[0].bc_value, 0.40 * crate::constants::ASM_TO_ICAO_BC);
7587 assert_eq!(segments[1].bc_value, 0.45 * crate::constants::ASM_TO_ICAO_BC);
7588 assert_eq!(segments[0].velocity_min, 0.0);
7590 assert_eq!(segments[1].velocity_max, 900.0);
7591 }
7592
7593 #[test]
7598 fn army_standard_metro_moves_impact_in_the_more_drag_direction() {
7599 let solve_at = |standard: BcReferenceStandard| {
7600 let inputs = BallisticInputs {
7601 bc_value: 0.475,
7602 bc_reference_standard: standard,
7603 ..base_inputs()
7604 };
7605 let mut solver =
7606 TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
7607 solver.set_max_range(500.0);
7608 solver.solve().expect("solve")
7609 };
7610
7611 let icao = solve_at(BcReferenceStandard::Icao);
7612 let asm = solve_at(BcReferenceStandard::ArmyStandardMetro);
7613
7614 let (y_icao, v_icao) = y_and_speed_at(&icao, 400.0);
7615 let (y_asm, v_asm) = y_and_speed_at(&asm, 400.0);
7616
7617 assert!(
7618 y_asm < y_icao,
7619 "ArmyStandardMetro must drop MORE (lower y) at 400m than Icao for the same raw \
7620 bc_value: icao_y={y_icao}, asm_y={y_asm}"
7621 );
7622 assert!(
7623 v_asm < v_icao,
7624 "ArmyStandardMetro must retain LESS velocity at 400m than Icao for the same raw \
7625 bc_value: icao_v={v_icao}, asm_v={v_asm}"
7626 );
7627 }
7628
7629 #[test]
7636 fn monte_carlo_inherits_the_normalized_bc_reference() {
7637 let base_inputs_asm = BallisticInputs {
7638 bc_value: 0.475,
7639 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
7640 ..base_inputs()
7641 };
7642 let wind = WindConditions::default();
7643
7644 let mut direct_solver =
7650 TrajectorySolver::new(base_inputs_asm.clone(), wind.clone(), AtmosphericConditions::default());
7651 direct_solver.set_max_range(base_inputs_asm.target_distance.max(1000.0) * 2.0);
7652 let direct = direct_solver.solve().expect("direct solve");
7653
7654 let mc_params = MonteCarloParams {
7655 num_simulations: 1,
7656 velocity_std_dev: 0.0,
7657 angle_std_dev: 0.0,
7658 bc_std_dev: 0.0,
7659 wind_speed_std_dev: 0.0,
7660 target_distance: None,
7661 base_wind_speed: 0.0,
7662 base_wind_direction: 0.0,
7663 azimuth_std_dev: 0.0,
7664 };
7665 let mc = run_monte_carlo_with_wind_and_direction_std_dev_seeded(
7666 base_inputs_asm,
7667 wind,
7668 mc_params,
7669 0.0,
7670 42,
7671 )
7672 .expect("monte carlo");
7673
7674 assert_eq!(mc.ranges.len(), 1);
7675 assert_eq!(
7676 mc.ranges[0].to_bits(),
7677 direct.max_range.to_bits(),
7678 "a zero-dispersion single MC sample must match a plain solve of the same \
7679 ASM-referenced inputs bit-for-bit"
7680 );
7681 assert_eq!(
7682 mc.impact_velocities[0].to_bits(),
7683 direct.impact_velocity.to_bits()
7684 );
7685 }
7686
7687 #[test]
7695 fn estimate_bc_fit_recovers_an_icao_referenced_bc() {
7696 let known_bc = 0.475;
7697 let velocity = 800.0;
7698 let mass = 0.0109;
7699 let diameter = 0.00782;
7700 let atmosphere = AtmosphericConditions::default();
7701
7702 let synth_inputs = BallisticInputs {
7703 muzzle_velocity: velocity,
7704 bc_value: known_bc,
7705 bc_type: DragModel::G7,
7706 bullet_mass: mass,
7707 bullet_diameter: diameter,
7708 bullet_length: 0.0309,
7709 sight_height: 0.05,
7710 twist_rate: 10.0,
7711 use_rk4: true,
7712 bc_reference_standard: BcReferenceStandard::Icao,
7713 ..BallisticInputs::default()
7714 };
7715 let mut solver = TrajectorySolver::new(synth_inputs, WindConditions::default(), atmosphere.clone());
7716 solver.set_max_range(500.0);
7717 let trajectory = solver.solve().expect("synthetic solve");
7718
7719 let points: Vec<(f64, f64)> = [100.0, 200.0, 300.0, 400.0]
7720 .iter()
7721 .map(|&d| {
7722 let (y, _) = {
7723 let pts = &trajectory.points;
7724 let mut found = None;
7725 for i in 1..pts.len() {
7726 if pts[i].position.x >= d {
7727 let (p1, p2) = (&pts[i - 1], &pts[i]);
7728 let dx = p2.position.x - p1.position.x;
7729 let t = if dx.abs() < 1e-12 {
7730 0.0
7731 } else {
7732 (d - p1.position.x) / dx
7733 };
7734 found = Some((
7735 p1.position.y + t * (p2.position.y - p1.position.y),
7736 0.0,
7737 ));
7738 break;
7739 }
7740 }
7741 found.expect("trajectory reached observation distance")
7742 };
7743 (d, -y) })
7745 .collect();
7746
7747 let estimate = estimate_bc_fit(
7748 velocity,
7749 mass,
7750 diameter,
7751 &points,
7752 DragModel::G7,
7753 BcFitMode::Drop,
7754 atmosphere,
7755 None,
7756 0.05,
7757 )
7758 .expect("fit should converge");
7759
7760 assert!(
7761 (estimate.bc - known_bc).abs() < 0.02,
7762 "fit should recover the known ICAO-referenced bc={known_bc}, got {}",
7763 estimate.bc
7764 );
7765 }
7766
7767 #[test]
7770 fn custom_drag_table_makes_bc_reference_standard_numerically_inert() {
7771 let table = crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0, 3.0], vec![0.3, 0.4, 0.3, 0.2])
7772 .expect("valid table");
7773
7774 let solve_with = |standard: BcReferenceStandard| {
7775 let inputs = BallisticInputs {
7776 bc_value: 0.5, bc_reference_standard: standard,
7778 custom_drag_table: Some(table.clone()),
7779 ..base_inputs()
7780 };
7781 let mut solver =
7782 TrajectorySolver::new(inputs, WindConditions::default(), AtmosphericConditions::default());
7783 solver.set_max_range(500.0);
7784 solver.solve().expect("solve")
7785 };
7786
7787 let icao = solve_with(BcReferenceStandard::Icao);
7788 let asm = solve_with(BcReferenceStandard::ArmyStandardMetro);
7789
7790 assert_eq!(
7791 icao.impact_velocity.to_bits(),
7792 asm.impact_velocity.to_bits(),
7793 "a custom drag table must make bc_reference_standard fully inert"
7794 );
7795 assert_eq!(icao.max_range.to_bits(), asm.max_range.to_bits());
7796 }
7797
7798 #[test]
7799 fn custom_drag_table_inert_warning_fires_only_for_army_standard_metro_with_a_table() {
7800 let table = crate::drag::DragTable::try_new(vec![0.5, 1.0, 2.0], vec![0.3, 0.4, 0.3])
7801 .expect("valid table");
7802
7803 let no_table_icao = base_inputs();
7805 assert!(no_table_icao.bc_reference_standard_inert_warning().is_none());
7806 let no_table_asm = BallisticInputs {
7807 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
7808 ..base_inputs()
7809 };
7810 assert!(no_table_asm.bc_reference_standard_inert_warning().is_none());
7811
7812 let table_icao = BallisticInputs {
7814 custom_drag_table: Some(table.clone()),
7815 ..base_inputs()
7816 };
7817 assert!(table_icao.bc_reference_standard_inert_warning().is_none());
7818
7819 let table_asm = BallisticInputs {
7821 custom_drag_table: Some(table),
7822 bc_reference_standard: BcReferenceStandard::ArmyStandardMetro,
7823 ..base_inputs()
7824 };
7825 let warning = table_asm
7826 .bc_reference_standard_inert_warning()
7827 .expect("must warn");
7828 assert!(warning.contains("--bc-reference"));
7829 assert!(warning.contains("--drag-table"));
7830 }
7831}
7832
7833#[cfg(test)]
7839mod effective_drag_coefficient_tests {
7840 use super::*;
7841
7842 fn inputs_175gr_g7() -> BallisticInputs {
7843 let mut inputs = BallisticInputs {
7844 bc_value: 0.243,
7845 bc_type: DragModel::G7,
7846 muzzle_velocity: 823.0,
7847 ..Default::default()
7848 };
7849 inputs.bullet_mass = 175.0 * crate::constants::GRAINS_TO_KG;
7853 inputs.bullet_diameter = 0.308 * 0.0254;
7854 inputs.weight_grains = 175.0;
7855 inputs.caliber_inches = 0.308;
7856 inputs
7857 }
7858
7859 fn solver(inputs: BallisticInputs) -> TrajectorySolver {
7860 TrajectorySolver::new(
7861 inputs,
7862 WindConditions::default(),
7863 AtmosphericConditions::default(),
7864 )
7865 }
7866
7867 #[test]
7871 fn reports_the_projectiles_own_cd_not_the_reference_tables() {
7872 let inputs = inputs_175gr_g7();
7873 let sd = inputs.sectional_density_lb_in2().expect("SD");
7874 let solver = solver(inputs);
7875
7876 let sos = 340.0;
7877 let velocity = 800.0;
7878 let mach = velocity / sos;
7879
7880 let reference = crate::drag::get_drag_coefficient(mach, &DragModel::G7);
7881 let reported = solver
7882 .effective_drag_coefficient(velocity, sos)
7883 .expect("mass and diameter are set");
7884
7885 let expected = reference * sd / 0.243;
7886 assert!(
7887 (reported - expected).abs() < 1e-12,
7888 "reported {reported} != Cd_ref * SD / BC {expected}"
7889 );
7890 assert!(
7893 (reported - reference).abs() > 1e-6,
7894 "form factor collapsed to 1; this fixture no longer distinguishes the two values"
7895 );
7896 }
7897
7898 #[test]
7901 fn a_custom_drag_table_passes_through_unscaled() {
7902 let mut inputs = inputs_175gr_g7();
7903 inputs.custom_drag_table = Some(crate::drag::DragTable::new(
7904 vec![0.5, 3.0],
7905 vec![0.15, 0.40],
7906 ));
7907 let solver = solver(inputs);
7908
7909 let sos = 340.0;
7910 let velocity = 0.9 * sos;
7911 let table_value = solver
7912 .inputs
7913 .custom_drag_table
7914 .as_ref()
7915 .expect("table")
7916 .interpolate(0.9);
7917
7918 let reported = solver
7919 .effective_drag_coefficient(velocity, sos)
7920 .expect("mass and diameter are set");
7921 assert!(
7922 (reported - table_value).abs() < 1e-12,
7923 "custom table Cd {table_value} was rescaled to {reported}"
7924 );
7925 }
7926
7927 #[test]
7930 fn a_velocity_segmented_bc_steps_the_reported_cd() {
7931 let mut inputs = inputs_175gr_g7();
7932 inputs.use_bc_segments = true;
7933 inputs.bc_segments_data = Some(vec![
7934 crate::BCSegmentData { velocity_min: 2400.0, velocity_max: 4000.0, bc_value: 0.243 },
7935 crate::BCSegmentData { velocity_min: 0.0, velocity_max: 2400.0, bc_value: 0.200 },
7936 ]);
7937 let solver = solver(inputs);
7938
7939 let sos = 340.0;
7940 let above = solver.effective_drag_coefficient(2500.0 / 3.28084, sos).expect("cd");
7942 let below = solver.effective_drag_coefficient(2300.0 / 3.28084, sos).expect("cd");
7943
7944 assert!(
7946 below > above,
7947 "expected the 0.200 band to report a higher Cd than the 0.243 band; got {below} vs {above}"
7948 );
7949 }
7950
7951 #[test]
7954 fn is_absent_when_sectional_density_is_unknown() {
7955 let mut inputs = inputs_175gr_g7();
7956 inputs.weight_grains = 0.0;
7957 inputs.bullet_mass = 0.0;
7958 let solver = solver(inputs);
7959 assert!(solver.effective_drag_coefficient(800.0, 340.0).is_none());
7960 }
7961
7962 #[test]
7967 fn the_json_emit_rule_is_flag_gated_and_absent_when_cd_is_unknown() {
7968 let mut point = TrajectoryPoint {
7969 time: 0.0,
7970 position: nalgebra::Vector3::new(0.0, 0.0, 0.0),
7971 velocity_magnitude: 800.0,
7972 kinetic_energy: 3000.0,
7973 drag_coefficient: Some(0.31),
7974 };
7975 assert_eq!(point.drag_coefficient_json_value(true), Some(0.31));
7976 assert_eq!(
7977 point.drag_coefficient_json_value(false),
7978 None,
7979 "without the flag the key must not exist, so default JSON stays byte-identical"
7980 );
7981 point.drag_coefficient = None;
7982 assert_eq!(
7983 point.drag_coefficient_json_value(true),
7984 None,
7985 "unknown sectional density must yield an ABSENT key, not null"
7986 );
7987 }
7988
7989 #[test]
7991 fn every_point_of_a_solved_trajectory_carries_the_value() {
7992 let mut solver = solver(inputs_175gr_g7());
7993 solver.set_max_range(300.0);
7994 let result = solver.solve().expect("solve");
7995 assert!(!result.points.is_empty());
7996 assert!(
7997 result.points.iter().all(|p| p.drag_coefficient.is_some()),
7998 "the post-integration pass missed at least one point"
7999 );
8000 }
8001}