hpr_motor/
lib.rs

1//! Solid rocket motor model (thrust, propellant mass, CG and inertia over time), `.eng`/`.rse`
2//! reading and writing, and catalog types.
3//!
4//! **Guide:** [Solid motors][guide-motor], [`.eng` files][guide-eng] and [`.rse` files][guide-rse]:
5//! the models, their sources, how well they are validated and what they leave out.
6//!
7//! [guide-motor]: https://nrdptel.github.io/hpr-sim/physics/motor.html
8//! [guide-eng]: https://nrdptel.github.io/hpr-sim/format/eng.html
9//! [guide-rse]: https://nrdptel.github.io/hpr-sim/format/rse.html
10//!
11//! - [`curve`]: thrust curves, total impulse, NFPA 1125 burn time and average thrust.
12//! - [`class`]: impulse classes (`1/8A` to `O` and beyond).
13//! - [`motor`]: the solid motor: consumption, mass properties and ambient-pressure thrust, with
14//!   [`grains`] (BATES grains) and [`mass`] (axisymmetric mass elements).
15//! - [`eng`] and [`rse`]: RASP and RockSim motor files, with [`delay`] strings.
16//! - [`catalog`]: the offline catalog and its bundled ThrustCurve.org curves.
17
18mod bundled;
19pub mod catalog;
20pub mod class;
21pub mod convert;
22pub mod curve;
23pub mod delay;
24pub mod eng;
25pub mod error;
26pub mod grains;
27pub mod mass;
28pub mod motor;
29pub mod rse;
30pub mod text;
31
32pub use catalog::{Catalog, CatalogCurve, CatalogMotor};
33pub use class::ImpulseClass;
34pub use curve::ThrustCurve;
35pub use delay::{Delay, DelayList, DelayWarning};
36pub use error::MotorError;
37pub use grains::{BatesGrains, GrainShape};
38pub use mass::MassElement;
39pub use motor::{
40    EXHAUST_VELOCITY_RANGE_M_S, MotorState, Nozzle, Propellant, PropellantColumn, SolidMotor,
41};
42pub use text::{ParseWarning, Parsed, WarningKind};
43
44#[cfg(test)]
45mod tests {
46    use super::*;
47
48    /// Loft lesson L38: Loft's class letter was off by one at band tops (2.5 N·s gave `B`) and
49    /// had no `1/8A`.
50    #[test]
51    fn impulse_class_upper_bound_inclusive() {
52        let class = |ns: f64| ImpulseClass::from_total_impulse(ns).unwrap().label();
53        assert_eq!(class(2.5), "A");
54        assert_eq!(class(2.500_000_000_000_001), "B");
55        assert_eq!(class(5.0), "B");
56        assert_eq!(class(10.0), "C");
57        assert_eq!(class(160.0), "G");
58        assert_eq!(class(160.000_1), "H");
59        assert_eq!(class(40960.0), "O");
60        assert_eq!(class(0.3125), "1/8A");
61        assert_eq!(class(0.1), "1/8A");
62        assert_eq!(class(0.3126), "1/4A");
63        assert_eq!(class(1.25), "1/2A");
64        assert_eq!(class(1.2501), "A");
65    }
66
67    /// Loft lesson L40: Loft fixed the motor CG at the casing midpoint with no inertia of its own,
68    /// and its impulse-fraction model was uncited.
69    #[test]
70    fn cg_and_inertia_move_from_loaded_to_burnout() {
71        let curve =
72            ThrustCurve::new(vec![0.05, 0.2, 1.8, 2.0], vec![900.0, 800.0, 700.0, 0.0]).unwrap();
73        // A 38 mm reload: a heavy nozzle and aft closure put the dry center at 0.10 m, and three
74        // grains are centered further forward, at 0.20 m.
75        let grains = BatesGrains {
76            count: 3,
77            density_kg_m3: 1815.0,
78            outer_radius_m: 0.0165,
79            initial_inner_radius_m: 0.006,
80            initial_height_m: 0.09,
81            separation_m: 0.005,
82            center_m: 0.2,
83            inhibited_ends: false,
84        };
85        let dry = MassElement {
86            mass_kg: 0.45,
87            cg_m: 0.10,
88            axial_inertia_kg_m2: 1.6e-4,
89            transverse_inertia_kg_m2: 7.6e-3,
90        };
91        let motor = SolidMotor::new(curve, Propellant::Grains(grains), dry, None).unwrap();
92        let loaded = motor.state(0.0);
93        let burnout = motor.state(motor.burnout_time_s());
94
95        // Loaded: the parallel-axis combination of the dry mass and the full grains.
96        let m_p = grains.initial_mass_kg();
97        assert_eq!(loaded.propellant.mass_kg, m_p);
98        let cg = (0.45 * 0.10 + m_p * 0.2) / (0.45 + m_p);
99        assert!((loaded.total.cg_m - cg).abs() < 1e-15);
100        let full = grains.mass_element(m_p);
101        let transverse = dry.transverse_inertia_kg_m2
102            + 0.45 * (0.10 - cg).powi(2)
103            + full.transverse_inertia_kg_m2
104            + m_p * (0.2 - cg).powi(2);
105        assert!((loaded.total.transverse_inertia_kg_m2 - transverse).abs() < 1e-15);
106
107        // Burnout: only the dry mass is left, with its own inertia.
108        assert_eq!(burnout.propellant.mass_kg, 0.0);
109        assert_eq!(burnout.total.mass_kg, 0.45);
110        assert!((burnout.total.cg_m - 0.10).abs() < 1e-15);
111        assert!((burnout.total.axial_inertia_kg_m2 - 1.6e-4).abs() < 1e-18);
112        assert!((burnout.total.transverse_inertia_kg_m2 - 7.6e-3).abs() < 1e-15);
113
114        // In between, the center moves aft and both inertias fall, monotonically.
115        let states: Vec<MotorState> = (0..=200)
116            .map(|i| motor.state(f64::from(i) * 0.01))
117            .collect();
118        for pair in states.windows(2) {
119            assert!(pair[1].total.cg_m <= pair[0].total.cg_m + 1e-15);
120            assert!(pair[1].total.mass_kg <= pair[0].total.mass_kg);
121            assert!(pair[1].total.axial_inertia_kg_m2 <= pair[0].total.axial_inertia_kg_m2 + 1e-18);
122        }
123        assert!(loaded.total.cg_m - burnout.total.cg_m > 0.03);
124        assert!(
125            loaded.total.transverse_inertia_kg_m2 > 1.3 * burnout.total.transverse_inertia_kg_m2
126        );
127    }
128
129    /// Loft lesson L39: Loft took the last sample as the burn time. ThrustCurve's glossary uses
130    /// NFPA 1125: from 5% of peak thrust on the way up to 5% of peak on the way down.
131    #[test]
132    fn burn_time_uses_the_nfpa_1125_definition() {
133        // A 1 s ramp to a 200 N peak, a 1 s plateau at 100 N, and a 1 s tail-off to zero; then a
134        // long trickle at 2 N (1% of peak) that NFPA 1125 leaves out.
135        let curve = ThrustCurve::new(
136            vec![0.0, 1.0, 1.001, 2.0, 3.0, 3.001, 6.0, 6.001],
137            vec![0.0, 200.0, 100.0, 100.0, 0.0, 2.0, 2.0, 0.0],
138        )
139        .unwrap();
140        assert_eq!(curve.end_time_s(), 6.001);
141        // Up through 10 N at 0.05 s; down through 10 N at 2.9 s.
142        let (start, end) = curve.burn_window_s();
143        assert!((start - 0.05).abs() < 1e-12, "{start}");
144        assert!((end - 2.9).abs() < 1e-12, "{end}");
145        assert!((curve.burn_time_s() - 2.85).abs() < 1e-12);
146        assert!(
147            (curve.average_thrust_n() - curve.total_impulse_ns() / 2.85).abs() < 1e-9,
148            "average thrust is the total impulse over the NFPA burn time"
149        );
150    }
151}
hpr_motor/lib.rs

hpr_motor/
lib.rs
