rockets_engine/
ballistics.rs

1//! Adapter to ballistics-engine 0.48's public derivative kernel.
2//!
3//! Uses the body's own Cd-versus-Mach curve, fixed diameter, standard gravity and
4//! a level world frame. SI mass and its imperial mirror are updated at every stage,
5//! so aerodynamic deceleration responds to fuel consumption. Modified point-mass
6//! effects — spin/Magnus with yaw of repose, tip-off yaw and Coriolis — are
7//! opt-in via the `with_*` builders; reference-BC and postprocessing effects
8//! (Litz spin drift) stay disabled.
9
10use crate::{finite_vec, positive, Error, FlightState, ForceModel, Result, Vec3};
11use ballistics_engine::{
12    atmosphere::calculate_air_density_cimp, constants::GRAINS_TO_KG,
13    derivatives::compute_derivatives, drag::DragTable, BallisticInputs,
14};
15
16/// Atmospheric inputs supported by the upstream derivative kernel.
17#[derive(Debug, Clone, Copy)]
18pub enum Atmosphere {
19    /// Station measurements; the upstream model applies its local altitude lapse.
20    /// Humidity affects station density; this derivative path uses dry-air sound speed.
21    Standard {
22        /// Elevation of the weather station above sea level, in metres.
23        station_altitude_m: f64,
24        /// Station temperature in degrees Celsius (greater than absolute zero).
25        temperature_c: f64,
26        /// Positive station pressure in hectopascals.
27        pressure_hpa: f64,
28        /// Relative humidity as a percentage in [0, 100].
29        humidity_percent: f64,
30    },
31    /// Constant density and sound speed, independent of position.
32    /// Upstream's direct-mode sentinel requires 0 <= density < 2 and sound speed > 200.
33    Uniform {
34        /// Air density in kilograms per cubic metre; zero disables drag.
35        density_kg_m3: f64,
36        /// Sound speed in metres per second.
37        speed_of_sound_mps: f64,
38    },
39}
40
41impl Default for Atmosphere {
42    fn default() -> Self {
43        Self::Standard {
44            station_altitude_m: 0.0,
45            temperature_c: 15.0,
46            pressure_hpa: 1013.25,
47            humidity_percent: 0.0,
48        }
49    }
50}
51
52impl Atmosphere {
53    fn parameters(self) -> Result<(f64, f64, f64, f64)> {
54        match self {
55            Self::Standard {
56                station_altitude_m,
57                temperature_c,
58                pressure_hpa,
59                humidity_percent,
60            } => {
61                if !station_altitude_m.is_finite()
62                    || !temperature_c.is_finite()
63                    || temperature_c <= -273.15
64                    || !humidity_percent.is_finite()
65                    || !(0.0..=100.0).contains(&humidity_percent)
66                {
67                    return Err(Error::invalid("standard atmosphere", "requires finite altitude, temperature above absolute zero and humidity in [0, 100]"));
68                }
69                positive(pressure_hpa, "pressure_hpa")?;
70                let ratio =
71                    calculate_air_density_cimp(temperature_c, pressure_hpa, humidity_percent)
72                        / 1.225;
73                positive(ratio, "station density ratio")?;
74                Ok((station_altitude_m, temperature_c, pressure_hpa, ratio))
75            }
76            Self::Uniform {
77                density_kg_m3,
78                speed_of_sound_mps,
79            } => {
80                if !density_kg_m3.is_finite()
81                    || !(0.0..2.0).contains(&density_kg_m3)
82                    || !speed_of_sound_mps.is_finite()
83                    || speed_of_sound_mps <= 200.0
84                {
85                    return Err(Error::invalid(
86                        "uniform atmosphere",
87                        "upstream requires density in [0, 2) and sound speed > 200 m/s",
88                    ));
89                }
90                Ok((density_kg_m3, speed_of_sound_mps, 0.0, 0.0))
91            }
92        }
93    }
94}
95
96/// Gravity and aerodynamic drag delegated to ballistics-engine's force kernel.
97#[derive(Debug, Clone)]
98pub struct BallisticForces {
99    inputs: BallisticInputs,
100    atmos_params: (f64, f64, f64, f64),
101    wind_mps: Vec3,
102    omega_rps: Option<Vec3>,
103}
104
105/// Earth rotation rate, rad/s (McCoy-level precision, matching the upstream kernel).
106const OMEGA_EARTH_RPS: f64 = 7.2921159e-5;
107
108impl BallisticForces {
109    /// Build using the vehicle's actual Cd curve, not a G-model reference curve or BC.
110    /// At least two ascending, nonnegative Mach values and positive Cd values are required.
111    /// Wind is a world-frame air velocity vector (positive X is a tailwind).
112    pub fn new(
113        diameter_m: f64,
114        cd_vs_mach: Vec<(f64, f64)>,
115        atmosphere: Atmosphere,
116        wind_mps: Vec3,
117    ) -> Result<Self> {
118        positive(diameter_m, "diameter_m")?;
119        let diameter_in = diameter_m / 0.0254;
120        positive(diameter_in * diameter_in, "diameter area")?;
121        if !finite_vec(&wind_mps) {
122            return Err(Error::invalid("wind_mps", "must be finite"));
123        }
124        let (mach, cd) = cd_vs_mach.into_iter().unzip();
125        let drag = DragTable::try_new(mach, cd).map_err(Error::Model)?;
126        let inputs = BallisticInputs {
127            bullet_diameter: diameter_m,
128            caliber_inches: diameter_in,
129            bullet_length: diameter_m,
130            custom_drag_table: Some(drag),
131            custom_drag_reference: false,
132            altitude: 0.0,
133            shooting_angle: 0.0,
134            enable_advanced_effects: false,
135            enable_magnus: false,
136            enable_coriolis: false,
137            use_enhanced_spin_drift: false,
138            tipoff_yaw: 0.0,
139            ..BallisticInputs::default()
140        };
141        Ok(Self {
142            inputs,
143            atmos_params: atmosphere.parameters()?,
144            wind_mps,
145            omega_rps: None,
146        })
147    }
148
149    /// Apply a positive multiplier to the entire Cd curve. Default is 1.0.
150    pub fn with_drag_scale(mut self, scale: f64) -> Result<Self> {
151        positive(scale, "drag scale")?;
152        if self
153            .inputs
154            .custom_drag_table
155            .as_ref()
156            .unwrap()
157            .cd_values
158            .iter()
159            .any(|cd| !(cd * scale).is_finite())
160        {
161            return Err(Error::invalid("drag scale", "scaled Cd must remain finite"));
162        }
163        self.inputs.cd_scale = scale;
164        Ok(self)
165    }
166
167    /// Change uniform world-frame wind between public integration steps.
168    pub fn set_wind(&mut self, wind_mps: Vec3) -> Result<()> {
169        if !finite_vec(&wind_mps) {
170            return Err(Error::invalid("wind_mps", "must be finite"));
171        }
172        self.wind_mps = wind_mps;
173        Ok(())
174    }
175
176    /// Modified point-mass spin: enable the kernel's Magnus force and yaw of
177    /// repose for a spinning body. `twist_rate_in` is the effective twist in
178    /// inches per turn — rifling language upstream; for a fin-canted rocket it
179    /// is the roll rate the cant produces per unit velocity. `body_length_m` is
180    /// the full airframe length (the stability correlations need it; [`Self::new`]
181    /// only knows the diameter). `right_hand` selects the roll direction.
182    ///
183    /// The correlations are McCoy bullet shapes, so treat magnitudes as
184    /// approximate for blunt or finned bodies; signs and scaling are physical.
185    pub fn with_spin(
186        mut self,
187        twist_rate_in: f64,
188        body_length_m: f64,
189        right_hand: bool,
190    ) -> Result<Self> {
191        positive(twist_rate_in, "twist_rate_in")?;
192        positive(body_length_m, "body_length_m")?;
193        self.inputs.twist_rate = twist_rate_in;
194        self.inputs.is_twist_right = right_hand;
195        self.inputs.bullet_length = body_length_m;
196        self.inputs.enable_magnus = true;
197        Ok(self)
198    }
199
200    /// Modified point-mass Coriolis deflection, for planet-scale arcs.
201    /// `latitude_deg` is the launch latitude in [-90, 90]; `shot_azimuth_deg` is
202    /// the compass bearing fired along, 0 = north, 90 = east. Do not combine
203    /// with an inertial-frame [`crate::CentralGravity`] model: Coriolis is the
204    /// rotating-frame term and would double-count Earth rotation.
205    pub fn with_coriolis(mut self, latitude_deg: f64, shot_azimuth_deg: f64) -> Result<Self> {
206        if !latitude_deg.is_finite() || !(-90.0..=90.0).contains(&latitude_deg) {
207            return Err(Error::invalid(
208                "latitude_deg",
209                "must be finite and in [-90, 90]",
210            ));
211        }
212        if !shot_azimuth_deg.is_finite() {
213            return Err(Error::invalid("shot_azimuth_deg", "must be finite"));
214        }
215        let lat = latitude_deg.to_radians();
216        let az = shot_azimuth_deg.to_radians();
217        self.omega_rps = Some(Vec3::new(
218            OMEGA_EARTH_RPS * lat.cos() * az.cos(),  // X: downrange
219            OMEGA_EARTH_RPS * lat.sin(),             // Y: vertical
220            -OMEGA_EARTH_RPS * lat.cos() * az.sin(), // Z: lateral
221        ));
222        self.inputs.latitude = Some(latitude_deg);
223        self.inputs.shot_azimuth = az;
224        self.inputs.enable_coriolis = true;
225        Ok(self)
226    }
227
228    /// Launch tip-off: an initial yaw in radians that decays exponentially with
229    /// downrange distance over `decay_distance_m`. Models the rail/launcher
230    /// disturbance; it perturbs the effective angle of attack early in flight.
231    pub fn with_tipoff_yaw(mut self, yaw_rad: f64, decay_distance_m: f64) -> Result<Self> {
232        if !yaw_rad.is_finite() {
233            return Err(Error::invalid("tipoff yaw_rad", "must be finite"));
234        }
235        positive(decay_distance_m, "tipoff decay_distance_m")?;
236        self.inputs.tipoff_yaw = yaw_rad;
237        self.inputs.tipoff_decay_distance = decay_distance_m;
238        Ok(self)
239    }
240
241    /// Evaluate using externally sampled wind, for a spatial or time-varying force wrapper.
242    pub fn acceleration_in_wind(
243        &self,
244        state: &FlightState,
245        mass_kg: f64,
246        wind_mps: Vec3,
247    ) -> Result<Vec3> {
248        state.validate()?;
249        positive(mass_kg, "mass_kg")?;
250        if !finite_vec(&wind_mps) {
251            return Err(Error::invalid("wind_mps", "must be finite"));
252        }
253        let mut inputs = self.inputs.clone();
254        inputs.bullet_mass = mass_kg;
255        inputs.weight_grains = mass_kg / GRAINS_TO_KG;
256        positive(inputs.weight_grains, "mass in grains")?;
257        let derivatives = compute_derivatives(
258            state.position_m,
259            state.velocity_mps,
260            &inputs,
261            wind_mps,
262            self.atmos_params,
263            inputs.bc_value,
264            self.omega_rps,
265            state.time_s,
266            None,
267        );
268        let acceleration = Vec3::new(derivatives[3], derivatives[4], derivatives[5]);
269        if !finite_vec(&acceleration) {
270            return Err(Error::NumericalFailure(
271                "ballistics force model returned non-finite acceleration",
272            ));
273        }
274        Ok(acceleration)
275    }
276}
277
278impl ForceModel for BallisticForces {
279    fn acceleration(&self, state: &FlightState, mass_kg: f64) -> Result<Vec3> {
280        self.acceleration_in_wind(state, mass_kg, self.wind_mps)
281    }
282}
rockets_engine/ballistics.rs

rockets_engine/
ballistics.rs
