1use std::f64::consts::PI;
39use std::sync::{Arc, OnceLock};
40
41use hpr_design::{Layout, NoseShape, Part, PlacedComponent};
42use serde::{Deserialize, Serialize};
43
44use crate::afterbody::SEPARATION_ONSET_RAD;
45use crate::body::{BodyGeometry, sinc};
46use crate::crossflow::BodyLift;
47use crate::custom::{DragModel, DragQuery, SharedDragModel};
48use crate::drag::{
49 BUILDUP_MACH_LIMIT, ComponentDrag, ComponentDragTerms, Drag, DragConditions, MOTOR_POD_SETS,
50 axial_drag_alpha_factor, body_friction_form_factor, couple_afterbody,
51};
52use crate::error::{AeroError, check_dimension, check_mach};
53use crate::fins::{
54 FinAero, FinLoading, FinRollTerms, fin_count_factor, interference_factor,
55 roll_damping_interference, roll_forcing_interference, roll_sum, side_sum,
56};
57use crate::shock_expansion::{
58 BodySegment, DEFAULT_ELEMENTS_PER_CURVE, SegmentSlope, ShockExpansionBody,
59 flare_corner_limit_rad,
60};
61use crate::supersonic_boattail::{wp_center_fraction, wp_slope};
62use crate::table::{DragTable, NormalForceLookup, NormalForceTable, TableReference};
63use crate::tube_fins::{TubeFinSetAero, check_tube_fin_mach};
64
65pub const MAX_CANT_RAD: f64 = 15.0 * std::f64::consts::PI / 180.0;
68
69pub const NORMAL_FORCE_MACH_LIMIT: f64 = 5.0;
72
73pub const SUPERSONIC_JOIN_START_MACH: f64 = 1.2;
77
78const SUPERSONIC_STEPS_PER_MACH: f64 = 20.0;
80
81const SUPERSONIC_FIRST_STEP: usize = 24;
83
84const SUPERSONIC_LAST_STEP: usize = 100;
86
87const SUPERSONIC_JOIN_STEPS: usize = 6;
89
90const SUPERSONIC_JOIN_BISECTIONS: usize = 64;
97
98const SUPERSONIC_SHELTER_FLOOR: f64 = 1e-6;
102
103pub const SUPERSONIC_JOIN_WIDTH_MACH: f64 =
106 SUPERSONIC_JOIN_STEPS as f64 / SUPERSONIC_STEPS_PER_MACH;
107
108#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
115#[serde(rename_all = "snake_case", tag = "switch")]
116#[non_exhaustive]
117pub enum SupersonicFallback {
118 RadiusStep {
121 component: String,
123 },
124 LongLip {
127 component: String,
129 },
130 SteepTip,
133 Other,
135}
136
137#[derive(Debug, Clone, PartialEq, Serialize)]
170#[non_exhaustive]
171pub struct SupersonicBody {
172 pub covered: usize,
181 pub join_start_mach: f64,
184 first_step: usize,
186 rows: Vec<Vec<SegmentSlope>>,
189 stationed: Vec<bool>,
192 lead: Option<Vec<SegmentSlope>>,
195 pub shape_weight: f64,
205}
206
207#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, Deserialize)]
210#[serde(rename_all = "snake_case")]
211#[non_exhaustive]
212pub enum SupersonicBoattail {
213 #[default]
218 WashingtonPettis,
219 Footnote8,
223}
224
225#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, Deserialize)]
228#[serde(rename_all = "snake_case")]
229#[non_exhaustive]
230pub enum SupersonicFlare {
231 #[default]
239 Marched,
240 SlenderBody,
243}
244
245#[derive(Debug, Clone, Copy, Default, PartialEq, Serialize, Deserialize)]
252#[serde(default, deny_unknown_fields)]
253#[non_exhaustive]
254pub struct BodyModel {
255 pub body_lift: BodyLift,
257 pub supersonic_boattail: SupersonicBoattail,
259 pub supersonic_flare: SupersonicFlare,
261}
262
263impl BodyModel {
264 pub const CURRENT: Self = Self {
267 body_lift: BodyLift::JORGENSEN,
268 supersonic_boattail: SupersonicBoattail::WashingtonPettis,
269 supersonic_flare: SupersonicFlare::Marched,
270 };
271
272 pub const BEFORE_M1_8E6: Self = Self {
276 body_lift: BodyLift::GALEJS,
277 supersonic_boattail: SupersonicBoattail::Footnote8,
278 supersonic_flare: SupersonicFlare::SlenderBody,
279 };
280
281 #[must_use]
283 pub const fn with_body_lift(mut self, body_lift: BodyLift) -> Self {
284 self.body_lift = body_lift;
285 self
286 }
287
288 #[must_use]
290 pub const fn with_supersonic_boattail(
291 mut self,
292 supersonic_boattail: SupersonicBoattail,
293 ) -> Self {
294 self.supersonic_boattail = supersonic_boattail;
295 self
296 }
297
298 #[must_use]
300 pub const fn with_supersonic_flare(mut self, supersonic_flare: SupersonicFlare) -> Self {
301 self.supersonic_flare = supersonic_flare;
302 self
303 }
304}
305
306#[derive(Debug, Clone, PartialEq)]
310struct RunBoattail {
311 in_its_place: Vec<BodySegment>,
312 fore_radius_m: f64,
313 aft_radius_m: f64,
314 length_m: f64,
315}
316
317#[derive(Debug, Clone, PartialEq)]
322struct RunFlare {
323 index: usize,
325 ahead: Vec<BodySegment>,
326 fore_radius_m: f64,
327 aft_radius_m: f64,
328 length_m: f64,
329 clipped: bool,
330}
331
332#[derive(Debug, Clone, PartialEq)]
337struct SupersonicRun {
338 segments: Vec<BodySegment>,
339 vertex_m: f64,
340 bounds_m: Vec<Option<(f64, f64)>>,
341 fore_m: Vec<f64>,
343 boattails: Vec<Option<RunBoattail>>,
345 flare: Option<RunFlare>,
347 sheltered_lips: usize,
351 shape_weight: f64,
356 }
360
361impl SupersonicRun {
362 fn shares(
371 &self,
372 body: &ShockExpansionBody,
373 in_its_place: &[Option<ShockExpansionBody>],
374 ahead: Option<&ShockExpansionBody>,
375 mach: f64,
376 reference_area_m2: f64,
377 ) -> Option<Vec<SegmentSlope>> {
378 let vertex_m = self.vertex_m;
379 let held = match (&self.flare, ahead) {
384 (Some(_), None) => return None,
387 (Some(flare), Some(ahead)) => {
388 let aft = ahead.aft_flow(mach).ok()?;
389 let rise_m = flare.aft_radius_m - flare.fore_radius_m;
390 let angle_rad = (rise_m / flare.length_m).atan();
391 let turn_limit_rad = flare_corner_limit_rad(aft.surface_mach).ok()?;
395 let angle_limit_rad =
396 (turn_limit_rad + aft.angle_rad).min(crate::blunt_tip::CONE_TABLE_CAP_RAD);
397 if angle_rad > angle_limit_rad {
398 if angle_limit_rad <= 0.0 || !angle_limit_rad.is_finite() {
403 return None;
404 }
405 let length_m = rise_m / angle_limit_rad.tan();
406 let mut segments = flare.ahead.clone();
407 segments.push(BodySegment::Profile {
408 profile: hpr_design::Profile::transition(
409 NoseShape::Conical {},
410 length_m,
411 flare.fore_radius_m,
412 flare.aft_radius_m,
413 flare.clipped,
414 )
415 .ok()?,
416 });
417 Some((
418 ShockExpansionBody::new(&segments, DEFAULT_ELEMENTS_PER_CURVE).ok()?,
419 length_m,
420 ))
421 } else {
422 None
423 }
424 }
425 _ => None,
426 };
427 let mut shares = match &held {
428 Some((drawn, _)) => drawn.segment_slopes(mach, reference_area_m2).ok()?,
429 None => body.segment_slopes(mach, reference_area_m2).ok()?,
430 };
431 if let (Some(flare), Some((_, drawn_length_m))) = (&self.flare, &held) {
435 let index = flare.index;
436 let share = *shares.get(index)?;
437 if share.slope_per_rad > 0.0 {
438 let fore_m = self.fore_m[index] - vertex_m;
439 let along = (share.moment_slope_m / share.slope_per_rad - fore_m) / drawn_length_m;
440 shares[index].moment_slope_m =
441 share.slope_per_rad * (fore_m + along * flare.length_m);
442 }
443 }
444 shares.extend(std::iter::repeat_n(
446 SegmentSlope::default(),
447 self.sheltered_lips,
448 ));
449 for (index, (boattail, cylinder_body)) in
450 self.boattails.iter().zip(in_its_place).enumerate()
451 {
452 let (Some(boattail), Some(cylinder_body)) = (boattail, cylinder_body) else {
453 continue;
454 };
455 let cylinder = *cylinder_body
458 .segment_slopes(mach, reference_area_m2)
459 .ok()?
460 .last()?;
461 let fore_radius_m = boattail.fore_radius_m;
462 let drop_m = fore_radius_m - boattail.aft_radius_m;
472 let read = |length_m| wp_slope(mach, fore_radius_m, boattail.aft_radius_m, length_m);
473 let at_true_angle = read(boattail.length_m).ok()?;
474 let held = read(boattail.length_m.max(drop_m / SEPARATION_ONSET_RAD.tan())).ok()?;
475 let ratio = boattail.aft_radius_m / fore_radius_m;
485 let measured = held.max(at_true_angle.min(2.0 * (ratio * ratio - 1.0)));
486 let increment = measured * PI * fore_radius_m * fore_radius_m / reference_area_m2;
487 let center_m =
488 self.fore_m[index] + wp_center_fraction(mach) * boattail.length_m - vertex_m;
489 shares[index] = SegmentSlope {
490 slope_per_rad: cylinder.slope_per_rad + increment,
491 moment_slope_m: cylinder.moment_slope_m + increment * center_m,
492 };
493 }
494 let on_segment = shares.iter().zip(&self.bounds_m).all(|(s, bounds)| {
495 let Some((fore, aft)) = *bounds else {
496 return s.slope_per_rad.is_finite() && s.moment_slope_m.is_finite();
497 };
498 let station = (s.moment_slope_m + s.slope_per_rad * vertex_m) / s.slope_per_rad;
499 let slack = 1e-9 * (aft - vertex_m);
500 s.slope_per_rad > 0.0 && station >= fore - slack && station <= aft + slack
501 });
502 on_segment.then(|| {
503 shares
504 .into_iter()
505 .map(|s| SegmentSlope {
506 slope_per_rad: s.slope_per_rad,
507 moment_slope_m: s.moment_slope_m + s.slope_per_rad * vertex_m,
508 })
509 .collect()
510 })
511 }
512}
513
514#[derive(Debug, Clone, Default)]
517struct SupersonicTable(Arc<OnceLock<Option<SupersonicBody>>>);
518
519impl PartialEq for SupersonicTable {
520 fn eq(&self, _: &Self) -> bool {
521 true
522 }
523}
524
525impl SupersonicBody {
526 fn new(run: &SupersonicRun, reference_area_m2: f64) -> Option<Self> {
530 let body = ShockExpansionBody::new(&run.segments, DEFAULT_ELEMENTS_PER_CURVE).ok()?;
531 let ahead = match &run.flare {
532 Some(flare) => {
533 Some(ShockExpansionBody::new(&flare.ahead, DEFAULT_ELEMENTS_PER_CURVE).ok()?)
534 }
535 None => None,
536 };
537 let mut in_its_place = Vec::with_capacity(run.boattails.len());
538 for boattail in &run.boattails {
539 in_its_place.push(match boattail {
540 Some(b) => Some(
541 ShockExpansionBody::new(&b.in_its_place, DEFAULT_ELEMENTS_PER_CURVE).ok()?,
542 ),
543 None => None,
544 });
545 }
546 let at = |step: f64| {
547 run.shares(
548 &body,
549 &in_its_place,
550 ahead.as_ref(),
551 step / SUPERSONIC_STEPS_PER_MACH,
552 reference_area_m2,
553 )
554 };
555 let mut rows = Vec::new();
556 let mut first_step = SUPERSONIC_LAST_STEP + 1;
557 for step in (SUPERSONIC_FIRST_STEP..=SUPERSONIC_LAST_STEP).rev() {
559 let Some(row) = at(step as f64) else {
560 break;
561 };
562 rows.push(row);
563 first_step = step;
564 }
565 rows.reverse();
566 if first_step + SUPERSONIC_JOIN_STEPS > SUPERSONIC_LAST_STEP {
568 return None;
569 }
570 let mut lead = None;
575 let mut join_start_mach = first_step as f64 / SUPERSONIC_STEPS_PER_MACH;
576 if first_step > SUPERSONIC_FIRST_STEP {
577 let (mut low, mut high) = (first_step as f64 - 1.0, first_step as f64);
578 let mut held = None;
579 for _ in 0..SUPERSONIC_JOIN_BISECTIONS {
580 let mid = 0.5 * (low + high);
581 if mid <= low || mid >= high {
582 break;
583 }
584 match at(mid) {
585 Some(row) => {
586 high = mid;
587 held = Some(row);
588 }
589 None => low = mid,
590 }
591 }
592 if let Some(row) = held {
593 join_start_mach = high / SUPERSONIC_STEPS_PER_MACH;
594 lead = Some(row);
595 }
596 }
597 Some(Self {
598 covered: run.segments.len() + run.sheltered_lips,
599 join_start_mach,
600 first_step,
601 rows,
602 lead,
603 shape_weight: run.shape_weight,
604 stationed: run
605 .bounds_m
606 .iter()
607 .map(Option::is_some)
608 .chain(std::iter::repeat_n(false, run.sheltered_lips))
610 .collect(),
611 })
612 }
613
614 #[must_use]
619 pub fn weight(&self, mach: f64) -> f64 {
620 self.shape_weight
621 * ((mach - self.join_start_mach) / SUPERSONIC_JOIN_WIDTH_MACH).clamp(0.0, 1.0)
622 }
623
624 pub fn share(&self, index: usize, mach: f64) -> Option<(f64, f64)> {
629 let x = mach * SUPERSONIC_STEPS_PER_MACH - self.first_step as f64;
630 let (a, b, t) = match &self.lead {
631 Some(lead) if x < 0.0 => {
633 let lead_x =
634 self.join_start_mach * SUPERSONIC_STEPS_PER_MACH - self.first_step as f64;
635 let t = ((x - lead_x) / -lead_x).clamp(0.0, 1.0);
636 (lead.get(index)?, self.rows[0].get(index)?, t)
637 }
638 _ => {
639 let i = (x.floor().max(0.0) as usize).min(self.rows.len() - 2);
641 let t = (x - i as f64).clamp(0.0, 1.0);
642 (self.rows[i].get(index)?, self.rows[i + 1].get(index)?, t)
643 }
644 };
645 Some((
646 a.slope_per_rad + t * (b.slope_per_rad - a.slope_per_rad),
647 a.moment_slope_m + t * (b.moment_slope_m - a.moment_slope_m),
648 ))
649 }
650}
651
652#[derive(Debug, Clone, Copy, Default, PartialEq, Serialize, Deserialize)]
654#[non_exhaustive]
655pub struct Roll {
656 pub forcing: f64,
658 pub damping: f64,
660}
661
662#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
664#[serde(deny_unknown_fields)]
665#[non_exhaustive]
666pub struct Flow {
667 pub mach: f64,
669 pub alpha_rad: f64,
672 pub roll_rad: f64,
675}
676
677impl Flow {
678 pub fn new(mach: f64, alpha_rad: f64, roll_rad: f64) -> Self {
681 Self {
682 mach,
683 alpha_rad,
684 roll_rad,
685 }
686 }
687
688 pub fn axial(mach: f64) -> Self {
690 Self::new(mach, 0.0, 0.0)
691 }
692
693 pub fn validate(&self) -> Result<(), AeroError> {
700 check_mach(self.mach, NORMAL_FORCE_MACH_LIMIT, "the normal force")?;
701 self.validate_angles()
702 }
703
704 fn validate_for_buildup(&self) -> Result<(), AeroError> {
707 check_mach(self.mach, BUILDUP_MACH_LIMIT, "the drag buildup")?;
708 self.validate_angles()
709 }
710
711 fn validate_angles(&self) -> Result<(), AeroError> {
713 if !(0.0..=PI).contains(&self.alpha_rad) {
714 return Err(AeroError::Domain {
715 what: "angle of attack",
716 value: self.alpha_rad,
717 });
718 }
719 if !self.roll_rad.is_finite() {
720 return Err(AeroError::Domain {
721 what: "flow roll angle",
722 value: self.roll_rad,
723 });
724 }
725 Ok(())
726 }
727}
728
729#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
731#[non_exhaustive]
732pub struct NormalForce {
733 pub coefficient: f64,
735 pub slope_per_rad: f64,
737 pub moment_m: f64,
740 pub moment_slope_m: f64,
744 pub cp_station_m: Option<f64>,
747 pub side_coefficient: f64,
750 pub side_moment_m: f64,
753 #[serde(default, skip_serializing_if = "Option::is_none")]
757 pub table: Option<NormalForceLookup>,
758}
759
760#[derive(Clone, Copy, Default)]
763struct Term {
764 slope: f64,
765 moment: f64,
766 side: f64,
767 side_moment: f64,
768 scale: f64,
769}
770
771impl Term {
772 fn times(self, k: f64) -> Term {
774 Term {
775 slope: k * self.slope,
776 moment: k * self.moment,
777 side: k * self.side,
778 side_moment: k * self.side_moment,
779 scale: k * self.scale,
780 }
781 }
782
783 fn add(self, other: Term) -> Term {
784 Term {
785 slope: self.slope + other.slope,
786 moment: self.moment + other.moment,
787 side: self.side + other.side,
788 side_moment: self.side_moment + other.side_moment,
789 scale: self.scale + other.scale,
790 }
791 }
792}
793
794impl NormalForce {
795 fn new(term: Term, alpha_rad: f64) -> Self {
796 let cancelled = term.slope.abs() <= 1e-12 * term.scale;
797 Self {
798 coefficient: term.slope * alpha_rad,
799 slope_per_rad: term.slope,
800 moment_m: term.moment * alpha_rad,
801 moment_slope_m: term.moment,
802 cp_station_m: (term.slope != 0.0 && !cancelled).then(|| term.moment / term.slope),
803 side_coefficient: term.side * alpha_rad,
804 side_moment_m: term.side_moment * alpha_rad,
805 table: None,
806 }
807 }
808}
809
810#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
812#[non_exhaustive]
813pub struct ComponentNormalForce {
814 pub id: String,
816 pub normal_force: NormalForce,
818}
819
820#[derive(Debug, Clone, PartialEq, Serialize)]
824#[non_exhaustive]
825pub struct BodyAero {
826 pub id: String,
828 pub fore_station_m: f64,
830 pub geometry: BodyGeometry,
832 pub step_area_m2: f64,
835 pub slope_per_rad: f64,
839 pub moment_slope_m: f64,
842 pub planform_ratio: f64,
845 pub lift_station_m: f64,
847}
848
849#[derive(Debug, Clone, PartialEq, Serialize)]
853#[non_exhaustive]
854pub struct FinSetAero {
855 pub id: String,
857 pub count: u32,
859 pub base_angle_rad: f64,
861 pub fin: FinAero,
863 pub count_factor: f64,
865 pub interference: f64,
867 pub fore_station_m: f64,
869 pub cant_rad: f64,
871 pub body_radius_m: f64,
873 pub roll_forcing_interference: f64,
875 pub roll_damping_interference: f64,
877 pub roll: FinRollTerms,
879 pub pods: Option<PodFins>,
881}
882
883#[derive(Debug, Clone, PartialEq, Serialize)]
887#[non_exhaustive]
888pub struct PodFins {
889 pub roll_rad: Vec<f64>,
893 pub axis_roll: [FinRollTerms; 2],
904 pub fins: u32,
906}
907
908#[derive(Debug, Clone, PartialEq, Serialize)]
918#[non_exhaustive]
919pub struct PodSetAero {
920 pub id: String,
922 pub copies: u32,
924 pub offset_squares_m2: f64,
926 pub fineness: f64,
929 crossflow_eta_low: f64,
931 pub bodies: Vec<BodyAero>,
935}
936
937impl FinSetAero {
938 pub fn cp_station_m(&self, mach: f64) -> Result<f64, AeroError> {
944 Ok(self.fore_station_m + self.fin.loading(mach)?.cp_m)
945 }
946}
947
948#[derive(Debug, Clone, PartialEq, Serialize)]
953pub struct AeroModel {
954 reference_area_m2: f64,
955 reference_diameter_m: f64,
956 length_m: f64,
957 max_body_radius_m: f64,
959 fineness: f64,
961 crossflow_eta_low: f64,
963 body_model: BodyModel,
965 bodies: Vec<BodyAero>,
966 pods: Vec<PodSetAero>,
968 #[serde(skip)]
970 supersonic_run: Option<SupersonicRun>,
971 #[serde(skip)]
973 supersonic_stop: Option<SupersonicFallback>,
974 #[serde(skip)]
978 supersonic: SupersonicTable,
979 fin_sets: Vec<FinSetAero>,
980 tube_fin_sets: Vec<TubeFinSetAero>,
981 drag_terms: Vec<ComponentDragTerms>,
982 drag_table: Option<DragTable>,
983 #[serde(skip_serializing_if = "Option::is_none")]
986 drag_model: Option<SharedDragModel>,
987 normal_force_table: Option<NormalForceTable>,
988 full_base_drag_under_power: bool,
991 #[serde(skip_serializing_if = "is_one")]
994 drag_scale: f64,
995}
996
997fn is_one(scale: &f64) -> bool {
1000 *scale == 1.0
1001}
1002
1003fn scaled(mut drag: Drag, scale: f64) -> Drag {
1006 if scale != 1.0 {
1007 for part in [
1008 &mut drag.zero_lift_coefficient,
1009 &mut drag.friction,
1010 &mut drag.pressure,
1011 &mut drag.base,
1012 &mut drag.parasitic,
1013 &mut drag.stated,
1014 ] {
1015 *part *= scale;
1016 }
1017 }
1018 drag
1019}
1020
1021fn apply_drag_overrides(
1034 layout: &Layout,
1035 drag_terms: &mut Vec<ComponentDragTerms>,
1036 body_terms: &[(usize, BodyGeometry)],
1037) -> Result<(), AeroError> {
1038 let check = |id: &str, coefficient: f64| {
1039 if coefficient.is_finite() && coefficient >= 0.0 {
1040 Ok(coefficient)
1041 } else {
1042 Err(AeroError::InComponent {
1043 id: id.to_owned(),
1044 source: Box::new(AeroError::Domain {
1045 what: "stated drag coefficient",
1046 value: coefficient,
1047 }),
1048 })
1049 }
1050 };
1051 let unsupported = |id: &str, what: &str| AeroError::InComponent {
1052 id: id.to_owned(),
1053 source: Box::new(AeroError::Unsupported(what.to_owned())),
1054 };
1055 for (k, stage) in layout.stages.iter().enumerate() {
1056 if let Some(stated) = stage.drag_override {
1057 check(&stage.id, stated.coefficient)?;
1058 if stage.hung_on.is_some() {
1061 return Err(unsupported(
1062 &stage.id,
1063 "a drag override on a parallel stage",
1064 ));
1065 }
1066 if stated.include_children && layout.stages.iter().any(|s| s.hung_on == Some(k)) {
1067 return Err(unsupported(
1068 &stage.id,
1069 "a drag override covering a stage a parallel stage hangs on",
1070 ));
1071 }
1072 }
1073 }
1074 let mut covered = Vec::with_capacity(layout.components.len());
1078 for (index, component) in layout.components.iter().enumerate() {
1079 if let Some(stated) = component.drag_override {
1080 check(&component.id, stated.coefficient)?;
1081 }
1082 let mut covers = layout
1083 .stages
1084 .get(component.stage)
1085 .and_then(|s| s.drag_override)
1086 .is_some_and(|o| o.include_children);
1087 let mut parent = component.parent;
1088 for _ in 0..layout.components.len() {
1089 let Some(up) = parent.and_then(|at| layout.components.get(at)) else {
1090 break;
1091 };
1092 covers |= up.drag_override.is_some_and(|o| o.include_children);
1093 parent = up.parent;
1094 }
1095 let in_pod =
1096 layout.pod_set_of(index).is_some() || matches!(component.part, Part::PodSet(_));
1097 if in_pod && component.drag_override.is_some() {
1098 return Err(unsupported(
1099 &component.id,
1100 "a drag override on a pod set or in a pod",
1101 ));
1102 }
1103 if in_pod && covers {
1104 return Err(unsupported(
1105 &component.id,
1106 "a drag override covering a pod set",
1107 ));
1108 }
1109 if component.drag_override.is_some() && matches!(component.part, Part::TubeFinSet(_)) {
1110 return Err(unsupported(
1111 &component.id,
1112 "a drag override on a tube fin set",
1113 ));
1114 }
1115 covered.push(covers);
1116 }
1117 for terms in drag_terms.iter_mut() {
1118 let Some((index, component)) = layout.find(&terms.id) else {
1119 continue;
1120 };
1121 terms.stated = if covered[index] {
1122 Some(0.0)
1123 } else if let Some(stated) = component.drag_override {
1124 let instances = match &component.part {
1125 Part::FinSet(set) => set.count,
1126 Part::LaunchLug(lug) => lug.count,
1127 Part::RailButton(button) => button.count,
1128 _ => 1,
1129 };
1130 Some(stated.coefficient * f64::from(instances))
1131 } else {
1132 None
1133 };
1134 }
1135 for pair in body_terms.windows(2) {
1137 let (ahead, behind) = (pair[0].0, pair[1].0);
1138 if drag_terms[ahead].stated.is_some() {
1139 let terms = &mut drag_terms[behind];
1140 terms.boattail_area_ratio -= terms.fore_step_down_area_ratio;
1141 terms.fore_step_down_area_ratio = 0.0;
1142 }
1143 }
1144 for stage in &layout.stages {
1145 if let Some(stated) = stage.drag_override {
1146 drag_terms.push(ComponentDragTerms::stage(&stage.id, stated.coefficient));
1147 }
1148 }
1149 Ok(())
1150}
1151
1152fn with_axial(mut drag: Drag, factor: f64) -> Result<Drag, AeroError> {
1159 drag.axial_coefficient = drag.zero_lift_coefficient * factor;
1160 if !(drag.zero_lift_coefficient.is_finite() && drag.axial_coefficient.is_finite()) {
1161 return Err(AeroError::Domain {
1162 what: "drag coefficient",
1163 value: drag.zero_lift_coefficient,
1164 });
1165 }
1166 Ok(drag)
1167}
1168
1169impl AeroModel {
1170 pub fn new(layout: &Layout) -> Result<Self, AeroError> {
1177 Self::with_body_model(layout, BodyModel::default())
1178 }
1179
1180 pub fn with_body_model(layout: &Layout, body_model: BodyModel) -> Result<Self, AeroError> {
1203 body_model.body_lift.validate()?;
1204 check_dimension("reference diameter", layout.reference_diameter_m, false)?;
1205 let reference_area_m2 = layout.reference_area_m2();
1206 let length_m = layout.length_m;
1207 check_dimension("rocket length", length_m, false)?;
1208 let mut max_radius: f64 = 0.0;
1209 for component in layout.body() {
1210 if let Some(radius) = component.part.max_radius_m()? {
1211 max_radius = max_radius.max(radius);
1212 }
1213 }
1214 check_dimension("maximum body radius", max_radius, false)?;
1215 let fineness = length_m / (2.0 * max_radius);
1216 let form_factor = body_friction_form_factor(fineness)?;
1217 let mut bodies = Vec::new();
1218 let mut fin_sets = Vec::new();
1219 let mut tube_fin_sets = Vec::new();
1220 let mut drag_terms = Vec::new();
1221 let mut previous_aft_area: Option<f64> = None;
1222 let mut body_terms_at = Vec::new();
1223 let mut last_body_terms: Option<usize> = None;
1224 let mut supersonic_segments = Vec::new();
1227 let mut supersonic_bounds = Vec::new();
1228 let mut supersonic_fore = Vec::new();
1229 let mut supersonic_boattails = Vec::new();
1230 let mut flare: Option<RunFlare> = None;
1231 let mut proposed_flare: Option<&hpr_design::Transition> = None;
1232 let mut supersonic_open = true;
1233 let mut behind_boattail = false;
1234 let (mut vertex_m, mut supersonic_end_m) = (0.0, 0.0);
1235 let mut pods: Vec<PodBuild> = Vec::new();
1237 let mut pod_at: Vec<Option<usize>> = vec![None; layout.components.len()];
1238 let motor_pod_sets = motor_pod_sets(layout)?;
1241 for (index, component) in layout.components.iter().enumerate() {
1242 let in_component = |e: AeroError| AeroError::InComponent {
1243 id: component.id.clone(),
1244 source: Box::new(e),
1245 };
1246 if !component.fore_station_m.is_finite() {
1247 return Err(in_component(AeroError::Domain {
1248 what: "component station",
1249 value: component.fore_station_m,
1250 }));
1251 }
1252 let pod = layout.pod_set_of(index).map(|set| pod_at[set]);
1254 let first_new_drag_term = drag_terms.len();
1255 if let Some(at) = pod
1256 && matches!(
1257 component.part,
1258 Part::NoseCone(_) | Part::Transition(_) | Part::BodyTube(_)
1259 )
1260 {
1261 let build = at.and_then(|at| pods.get_mut(at)).ok_or_else(|| {
1263 in_component(AeroError::Layout(
1264 "a pod's body component before its pod set".to_owned(),
1265 ))
1266 })?;
1267 build
1268 .add_body(component, &mut drag_terms, length_m, reference_area_m2)
1269 .map_err(in_component)?;
1270 continue;
1271 }
1272 let body = match &component.part {
1273 Part::NoseCone(nose) => Some(
1274 nose.profile()
1275 .map_err(AeroError::from)
1276 .and_then(|p| BodyGeometry::from_profile(&p)),
1277 ),
1278 Part::Transition(transition) => Some(
1279 transition
1280 .profile()
1281 .map_err(AeroError::from)
1282 .and_then(|p| BodyGeometry::from_profile(&p)),
1283 ),
1284 Part::BodyTube(tube) => {
1285 Some(BodyGeometry::cylinder(tube.length_m, tube.outer_radius_m))
1286 }
1287 Part::FinSet(set) => {
1288 let mut terms = (|| {
1289 let fin = FinAero::new(&set.planform, reference_area_m2)?;
1290 let body_radius = component.body_radius_m.ok_or_else(|| {
1291 AeroError::Layout(
1292 "a fin set needs the radius of the body tube it is on".to_owned(),
1293 )
1294 })?;
1295 if !set.cant_rad.is_finite() || set.cant_rad.abs() > MAX_CANT_RAD {
1298 return Err(AeroError::Domain {
1299 what: "fin cant",
1300 value: set.cant_rad,
1301 });
1302 }
1303 if set.count < 2 && set.cant_rad != 0.0 {
1304 return Err(AeroError::Unsupported(
1305 "cant on a single fin, whose side force isn't modeled".to_owned(),
1306 ));
1307 }
1308 let span = fin.geometry().span_m;
1309 let taper = fin.outline().tip_chord_m() / set.planform.root_chord_m();
1310 let roll = fin.roll_terms(body_radius, layout.reference_diameter_m)?;
1311 Ok(FinSetAero {
1312 id: component.id.clone(),
1313 count: set.count,
1314 base_angle_rad: set.base_angle_rad,
1315 count_factor: fin_count_factor(set.count)?,
1316 interference: interference_factor(span, body_radius)?,
1317 fore_station_m: component.fore_station_m,
1318 cant_rad: set.cant_rad,
1319 body_radius_m: body_radius,
1320 roll_forcing_interference: roll_forcing_interference(
1321 span,
1322 body_radius,
1323 )?,
1324 roll_damping_interference: roll_damping_interference(
1325 span,
1326 body_radius,
1327 taper,
1328 )?,
1329 roll,
1330 fin,
1331 pods: None,
1332 })
1333 })()
1334 .map_err(in_component)?;
1335 if pod.is_some() {
1336 terms.pods = Some(
1337 pod_fins(set, component, &terms, layout.reference_diameter_m)
1338 .map_err(in_component)?,
1339 );
1340 }
1341 drag_terms.push(
1342 ComponentDragTerms::fins(
1343 component,
1344 set,
1345 terms.fin.geometry(),
1346 length_m,
1347 reference_area_m2,
1348 )
1349 .map_err(in_component)?,
1350 );
1351 fin_sets.push(terms);
1352 None
1353 }
1354 Part::TubeFinSet(set) => {
1355 if pod.is_some() {
1356 return Err(in_component(AeroError::Unsupported(
1357 "tube fins on a pod".to_owned(),
1358 )));
1359 }
1360 let terms = TubeFinSetAero::new(component, set, reference_area_m2)
1361 .map_err(in_component)?;
1362 drag_terms.push(
1363 ComponentDragTerms::tube_fins(component, set, length_m, reference_area_m2)
1364 .map_err(in_component)?,
1365 );
1366 tube_fin_sets.push(terms);
1367 None
1368 }
1369 Part::PodSet(_) if !layout.components.iter().any(|c| c.parent == Some(index)) => {
1371 None
1372 }
1373 Part::PodSet(_) => {
1374 pod_at[index] = Some(pods.len());
1375 let mut build =
1376 PodBuild::new(layout, index, component).map_err(in_component)?;
1377 build.motor_pod_set = motor_pod_sets.iter().position(|&set| set == index);
1378 pods.push(build);
1379 None
1380 }
1381 Part::LaunchLug(lug) => {
1383 drag_terms.push(
1384 ComponentDragTerms::launch_lugs(
1385 component,
1386 lug,
1387 length_m,
1388 reference_area_m2,
1389 )
1390 .map_err(in_component)?,
1391 );
1392 None
1393 }
1394 Part::RailButton(button) => {
1395 drag_terms.push(
1396 ComponentDragTerms::rail_buttons(
1397 component,
1398 button,
1399 length_m,
1400 reference_area_m2,
1401 )
1402 .map_err(in_component)?,
1403 );
1404 None
1405 }
1406 Part::InnerTube(_)
1408 | Part::CenteringRing(_)
1409 | Part::MassComponent(_)
1410 | Part::Parachute(_)
1411 | Part::Streamer(_)
1412 | Part::ShockCord(_) => None,
1413 other => {
1414 return Err(in_component(AeroError::Unsupported(format!(
1415 "a {} part",
1416 other.kind_name()
1417 ))));
1418 }
1419 };
1420 if pod.is_some() {
1422 let copies = copy_count(component).map_err(in_component)?;
1423 for terms in &mut drag_terms[first_new_drag_term..] {
1424 terms.copies = copies;
1425 terms.in_pod = true;
1426 }
1427 }
1428 if let Some(geometry) = body {
1429 let geometry = geometry.map_err(in_component)?;
1430 let step = previous_aft_area.map_or(0.0, |aft| geometry.fore_area_m2 - aft);
1431 last_body_terms = Some(drag_terms.len());
1432 drag_terms.push(
1433 ComponentDragTerms::body(
1434 component,
1435 &geometry,
1436 match &component.part {
1437 Part::NoseCone(nose) => Some(nose.shape),
1438 Part::Transition(transition) => Some(transition.shape),
1439 _ => None,
1440 },
1441 previous_aft_area,
1442 form_factor,
1443 length_m,
1444 reference_area_m2,
1445 )
1446 .map_err(in_component)?,
1447 );
1448 body_terms_at.push((drag_terms.len() - 1, geometry));
1449 let segment = match &component.part {
1450 Part::NoseCone(nose) if bodies.is_empty() => nose
1451 .profile()
1452 .ok()
1453 .map(|profile| BodySegment::Profile { profile }),
1454 Part::BodyTube(tube)
1455 if !bodies.is_empty()
1456 && step.abs() <= 1e-6 * geometry.fore_area_m2
1457 && (component.fore_station_m - supersonic_end_m).abs()
1458 <= 1e-9 * length_m =>
1459 {
1460 Some(BodySegment::Cylinder {
1461 length_m: tube.length_m,
1462 radius_m: tube.outer_radius_m,
1463 })
1464 }
1465 Part::Transition(transition)
1467 if !bodies.is_empty()
1468 && transition.aft_radius_m < transition.fore_radius_m
1469 && step.abs() <= 1e-6 * geometry.fore_area_m2
1470 && (component.fore_station_m - supersonic_end_m).abs()
1471 <= 1e-9 * length_m =>
1472 {
1473 transition
1474 .profile()
1475 .ok()
1476 .map(|profile| BodySegment::Profile { profile })
1477 }
1478 Part::Transition(transition)
1486 if !bodies.is_empty()
1487 && !behind_boattail
1488 && transition.aft_radius_m > transition.fore_radius_m
1489 && matches!(transition.shape, NoseShape::Conical {})
1490 && body_model.supersonic_flare == SupersonicFlare::Marched
1491 && step.abs() <= 1e-6 * geometry.fore_area_m2
1492 && (component.fore_station_m - supersonic_end_m).abs()
1493 <= 1e-9 * length_m =>
1494 {
1495 proposed_flare = Some(transition);
1496 transition
1497 .profile()
1498 .ok()
1499 .map(|profile| BodySegment::Profile { profile })
1500 }
1501 _ => None,
1502 };
1503 match segment {
1504 Some(segment) if supersonic_open => {
1505 if supersonic_segments.is_empty() {
1506 vertex_m = component.fore_station_m;
1507 }
1508 if let Some(transition) = proposed_flare.take() {
1511 flare = Some(RunFlare {
1512 index: supersonic_segments.len(),
1513 ahead: supersonic_segments.clone(),
1514 fore_radius_m: transition.fore_radius_m,
1515 aft_radius_m: transition.aft_radius_m,
1516 length_m: transition.length_m,
1517 clipped: transition.clipped,
1518 });
1519 }
1520 supersonic_boattails.push(match (&component.part, body_model) {
1523 (
1524 Part::Transition(transition),
1525 BodyModel {
1526 supersonic_boattail: SupersonicBoattail::WashingtonPettis,
1527 ..
1528 },
1529 ) if transition.aft_radius_m < transition.fore_radius_m => {
1530 let mut in_its_place = supersonic_segments.clone();
1531 in_its_place.push(BodySegment::Cylinder {
1532 length_m: transition.length_m,
1533 radius_m: transition.fore_radius_m,
1534 });
1535 Some(RunBoattail {
1536 in_its_place,
1537 fore_radius_m: transition.fore_radius_m,
1538 aft_radius_m: transition.aft_radius_m,
1539 length_m: transition.length_m,
1540 })
1541 }
1542 _ => None,
1543 });
1544 supersonic_segments.push(segment);
1545 let fore_m = component.fore_station_m;
1546 supersonic_fore.push(fore_m);
1547 supersonic_end_m = fore_m + geometry.length_m;
1548 behind_boattail |= matches!(&component.part, Part::Transition(t)
1553 if t.aft_radius_m < t.fore_radius_m);
1554 supersonic_bounds
1555 .push((!behind_boattail).then_some((fore_m, supersonic_end_m)));
1556 supersonic_open &= flare.is_none();
1558 }
1559 _ => supersonic_open = false,
1560 }
1561 previous_aft_area = Some(geometry.aft_area_m2);
1562 bodies.push(body_terms(component, geometry, step, reference_area_m2));
1563 }
1564 }
1565 if let (Some(index), Some(last)) = (last_body_terms, bodies.last()) {
1567 drag_terms[index].base_area_m2 = last.geometry.aft_area_m2;
1568 }
1569 couple_afterbody(&mut drag_terms, &body_terms_at, reference_area_m2)?;
1571 let mut pod_sets = Vec::new();
1573 for build in pods {
1574 if let Some((index, geometry)) = build.terms_at.last() {
1575 drag_terms[*index].base_area_m2 = geometry.aft_area_m2;
1576 }
1577 couple_afterbody(&mut drag_terms, &build.terms_at, reference_area_m2)?;
1578 if !build.aero.bodies.is_empty() {
1579 pod_sets.push(build.aero);
1580 }
1581 }
1582 apply_drag_overrides(layout, &mut drag_terms, &body_terms_at)?;
1584 let marched = supersonic_segments.len();
1588 let mut shelter_weight = 1.0_f64;
1592 let mut lip_too_long = false;
1595 let sheltered_lips = (marched..bodies.len())
1596 .take_while(|&index| {
1597 let (term, geometry) = body_terms_at[index];
1598 if bodies[index].slope_per_rad <= 0.0 {
1602 return false;
1603 }
1604 let Some(wake) = drag_terms[term].in_wake_of else {
1607 return false;
1608 };
1609 let boattail = &wake.boattail;
1610 if geometry.length_m > boattail.fore_diameter_m - boattail.aft_diameter_m {
1611 lip_too_long = true;
1612 return false;
1613 }
1614 let steps = body_terms_at[index.saturating_sub(1)]
1618 .1
1619 .aft_area_m2
1620 .lt(&geometry.fore_area_m2);
1621 let shoulders = geometry.aft_area_m2 > geometry.fore_area_m2;
1622 let mut covered = 1.0_f64;
1630 if steps {
1631 covered = covered.min(wake.step_fraction);
1632 }
1633 if shoulders {
1634 covered = covered.min(wake.shoulder_fraction);
1635 }
1636 if !covered.is_finite() || covered <= SUPERSONIC_SHELTER_FLOOR {
1639 return false;
1640 }
1641 shelter_weight = shelter_weight.min(covered);
1642 true
1643 })
1644 .count();
1645 let covered = marched + sheltered_lips;
1646 let rest_carries_nothing = bodies[covered..]
1647 .iter()
1648 .all(|body| body.slope_per_rad.abs() <= 1e-9);
1649 let carrying = (covered..bodies.len()).find(|&index| {
1653 let nothing = bodies[index].slope_per_rad.abs() <= 1e-9;
1654 !nothing
1655 });
1656 let supersonic_stop = if marched > 0 && rest_carries_nothing {
1657 None
1658 } else {
1659 Some(match carrying.map(|index| (index, &bodies[index])) {
1660 Some((index, body)) if lip_too_long && index == covered => {
1661 SupersonicFallback::LongLip {
1662 component: body.id.clone(),
1663 }
1664 }
1665 Some((_, body))
1667 if marched > 0
1668 && (body.step_area_m2.abs() > 1e-6 * body.geometry.fore_area_m2
1669 || body.step_area_m2.is_nan()) =>
1670 {
1671 SupersonicFallback::RadiusStep {
1672 component: body.id.clone(),
1673 }
1674 }
1675 _ => SupersonicFallback::Other,
1676 })
1677 };
1678 let supersonic_run = (marched > 0 && rest_carries_nothing).then_some(SupersonicRun {
1679 segments: supersonic_segments,
1680 vertex_m,
1681 bounds_m: supersonic_bounds,
1682 fore_m: supersonic_fore,
1683 boattails: supersonic_boattails,
1684 flare,
1685 sheltered_lips,
1686 shape_weight: shelter_weight,
1687 });
1688 Ok(Self {
1689 reference_area_m2,
1690 reference_diameter_m: layout.reference_diameter_m,
1691 length_m,
1692 max_body_radius_m: max_radius,
1693 fineness,
1694 crossflow_eta_low: crate::crossflow::crossflow_eta_low(fineness),
1695 body_model,
1696 bodies,
1697 pods: pod_sets,
1698 supersonic_run,
1699 supersonic_stop,
1700 supersonic: SupersonicTable::default(),
1701 fin_sets,
1702 tube_fin_sets,
1703 drag_terms,
1704 drag_table: None,
1705 drag_model: None,
1706 normal_force_table: None,
1707 full_base_drag_under_power: false,
1708 drag_scale: 1.0,
1709 })
1710 }
1711
1712 #[must_use]
1715 pub fn with_drag_table(mut self, table: DragTable) -> Self {
1716 self.drag_table = Some(table);
1717 self.drag_model = None;
1718 self
1719 }
1720
1721 pub fn drag_table(&self) -> Option<&DragTable> {
1723 self.drag_table.as_ref()
1724 }
1725
1726 #[must_use]
1729 pub fn with_drag_model(self, model: impl DragModel + 'static) -> Self {
1730 self.with_shared_drag_model(Arc::new(model))
1731 }
1732
1733 #[must_use]
1736 pub fn with_shared_drag_model(mut self, model: Arc<dyn DragModel>) -> Self {
1737 self.drag_model = Some(SharedDragModel(model));
1738 self.drag_table = None;
1739 self
1740 }
1741
1742 pub fn drag_model(&self) -> Option<&Arc<dyn DragModel>> {
1744 self.drag_model.as_ref().map(|shared| &shared.0)
1745 }
1746
1747 #[must_use]
1761 pub fn with_full_base_drag_under_power(mut self) -> Self {
1762 self.full_base_drag_under_power = true;
1763 self
1764 }
1765
1766 pub fn full_base_drag_under_power(&self) -> bool {
1769 self.full_base_drag_under_power
1770 }
1771
1772 pub fn with_drag_scale(mut self, scale: f64) -> Result<Self, AeroError> {
1784 if !(scale.is_finite() && scale >= 0.0) {
1785 return Err(AeroError::Domain {
1786 what: "drag scale",
1787 value: scale,
1788 });
1789 }
1790 self.drag_scale = scale;
1791 Ok(self)
1792 }
1793
1794 pub fn drag_scale(&self) -> f64 {
1797 self.drag_scale
1798 }
1799
1800 pub fn supersonic_table_built(&self) -> bool {
1804 self.supersonic.0.get().is_some()
1805 }
1806
1807 pub fn share_supersonic_table(&mut self, other: &AeroModel) -> bool {
1821 let shared = self.supersonic_run.is_some()
1822 && self.supersonic_run == other.supersonic_run
1823 && self.reference_area_m2 == other.reference_area_m2;
1824 if shared {
1825 self.supersonic = other.supersonic.clone();
1826 }
1827 shared
1828 }
1829
1830 fn read(&self, conditions: &DragConditions) -> DragConditions {
1833 if self.full_base_drag_under_power {
1834 DragConditions {
1835 thrusting_motor_area_m2: 0.0,
1836 thrusting_pod_motor_areas_m2: [0.0; MOTOR_POD_SETS],
1837 ..*conditions
1838 }
1839 } else {
1840 *conditions
1841 }
1842 }
1843
1844 pub fn with_normal_force_table(mut self, table: NormalForceTable) -> Result<Self, AeroError> {
1856 for column in table.columns() {
1857 let knots = column.cp_station_m.xs();
1858 let used = knots.partition_point(|&mach| mach <= NORMAL_FORCE_MACH_LIMIT) + 1;
1859 for &cp in column.cp_station_m.ys().iter().take(used) {
1860 if !(0.0..=self.length_m).contains(&cp) {
1861 return Err(AeroError::Domain {
1862 what: "normal-force table center of pressure, m aft of the nose tip",
1863 value: cp,
1864 });
1865 }
1866 }
1867 }
1868 self.normal_force_table = Some(table);
1869 Ok(self)
1870 }
1871
1872 pub fn normal_force_table(&self) -> Option<&NormalForceTable> {
1874 self.normal_force_table.as_ref()
1875 }
1876
1877 pub fn drag_terms(&self) -> &[ComponentDragTerms] {
1879 &self.drag_terms
1880 }
1881
1882 pub fn length_m(&self) -> f64 {
1885 self.length_m
1886 }
1887
1888 pub fn drag(&self, flow: &Flow, conditions: &DragConditions) -> Result<Drag, AeroError> {
1908 conditions.validate()?;
1909 let conditions = &self.read(conditions);
1910 let factor = axial_drag_alpha_factor(flow.alpha_rad)?;
1911 let drag = if let Some(custom) = &self.drag_model {
1912 flow.validate_angles()?;
1913 if !(flow.mach.is_finite() && flow.mach >= 0.0) {
1914 return Err(AeroError::Domain {
1915 what: "Mach number for a drag model",
1916 value: flow.mach,
1917 });
1918 }
1919 let coefficient = custom
1920 .0
1921 .zero_lift_drag(&DragQuery::new(flow, conditions, self))
1922 .map_err(|source| AeroError::DragModel {
1923 source: Box::new(source),
1924 })?;
1925 if !(coefficient.is_finite() && coefficient >= 0.0) {
1926 return Err(AeroError::Domain {
1927 what: "zero-lift drag coefficient from a drag model",
1928 value: coefficient,
1929 });
1930 }
1931 Drag {
1932 zero_lift_coefficient: coefficient,
1933 ..Drag::default()
1934 }
1935 } else if let Some(table) = &self.drag_table {
1936 flow.validate_angles()?;
1937 let lookup = table.lookup(flow.mach, conditions.thrusting)?;
1938 if lookup.value < 0.0 {
1942 return Err(AeroError::Domain {
1943 what: "zero-lift drag coefficient from a drag table",
1944 value: lookup.value,
1945 });
1946 }
1947 let scale = match table.reference_diameter_m {
1948 Some(d) => {
1949 check_dimension("drag table reference diameter", d, false)?;
1950 0.25 * PI * d * d / self.reference_area_m2
1951 }
1952 None => 1.0,
1953 };
1954 Drag {
1955 zero_lift_coefficient: lookup.value * scale,
1956 table: Some(lookup),
1957 ..Drag::default()
1958 }
1959 } else {
1960 self.buildup_sum(flow, conditions)?
1961 };
1962 with_axial(scaled(drag, self.drag_scale), factor)
1963 }
1964
1965 pub fn buildup_drag(
1974 &self,
1975 flow: &Flow,
1976 conditions: &DragConditions,
1977 ) -> Result<Drag, AeroError> {
1978 conditions.validate()?;
1979 let conditions = &self.read(conditions);
1980 let factor = axial_drag_alpha_factor(flow.alpha_rad)?;
1981 with_axial(self.buildup_sum(flow, conditions)?, factor)
1982 }
1983
1984 fn buildup_sum(&self, flow: &Flow, conditions: &DragConditions) -> Result<Drag, AeroError> {
1987 flow.validate_for_buildup()?;
1988 self.check_tube_fins(flow.mach)?;
1990 let reynolds = conditions.reynolds_per_m * self.length_m;
1991 let mut sum = Drag::default();
1992 for terms in &self.drag_terms {
1993 let d = terms.evaluate(reynolds, flow.mach, conditions, self.reference_area_m2)?;
1994 sum.friction += d.friction;
1995 sum.pressure += d.pressure;
1996 sum.base += d.base;
1997 sum.parasitic += d.parasitic;
1998 sum.stated += d.stated;
1999 }
2000 sum.zero_lift_coefficient =
2001 sum.friction + sum.pressure + sum.base + sum.parasitic + sum.stated;
2002 Ok(sum)
2003 }
2004
2005 pub fn buildup_components(
2015 &self,
2016 flow: &Flow,
2017 conditions: &DragConditions,
2018 ) -> Result<Vec<ComponentDrag>, AeroError> {
2019 flow.validate_for_buildup()?;
2020 self.check_tube_fins(flow.mach)?;
2021 conditions.validate()?;
2022 let conditions = &self.read(conditions);
2023 let factor = axial_drag_alpha_factor(flow.alpha_rad)?;
2024 let reynolds = conditions.reynolds_per_m * self.length_m;
2025 self.drag_terms
2026 .iter()
2027 .map(|terms| {
2028 let mut drag =
2029 terms.evaluate(reynolds, flow.mach, conditions, self.reference_area_m2)?;
2030 drag.axial_coefficient *= factor;
2031 Ok(ComponentDrag {
2032 id: terms.id.clone(),
2033 drag,
2034 })
2035 })
2036 .collect()
2037 }
2038
2039 pub fn reference_area_m2(&self) -> f64 {
2041 self.reference_area_m2
2042 }
2043
2044 pub fn reference_diameter_m(&self) -> f64 {
2046 self.reference_diameter_m
2047 }
2048
2049 pub fn roll(&self, mach: f64) -> Result<Roll, AeroError> {
2064 check_mach(mach, NORMAL_FORCE_MACH_LIMIT, "the roll moment")?;
2065 let mut roll = Roll::default();
2066 for set in &self.fin_sets {
2067 if let Some(pods) = &set.pods {
2068 let [a, b] = &pods.axis_roll;
2070 let damping =
2071 0.5 * (set.fin.roll_with(a, mach).damping + set.fin.roll_with(b, mach).damping);
2072 roll.damping += f64::from(pods.fins) * damping * set.roll_damping_interference;
2073 continue;
2074 }
2075 let fin = set.fin.roll_with(&set.roll, mach);
2076 let n = f64::from(set.count);
2077 roll.forcing -= n * fin.forcing_per_rad * set.roll_forcing_interference * set.cant_rad;
2078 roll.damping += n * fin.damping * set.roll_damping_interference;
2079 }
2080 let d = self.reference_diameter_m;
2084 for pod in &self.pods {
2085 for body in &pod.bodies {
2086 roll.damping -= 2.0 * body.slope_per_rad * pod.offset_squares_m2 / (d * d);
2087 }
2088 }
2089 for set in &self.tube_fin_sets {
2091 roll.damping += set.roll_damping(mach, d)?;
2092 }
2093 Ok(roll)
2094 }
2095
2096 pub fn steady_roll_rate_rad_s(&self, mach: f64, speed_m_s: f64) -> Result<f64, AeroError> {
2103 check_dimension("airspeed", speed_m_s, true)?;
2104 let roll = self.roll(mach)?;
2105 Ok(if roll.damping < 0.0 {
2106 -roll.forcing / roll.damping * 2.0 * speed_m_s / self.reference_diameter_m
2107 } else {
2108 0.0
2109 })
2110 }
2111
2112 pub fn bodies(&self) -> &[BodyAero] {
2114 &self.bodies
2115 }
2116
2117 pub fn supersonic_body(&self) -> Option<&SupersonicBody> {
2126 let run = self.supersonic_run.as_ref()?;
2127 self.supersonic
2128 .0
2129 .get_or_init(|| SupersonicBody::new(run, self.reference_area_m2))
2130 .as_ref()
2131 }
2132
2133 pub fn supersonic_fallback(&self) -> Option<SupersonicFallback> {
2137 if self.supersonic_run.is_none() {
2138 return Some(
2139 self.supersonic_stop
2140 .clone()
2141 .unwrap_or(SupersonicFallback::Other),
2142 );
2143 }
2144 if self.supersonic_body().is_some() {
2145 return None;
2146 }
2147 let steep_tip = self.supersonic_run.as_ref().is_some_and(|run| {
2149 ShockExpansionBody::new(&run.segments, DEFAULT_ELEMENTS_PER_CURVE).is_ok_and(|body| {
2150 !body.has_blunt_tip()
2151 && crate::shock_expansion::past_cone_tables(body.vertex_angle_rad())
2152 })
2153 });
2154 if steep_tip {
2155 return Some(SupersonicFallback::SteepTip);
2156 }
2157 let marched = self
2161 .supersonic_run
2162 .as_ref()
2163 .map_or(0, |run| run.segments.len());
2164 let stepped = self
2165 .bodies
2166 .iter()
2167 .take(marched)
2168 .find(|body| body.step_area_m2 != 0.0);
2169 Some(match stepped {
2170 Some(body) => SupersonicFallback::RadiusStep {
2171 component: body.id.clone(),
2172 },
2173 None => SupersonicFallback::Other,
2174 })
2175 }
2176
2177 fn supersonic_at(&self, mach: f64) -> Option<&SupersonicBody> {
2180 (mach > SUPERSONIC_JOIN_START_MACH)
2181 .then(|| self.supersonic_body())
2182 .flatten()
2183 }
2184
2185 fn body_potential(&self, index: usize, body: &BodyAero, mach: f64) -> (f64, f64) {
2189 let (slope, moment) = (body.slope_per_rad, body.moment_slope_m);
2190 match self.supersonic_at(mach) {
2191 Some(s) if index < s.covered && mach > s.join_start_mach => {
2192 let w = s.weight(mach);
2193 let Some((se_slope, se_moment)) = s.share(index, mach) else {
2194 return (slope, moment);
2195 };
2196 (
2197 slope + w * (se_slope - slope),
2198 moment + w * (se_moment - moment),
2199 )
2200 }
2201 _ => (slope, moment),
2202 }
2203 }
2204
2205 pub fn body_lift_factor(&self, flow: &Flow) -> f64 {
2209 self.lift_factor_at(flow.mach, flow.alpha_rad.sin())
2210 }
2211
2212 fn lift_factor_at(&self, mach: f64, sin_alpha: f64) -> f64 {
2214 self.lift_factor_of(self.fineness, self.crossflow_eta_low, mach, sin_alpha)
2215 }
2216
2217 fn lift_factor_of(&self, fineness: f64, eta_low: f64, mach: f64, sin_alpha: f64) -> f64 {
2220 let crossflow_mach = mach * sin_alpha.abs();
2221 match self.body_model.body_lift {
2222 BodyLift::Jorgensen {} => {
2223 crate::crossflow::crossflow_factor_from_eta_low(eta_low, crossflow_mach)
2224 }
2225 other => other.factor(fineness, crossflow_mach),
2226 }
2227 }
2228
2229 fn body_factors(&self, flow: &Flow) -> (f64, f64) {
2232 let (potential, lift, sin_alpha) = alpha_factors(flow.alpha_rad);
2233 if lift == 0.0 {
2234 (potential, 0.0)
2235 } else {
2236 (potential, lift * self.lift_factor_at(flow.mach, sin_alpha))
2237 }
2238 }
2239
2240 fn pod_factors(&self, pod: &PodSetAero, flow: &Flow) -> (f64, f64) {
2242 let (potential, lift, sin_alpha) = alpha_factors(flow.alpha_rad);
2243 if lift == 0.0 {
2244 (potential, 0.0)
2245 } else {
2246 let factor =
2247 self.lift_factor_of(pod.fineness, pod.crossflow_eta_low, flow.mach, sin_alpha);
2248 (potential, lift * factor)
2249 }
2250 }
2251
2252 fn pod_body(&self, index: usize) -> Option<(&PodSetAero, &BodyAero)> {
2255 let mut index = index;
2256 for pod in &self.pods {
2257 match pod.bodies.get(index) {
2258 Some(body) => return Some((pod, body)),
2259 None => index -= pod.bodies.len(),
2260 }
2261 }
2262 None
2263 }
2264
2265 fn pod_bodies(&self) -> impl Iterator<Item = (&PodSetAero, &BodyAero)> {
2267 self.pods
2268 .iter()
2269 .flat_map(|pod| pod.bodies.iter().map(move |body| (pod, body)))
2270 }
2271
2272 pub fn body_model(&self) -> BodyModel {
2274 self.body_model
2275 }
2276
2277 pub fn fineness(&self) -> f64 {
2279 self.fineness
2280 }
2281
2282 pub fn fin_sets(&self) -> &[FinSetAero] {
2284 &self.fin_sets
2285 }
2286
2287 pub fn rolls(&self) -> bool {
2291 self.normal_force_table.is_none() && self.fin_sets.iter().any(|set| set.count < 3)
2292 }
2293
2294 pub fn tube_fin_sets(&self) -> &[TubeFinSetAero] {
2296 &self.tube_fin_sets
2297 }
2298
2299 fn check_tube_fins(&self, mach: f64) -> Result<(), AeroError> {
2301 if self.tube_fin_sets.is_empty() {
2302 Ok(())
2303 } else {
2304 check_tube_fin_mach(mach)
2305 }
2306 }
2307
2308 pub fn pod_sets(&self) -> &[PodSetAero] {
2310 &self.pods
2311 }
2312
2313 fn terms<'a>(&'a self, flow: &Flow) -> impl Iterator<Item = (&'a str, Term)> + 'a {
2317 let (potential, lift) = self.body_factors(flow);
2318 let (mach, roll) = (flow.mach, flow.roll_rad);
2319 let bodies = self.bodies.iter().enumerate().map(move |(index, body)| {
2320 let (slope, moment) = self.body_potential(index, body, mach);
2321 (
2322 body.id.as_str(),
2323 body_term(body, slope, moment, potential, lift),
2324 )
2325 });
2326 let flow = *flow;
2327 let pods = self
2328 .pod_bodies()
2329 .map(move |(pod, body)| (body.id.as_str(), self.pod_body_term(pod, body, &flow)));
2330 let fins = self
2331 .fin_sets
2332 .iter()
2333 .map(move |set| (set.id.as_str(), fin_term(set, mach, roll)));
2334 let tubes = self
2335 .tube_fin_sets
2336 .iter()
2337 .map(move |set| (set.id.as_str(), tube_fin_term(set, mach)));
2338 bodies.chain(pods).chain(fins).chain(tubes)
2339 }
2340
2341 fn pod_body_term(&self, pod: &PodSetAero, body: &BodyAero, flow: &Flow) -> Term {
2344 let (potential, lift) = self.pod_factors(pod, flow);
2345 body_term(
2346 body,
2347 body.slope_per_rad,
2348 body.moment_slope_m,
2349 potential,
2350 lift,
2351 )
2352 .times(f64::from(pod.copies))
2353 }
2354
2355 pub fn component_count(&self) -> usize {
2358 self.tube_fin_set_start() + self.tube_fin_sets.len()
2359 }
2360
2361 pub fn tube_fin_set_start(&self) -> usize {
2364 self.fin_set_start() + self.fin_sets.len()
2365 }
2366
2367 pub fn fin_set_start(&self) -> usize {
2371 self.bodies.len() + self.pods.iter().map(|pod| pod.bodies.len()).sum::<usize>()
2372 }
2373
2374 pub fn component_normal_force(
2383 &self,
2384 index: usize,
2385 flow: &Flow,
2386 ) -> Result<NormalForce, AeroError> {
2387 flow.validate()?;
2388 let term = if let Some(body) = self.bodies.get(index) {
2389 let (potential, lift) = self.body_factors(flow);
2390 let (slope, moment) = self.body_potential(index, body, flow.mach);
2391 body_term(body, slope, moment, potential, lift)
2392 } else if let Some((pod, body)) = self.pod_body(index - self.bodies.len()) {
2393 self.pod_body_term(pod, body, flow)
2394 } else if let Some(set) = self.fin_sets.get(index - self.fin_set_start()) {
2395 fin_term(set, flow.mach, flow.roll_rad)
2396 } else if let Some(set) = self.tube_fin_sets.get(index - self.tube_fin_set_start()) {
2397 check_tube_fin_mach(flow.mach)?;
2398 tube_fin_term(set, flow.mach)
2399 } else {
2400 return Err(AeroError::Domain {
2401 what: "component index",
2402 value: index as f64,
2403 });
2404 };
2405 Ok(NormalForce::new(term, flow.alpha_rad))
2406 }
2407
2408 pub fn component_station_m(&self, index: usize, mach: f64) -> Result<f64, AeroError> {
2430 check_mach(mach, NORMAL_FORCE_MACH_LIMIT, "the normal force")?;
2431 if let Some(body) = self.bodies.get(index) {
2432 let station = slender_station_m(body, self.reference_area_m2);
2433 Ok(match self.supersonic_at(mach) {
2437 Some(s) if index < s.covered && s.stationed[index] && mach > s.join_start_mach => {
2438 match s.share(index, mach) {
2439 Some((slope, moment)) => {
2442 station + s.weight(mach) * (moment / slope - station)
2443 }
2444 None => station,
2445 }
2446 }
2447 _ => station,
2448 })
2449 } else if let Some((_, body)) = self.pod_body(index - self.bodies.len()) {
2450 Ok(slender_station_m(body, self.reference_area_m2))
2451 } else if let Some(set) = self.fin_sets.get(index - self.fin_set_start()) {
2452 Ok(set.fore_station_m + set.fin.loading_at(mach).cp_m)
2453 } else {
2454 let set = self
2455 .tube_fin_sets
2456 .get(index - self.tube_fin_set_start())
2457 .ok_or(AeroError::Domain {
2458 what: "component index",
2459 value: index as f64,
2460 })?;
2461 check_tube_fin_mach(mach)?;
2462 Ok(set.loading_at(mach).1)
2463 }
2464 }
2465
2466 pub fn normal_force(&self, flow: &Flow) -> Result<NormalForce, AeroError> {
2480 if let Some(table) = &self.normal_force_table {
2481 flow.validate_angles()?;
2482 let lookup = table.lookup_within(flow.mach, flow.alpha_rad, (0.0, self.length_m))?;
2483 let area_m2 = match table.reference() {
2484 TableReference::Diameter { diameter_m } => 0.25 * PI * diameter_m * diameter_m,
2485 TableReference::LargestBody => PI * self.max_body_radius_m * self.max_body_radius_m,
2486 TableReference::Rocket => self.reference_area_m2,
2488 };
2489 let scale = area_m2 / self.reference_area_m2;
2490 let coefficient = lookup.coefficient * scale;
2491 let slope = lookup.slope_per_rad * scale;
2492 return Ok(NormalForce {
2493 coefficient,
2494 slope_per_rad: slope,
2495 moment_m: coefficient * lookup.cp_station_m,
2496 moment_slope_m: slope * lookup.cp_station_m,
2497 cp_station_m: (slope != 0.0).then_some(lookup.cp_station_m),
2498 side_coefficient: 0.0,
2499 side_moment_m: 0.0,
2500 table: Some(lookup),
2501 });
2502 }
2503 flow.validate()?;
2504 self.check_tube_fins(flow.mach)?;
2505 let total = self
2506 .terms(flow)
2507 .fold(Term::default(), |sum, (_, term)| sum.add(term));
2508 Ok(NormalForce::new(total, flow.alpha_rad))
2509 }
2510
2511 pub fn components(&self, flow: &Flow) -> Result<Vec<ComponentNormalForce>, AeroError> {
2523 flow.validate()?;
2524 self.check_tube_fins(flow.mach)?;
2525 Ok(self
2526 .terms(flow)
2527 .map(|(id, term)| ComponentNormalForce {
2528 id: id.to_owned(),
2529 normal_force: NormalForce::new(term, flow.alpha_rad),
2530 })
2531 .collect())
2532 }
2533}
2534
2535fn slender_station_m(body: &BodyAero, reference_area_m2: f64) -> f64 {
2540 let step_slope = 2.0 * body.step_area_m2 / reference_area_m2;
2541 let scale = (body.slope_per_rad - step_slope).abs() + step_slope.abs();
2542 if body.slope_per_rad.abs() <= 1e-12 * scale {
2543 body.lift_station_m
2544 } else {
2545 body.moment_slope_m / body.slope_per_rad
2546 }
2547}
2548
2549fn fin_term(set: &FinSetAero, mach: f64, roll: f64) -> Term {
2552 let FinLoading {
2553 slope_per_rad,
2554 cp_m,
2555 } = set.fin.loading_at(mach);
2556 let per_set = slope_per_rad * set.count_factor * set.interference;
2557 let station = set.fore_station_m + cp_m;
2558 let (roll_share, side_share) = match &set.pods {
2559 None => (
2560 roll_sum(set.count, set.base_angle_rad, roll),
2561 side_sum(set.count, set.base_angle_rad, roll),
2562 ),
2563 Some(pods) => pods.roll_rad.iter().fold((0.0, 0.0), |(r, s), turn| {
2564 let base = set.base_angle_rad + turn;
2565 (
2566 r + roll_sum(set.count, base, roll),
2567 s + side_sum(set.count, base, roll),
2568 )
2569 }),
2570 };
2571 let slope = per_set * roll_share;
2572 let side = per_set * side_share;
2573 Term {
2574 slope,
2575 moment: slope * station,
2576 side,
2577 side_moment: side * station,
2578 scale: slope.abs(),
2579 }
2580}
2581
2582fn tube_fin_term(set: &TubeFinSetAero, mach: f64) -> Term {
2585 let (slope, station) = set.loading_at(mach);
2586 Term {
2587 slope,
2588 moment: slope * station,
2589 side: 0.0,
2590 side_moment: 0.0,
2591 scale: slope.abs(),
2592 }
2593}
2594
2595fn body_term(body: &BodyAero, slope: f64, moment: f64, potential: f64, lift: f64) -> Term {
2599 let (attached, lift) = (slope * potential, body.planform_ratio * lift);
2600 Term {
2601 slope: attached + lift,
2602 moment: moment * potential + lift * body.lift_station_m,
2603 scale: attached.abs() + lift.abs(),
2604 ..Term::default()
2605 }
2606}
2607
2608fn alpha_factors(alpha_rad: f64) -> (f64, f64, f64) {
2611 let (s, sin) = (sinc(alpha_rad), alpha_rad.sin());
2612 (s, sin * s, sin)
2613}
2614
2615fn motor_pod_sets(layout: &Layout) -> Result<Vec<usize>, AeroError> {
2623 let sets = layout.motor_pod_sets();
2624 if sets.len() > MOTOR_POD_SETS {
2625 return Err(AeroError::Unsupported(format!(
2626 "motor mounts in {} pod sets; the drag tells apart the thrusting motors' areas of at \
2627 most {MOTOR_POD_SETS}",
2628 sets.len()
2629 )));
2630 }
2631 Ok(sets)
2632}
2633
2634fn copy_count(component: &PlacedComponent) -> Result<u32, AeroError> {
2636 u32::try_from(component.copies.len()).map_err(|_| AeroError::Domain {
2637 what: "copies of a component",
2638 value: component.copies.len() as f64,
2639 })
2640}
2641
2642struct PodBuild {
2644 aero: PodSetAero,
2645 form_factor: Option<f64>,
2647 motor_pod_set: Option<usize>,
2650 previous_aft_area: Option<f64>,
2651 terms_at: Vec<(usize, BodyGeometry)>,
2653}
2654
2655impl PodBuild {
2656 fn new(layout: &Layout, index: usize, component: &PlacedComponent) -> Result<Self, AeroError> {
2658 let copies = component.contents_copies()?;
2659 let offset_squares_m2 = copies
2660 .iter()
2661 .map(|c| c.offset_m[0] * c.offset_m[0] + c.offset_m[1] * c.offset_m[1])
2662 .sum();
2663 let mut max_radius: f64 = 0.0;
2664 for child in layout.components.iter().filter(|c| c.parent == Some(index)) {
2665 if let Some(radius) = child.part.max_radius_m()? {
2666 max_radius = max_radius.max(radius);
2667 }
2668 }
2669 let fineness = if max_radius > 0.0 && component.length_m > 0.0 {
2671 component.length_m / (2.0 * max_radius)
2672 } else {
2673 0.0
2674 };
2675 let form_factor = if fineness > 0.0 {
2676 Some(body_friction_form_factor(fineness)?)
2677 } else {
2678 None
2679 };
2680 Ok(Self {
2681 aero: PodSetAero {
2682 id: component.id.clone(),
2683 copies: u32::try_from(copies.len()).map_err(|_| AeroError::Domain {
2684 what: "pods in a pod set",
2685 value: copies.len() as f64,
2686 })?,
2687 offset_squares_m2,
2688 fineness,
2689 crossflow_eta_low: if fineness > 0.0 {
2690 crate::crossflow::crossflow_eta_low(fineness)
2691 } else {
2692 0.0
2693 },
2694 bodies: Vec::new(),
2695 },
2696 form_factor,
2697 motor_pod_set: None,
2698 previous_aft_area: None,
2699 terms_at: Vec::new(),
2700 })
2701 }
2702
2703 fn add_body(
2707 &mut self,
2708 component: &PlacedComponent,
2709 drag_terms: &mut Vec<ComponentDragTerms>,
2710 length_m: f64,
2711 reference_area_m2: f64,
2712 ) -> Result<(), AeroError> {
2713 let (geometry, shape) = match &component.part {
2714 Part::NoseCone(nose) => (
2715 BodyGeometry::from_profile(&nose.profile()?)?,
2716 Some(nose.shape),
2717 ),
2718 Part::Transition(transition) => (
2719 BodyGeometry::from_profile(&transition.profile()?)?,
2720 Some(transition.shape),
2721 ),
2722 Part::BodyTube(tube) if tube.length_m == 0.0 => {
2723 if tube.outer_radius_m == 0.0 {
2724 return Ok(());
2725 }
2726 return Err(AeroError::Unsupported(
2727 "a pod's tube of no length with a radius: a flat disc, which the drag buildup \
2728 has no term for"
2729 .to_owned(),
2730 ));
2731 }
2732 Part::BodyTube(tube) => (
2733 BodyGeometry::cylinder(tube.length_m, tube.outer_radius_m)?,
2734 None,
2735 ),
2736 _ => return Ok(()),
2737 };
2738 let form_factor = self.form_factor.ok_or(AeroError::Domain {
2740 what: "pod fineness",
2741 value: self.aero.fineness,
2742 })?;
2743 let step = self
2744 .previous_aft_area
2745 .map_or(0.0, |aft| geometry.fore_area_m2 - aft);
2746 let mut terms = ComponentDragTerms::body(
2747 component,
2748 &geometry,
2749 shape,
2750 self.previous_aft_area,
2751 form_factor,
2752 length_m,
2753 reference_area_m2,
2754 )?;
2755 terms.copies = self.aero.copies;
2756 terms.in_pod = true;
2757 terms.motor_pod_set = self.motor_pod_set;
2758 drag_terms.push(terms);
2759 self.terms_at.push((drag_terms.len() - 1, geometry));
2760 self.previous_aft_area = Some(geometry.aft_area_m2);
2761 self.aero
2762 .bodies
2763 .push(body_terms(component, geometry, step, reference_area_m2));
2764 Ok(())
2765 }
2766}
2767
2768fn pod_fins(
2776 set: &hpr_design::FinSet,
2777 component: &PlacedComponent,
2778 terms: &FinSetAero,
2779 reference_diameter_m: f64,
2780) -> Result<PodFins, AeroError> {
2781 use std::f64::consts::TAU;
2782 if set.cant_rad != 0.0 {
2783 return Err(AeroError::Unsupported(
2784 "cant on a pod's fins, whose roll forcing about the rocket's axis isn't modeled"
2785 .to_owned(),
2786 ));
2787 }
2788 let mut offsets = Vec::new();
2789 for copy in &component.copies {
2790 for j in 0..set.count {
2791 let angle =
2792 set.base_angle_rad + TAU * f64::from(j) / f64::from(set.count) + copy.roll_rad;
2793 let (sin, cos) = angle.sin_cos();
2794 offsets.push(copy.offset_m[0] * cos + copy.offset_m[1] * sin + terms.body_radius_m);
2795 }
2796 }
2797 let fins = u32::try_from(offsets.len()).map_err(|_| AeroError::Domain {
2798 what: "fins over a pod set's pods",
2799 value: offsets.len() as f64,
2800 })?;
2801 let (mean, deviation) = if offsets.is_empty() {
2802 (0.0, 0.0)
2803 } else {
2804 let n = f64::from(fins);
2805 let mean = offsets.iter().sum::<f64>() / n;
2806 let variance = offsets.iter().map(|o| (o - mean) * (o - mean)).sum::<f64>() / n;
2807 (mean, variance.sqrt())
2808 };
2809 Ok(PodFins {
2810 roll_rad: component.copies.iter().map(|c| c.roll_rad).collect(),
2811 axis_roll: [
2812 terms
2813 .fin
2814 .roll_terms_about(mean + deviation, reference_diameter_m)?,
2815 terms
2816 .fin
2817 .roll_terms_about(mean - deviation, reference_diameter_m)?,
2818 ],
2819 fins,
2820 })
2821}
2822
2823fn body_terms(
2824 component: &PlacedComponent,
2825 geometry: BodyGeometry,
2826 step_area_m2: f64,
2827 a_ref: f64,
2828) -> BodyAero {
2829 let station = component.fore_station_m;
2830 let step_slope = 2.0 * step_area_m2 / a_ref;
2831 let slope = geometry.normal_force_slope(a_ref);
2832 BodyAero {
2833 id: component.id.clone(),
2834 fore_station_m: station,
2835 geometry,
2836 step_area_m2,
2837 slope_per_rad: slope + step_slope,
2838 moment_slope_m: (slope + step_slope) * station + geometry.moment_slope_m(a_ref),
2839 planform_ratio: geometry.planform_area_m2 / a_ref,
2840 lift_station_m: station + geometry.planform_centroid_m,
2841 }
2842}
2843
2844#[cfg(test)]
2845mod tests {
2846 use std::f64::consts::{FRAC_PI_2, FRAC_PI_4};
2847
2848 use hpr_design::{
2849 FinPlanform, LaunchLug, NoseShape, Part, PodSet, Position, RailButton, ReferenceDiameter,
2850 TubeFinSet,
2851 };
2852 use proptest::prelude::*;
2853
2854 use super::*;
2855 use crate::BODY_LIFT_K;
2856 use crate::drag::base_drag_coefficient;
2857 use crate::testing::{body_part, component, fin_set, finned_rocket, material, nose, one_stage};
2858
2859 fn close(got: f64, want: f64, rel: f64, what: &str) {
2860 let err = if want == 0.0 {
2861 got.abs()
2862 } else {
2863 ((got - want) / want).abs()
2864 };
2865 assert!(
2866 err <= rel,
2867 "{what}: got {got}, want {want}, rel err {err:e}"
2868 );
2869 }
2870
2871 fn model(rocket: &hpr_design::Rocket) -> AeroModel {
2872 AeroModel::new(&rocket.layout().unwrap()).unwrap()
2873 }
2874
2875 fn podded_rocket(count: u32) -> hpr_design::Rocket {
2878 let mut rocket = finned_rocket(4);
2879 let mut pods = component(
2880 "pods",
2881 Part::PodSet(PodSet {
2882 count,
2883 radial_offset_m: 0.04,
2884 angle_rad: 0.3,
2885 }),
2886 Some(Position::Top { aft_offset_m: 0.1 }),
2887 );
2888 pods.children = vec![
2889 component("pod-nose", nose(NoseShape::Conical {}, 0.05, 0.01), None),
2890 component("pod-tube", body_part(0.2, 0.01, 0.01), None),
2891 ];
2892 rocket.stages[0].components[1].children.push(pods);
2893 rocket
2894 }
2895
2896 fn pod_fin_planform() -> FinPlanform {
2898 FinPlanform::Trapezoidal {
2899 root_chord_m: 0.06,
2900 tip_chord_m: 0.03,
2901 span_m: 0.04,
2902 sweep_m: 0.02,
2903 }
2904 }
2905
2906 fn winglet_rocket(fins: u32, base_angle_rad: f64) -> hpr_design::Rocket {
2910 let mut rocket = finned_rocket(4);
2911 let mut pods = component(
2912 "pods",
2913 Part::PodSet(PodSet {
2914 count: 2,
2915 radial_offset_m: 0.05,
2916 angle_rad: 0.0,
2917 }),
2918 Some(Position::Top { aft_offset_m: 0.3 }),
2919 );
2920 let mut phantom = component("phantom", body_part(0.0, 0.0, 0.0), None);
2921 let mut set = component(
2922 "winglets",
2923 fin_set(fins, pod_fin_planform()),
2924 Some(Position::Top { aft_offset_m: 0.0 }),
2925 );
2926 if let Part::FinSet(fin_set) = &mut set.part {
2927 fin_set.base_angle_rad = base_angle_rad;
2928 }
2929 phantom.children = vec![set];
2930 pods.children = vec![phantom];
2931 rocket.stages[0].components[1].children.push(pods);
2932 rocket
2933 }
2934
2935 #[test]
2940 fn a_pod_adds_its_bodies_slopes_once_per_pod() {
2941 let bare = model(&finned_rocket(4));
2942 let area_ratio = (0.01_f64 / 0.027).powi(2);
2943 for count in [1, 2, 3] {
2944 let m = model(&podded_rocket(count));
2945 assert_eq!(m.pod_sets().len(), 1);
2946 let pod = &m.pod_sets()[0];
2947 assert_eq!(pod.copies, count);
2948 close(pod.fineness, 0.25 / 0.02, 1e-15, "pod fineness");
2949 close(
2950 pod.offset_squares_m2,
2951 f64::from(count) * 0.04 * 0.04,
2952 1e-14,
2953 "offsets",
2954 );
2955 let n = f64::from(count);
2956 for mach in [0.3, 0.95, 2.0, 3.5] {
2957 let what = format!("{count} pods, Mach {mach}");
2958 let parts = m.components(&Flow::axial(mach)).unwrap();
2959 let find = |id: &str| {
2960 parts
2961 .iter()
2962 .find(|c| c.id == id)
2963 .map(|c| c.normal_force)
2964 .unwrap()
2965 };
2966 let cone = find("pod-nose");
2967 close(cone.slope_per_rad, n * 2.0 * area_ratio, 1e-12, &what);
2968 close(
2970 cone.cp_station_m.unwrap(),
2971 0.35 + 2.0 / 3.0 * 0.05,
2972 1e-12,
2973 &what,
2974 );
2975 close(find("pod-tube").slope_per_rad, 0.0, 1e-15, &what);
2976 let with = m.normal_force(&Flow::axial(mach)).unwrap();
2977 let without = bare.normal_force(&Flow::axial(mach)).unwrap();
2978 close(
2979 with.slope_per_rad - without.slope_per_rad,
2980 n * 2.0 * area_ratio,
2981 1e-11,
2982 &what,
2983 );
2984 close(
2985 with.moment_slope_m - without.moment_slope_m,
2986 n * 2.0 * area_ratio * (0.35 + 2.0 / 3.0 * 0.05),
2987 1e-11,
2988 &what,
2989 );
2990 let index = parts.iter().position(|c| c.id == "pod-nose").unwrap();
2992 assert_eq!(index, bare.bodies().len());
2993 assert_eq!(m.fin_set_start(), bare.bodies().len() + 2);
2994 close(
2995 m.component_station_m(index, mach).unwrap(),
2996 0.35 + 2.0 / 3.0 * 0.05,
2997 1e-12,
2998 &what,
2999 );
3000 let one = m.component_normal_force(index, &Flow::axial(mach)).unwrap();
3001 assert_eq!(one, cone);
3002 }
3003 }
3004 }
3005
3006 #[test]
3009 fn a_pod_s_body_lift_takes_its_own_fineness() {
3010 let m = model(&podded_rocket(3));
3011 let pod = &m.pod_sets()[0];
3012 let (mach, alpha) = (0.5, 0.2_f64);
3013 let parts = m.components(&flow(mach, alpha, 0.0)).unwrap();
3014 let tube = parts.iter().find(|c| c.id == "pod-tube").unwrap();
3015 let factor = m.lift_factor_of(
3016 12.5,
3017 crate::crossflow::crossflow_eta_low(12.5),
3018 mach,
3019 alpha.sin(),
3020 );
3021 assert!(
3022 (factor - m.lift_factor_of(m.fineness(), m.crossflow_eta_low, mach, alpha.sin())).abs()
3023 > 1e-3,
3024 "the pod's fineness must matter"
3025 );
3026 let planform = 2.0 * 0.01 * 0.2 / m.reference_area_m2();
3027 close(
3028 tube.normal_force.coefficient,
3029 3.0 * factor * planform * alpha.sin().powi(2),
3030 1e-12,
3031 "pod tube's body lift",
3032 );
3033 close(pod.fineness, 12.5, 1e-15, "fineness");
3034 }
3035
3036 #[test]
3040 fn a_pod_drags_once_per_pod() {
3041 let (mach, reynolds_per_m) = (0.4, 5e6);
3042 for count in [1_u32, 2, 3] {
3043 let m = model(&podded_rocket(count));
3044 let n = f64::from(count);
3045 let a_ref = m.reference_area_m2();
3046 for conditions in [
3047 DragConditions::coasting(reynolds_per_m),
3048 DragConditions::thrusting(reynolds_per_m, 2e-4),
3049 ] {
3050 let parts = m
3051 .buildup_components(&Flow::axial(mach), &conditions)
3052 .unwrap();
3053 let find = |id: &str| parts.iter().find(|c| c.id == id).unwrap().drag;
3054 let (body, tube) = (find("body"), find("pod-tube"));
3055 let ratio = n * body_friction_form_factor(12.5).unwrap()
3057 / body_friction_form_factor(m.fineness()).unwrap()
3058 * (0.01 * 0.2)
3059 / (0.027 * 0.7);
3060 close(tube.friction, ratio * body.friction, 1e-12, "pod friction");
3061 close(
3062 tube.base,
3063 n * base_drag_coefficient(mach).unwrap() * PI * 0.01 * 0.01 / a_ref,
3064 1e-12,
3065 "pod base",
3066 );
3067 let cone = find("pod-nose");
3069 assert_eq!(cone.base, 0.0);
3070 assert!(cone.pressure > 0.0);
3071 let one = model(&podded_rocket(1));
3072 let single = one
3073 .buildup_components(&Flow::axial(mach), &conditions)
3074 .unwrap();
3075 let single_cone = single.iter().find(|c| c.id == "pod-nose").unwrap().drag;
3076 close(cone.pressure, n * single_cone.pressure, 1e-14, "pod nose");
3077 let total = m.drag(&Flow::axial(mach), &conditions).unwrap();
3078 let sum: f64 = parts.iter().map(|c| c.drag.zero_lift_coefficient).sum();
3079 close(total.zero_lift_coefficient, sum, 1e-14, "total");
3080 }
3081 }
3082 }
3083
3084 #[test]
3089 fn a_pod_s_base_takes_its_own_motors_area() {
3090 let (mach, reynolds_per_m) = (0.4, 5e6);
3091 let mut rocket = podded_rocket(2);
3092 let pods = rocket.stages[0].components[1].children.last_mut().unwrap();
3093 pods.children[1].motor_mount = Some(hpr_design::MotorMount::default());
3094 let m = model(&rocket);
3095 let a_ref = m.reference_area_m2();
3096 let base = |conditions: &DragConditions, id: &str| {
3097 m.buildup_components(&Flow::axial(mach), conditions)
3098 .unwrap()
3099 .into_iter()
3100 .find(|c| c.id == id)
3101 .unwrap()
3102 .drag
3103 .base
3104 };
3105 let pod_area = PI * 0.01 * 0.01;
3106 let c_b = base_drag_coefficient(mach).unwrap();
3107 let coasting = DragConditions::coasting(reynolds_per_m);
3108 close(
3109 base(&coasting, "pod-tube"),
3110 2.0 * c_b * pod_area / a_ref,
3111 1e-14,
3112 "coasting",
3113 );
3114 let motor = PI * 0.008 * 0.008;
3116 let thrusting = DragConditions::thrusting(reynolds_per_m, 0.0).with_pod_motors([
3117 2.0 * motor,
3118 0.0,
3119 0.0,
3120 0.0,
3121 ]);
3122 close(
3123 base(&thrusting, "pod-tube"),
3124 2.0 * c_b * (pod_area - motor) / a_ref,
3125 1e-14,
3126 "pod motors",
3127 );
3128 assert_eq!(base(&thrusting, "tail"), base(&coasting, "tail"));
3129 let core = DragConditions::thrusting(reynolds_per_m, 1e-4);
3131 assert_eq!(base(&core, "pod-tube"), base(&coasting, "pod-tube"));
3132 assert!(base(&core, "tail") < base(&coasting, "tail"));
3133 let mut second = component(
3135 "more-pods",
3136 Part::PodSet(PodSet {
3137 count: 2,
3138 radial_offset_m: 0.04,
3139 angle_rad: 1.8,
3140 }),
3141 Some(Position::Top { aft_offset_m: 0.4 }),
3142 );
3143 let mut tube = component("more-pod-tube", body_part(0.1, 0.01, 0.01), None);
3144 tube.motor_mount = Some(hpr_design::MotorMount::default());
3145 second.children = vec![tube];
3146 rocket.stages[0].components[1].children.push(second);
3147 let two = model(&rocket);
3148 let base = |conditions: &DragConditions, id: &str| {
3149 two.buildup_components(&Flow::axial(mach), conditions)
3150 .unwrap()
3151 .into_iter()
3152 .find(|c| c.id == id)
3153 .unwrap()
3154 .drag
3155 .base
3156 };
3157 let motor = PI * 0.004 * 0.004;
3158 let first = DragConditions::thrusting(reynolds_per_m, 0.0).with_pod_motors([
3159 2.0 * motor,
3160 0.0,
3161 0.0,
3162 0.0,
3163 ]);
3164 let second = DragConditions::thrusting(reynolds_per_m, 0.0).with_pod_motors([
3165 0.0,
3166 2.0 * motor,
3167 0.0,
3168 0.0,
3169 ]);
3170 for (conditions, relieved, whole) in [
3171 (&first, "pod-tube", "more-pod-tube"),
3172 (&second, "more-pod-tube", "pod-tube"),
3173 ] {
3174 close(
3175 base(conditions, relieved),
3176 2.0 * c_b * (pod_area - motor) / a_ref,
3177 1e-14,
3178 relieved,
3179 );
3180 assert_eq!(base(conditions, whole), base(&coasting, whole), "{whole}");
3181 }
3182 }
3183
3184 #[test]
3188 fn a_whole_base_under_power_keeps_the_motors_area() {
3189 let mach = 0.4;
3190 let mut rocket = podded_rocket(2);
3191 let pods = rocket.stages[0].components[1].children.last_mut().unwrap();
3192 pods.children[1].motor_mount = Some(hpr_design::MotorMount::default());
3193 let relieved = model(&rocket);
3194 assert!(!relieved.full_base_drag_under_power());
3195 let whole = relieved.clone().with_full_base_drag_under_power();
3196 assert!(whole.full_base_drag_under_power());
3197 let flow = Flow::axial(mach);
3198 let coasting = DragConditions::coasting(5e6);
3199 let thrusting = DragConditions::thrusting(5e6, 1e-4).with_pod_motors([2e-4, 0.0, 0.0, 0.0]);
3200 let base =
3201 |m: &AeroModel, conditions: &DragConditions| m.drag(&flow, conditions).unwrap().base;
3202 assert_eq!(base(&whole, &thrusting), base(&whole, &coasting));
3203 assert_eq!(base(&whole, &coasting), base(&relieved, &coasting));
3204 assert!(base(&relieved, &thrusting) < base(&relieved, &coasting));
3205 for id in ["tail", "pod-tube"] {
3206 let part = |m: &AeroModel, conditions: &DragConditions| {
3207 m.buildup_components(&flow, conditions)
3208 .unwrap()
3209 .into_iter()
3210 .find(|c| c.id == id)
3211 .unwrap()
3212 .drag
3213 .base
3214 };
3215 assert_eq!(part(&whole, &thrusting), part(&whole, &coasting), "{id}");
3216 assert!(
3217 part(&relieved, &thrusting) < part(&relieved, &coasting),
3218 "{id}"
3219 );
3220 }
3221 let bad = DragConditions::thrusting(5e6, -1e-4);
3223 assert!(whole.drag(&flow, &bad).is_err());
3224 assert!(whole.buildup_components(&flow, &bad).is_err());
3225 }
3226
3227 #[test]
3231 fn a_pod_s_fins_turn_with_their_pod() {
3232 let bare = model(&finned_rocket(4));
3233 let mut rocket = winglet_rocket(1, 0.0);
3234 if let Part::PodSet(pods) = &mut rocket.stages[0].components[1]
3235 .children
3236 .last_mut()
3237 .unwrap()
3238 .part
3239 {
3240 pods.angle_rad = FRAC_PI_4;
3241 }
3242 let m = model(&rocket);
3243 let mach = 0.3;
3244 let slope = FinAero::new(&pod_fin_planform(), m.reference_area_m2())
3245 .unwrap()
3246 .loading(mach)
3247 .unwrap()
3248 .slope_per_rad;
3249 for (roll, share) in [(FRAC_PI_4, 0.0), (-FRAC_PI_4, 2.0)] {
3250 let with = m.normal_force(&flow(mach, 0.0, roll)).unwrap();
3251 let without = bare.normal_force(&flow(mach, 0.0, roll)).unwrap();
3252 close(
3253 with.slope_per_rad - without.slope_per_rad,
3254 share * slope,
3255 1e-12,
3256 &format!("roll {roll}"),
3257 );
3258 }
3259 }
3260
3261 #[test]
3265 fn every_component_index_maps_to_its_own_terms_with_pods() {
3266 let mut rocket = podded_rocket(3);
3267 let winglets = winglet_rocket(2, 0.3);
3268 let mut pods = winglets.stages[0].components[1]
3269 .children
3270 .last()
3271 .unwrap()
3272 .clone();
3273 pods.id = "winglet-pods".to_owned();
3274 rocket.stages[0].components[1].children.push(pods);
3275 let mut second = component(
3276 "aft-pods",
3277 Part::PodSet(PodSet {
3278 count: 2,
3279 radial_offset_m: 0.05,
3280 angle_rad: 1.0,
3281 }),
3282 Some(Position::Top { aft_offset_m: 0.45 }),
3283 );
3284 second.children = vec![
3285 component(
3286 "aft-pod-nose",
3287 nose(NoseShape::Conical {}, 0.04, 0.012),
3288 None,
3289 ),
3290 component("aft-pod-tube", body_part(0.1, 0.012, 0.012), None),
3291 component("aft-pod-tail", body_part(0.03, 0.012, 0.008), None),
3292 ];
3293 rocket.stages[0].components[1].children.push(second);
3294 let m = model(&rocket);
3295 assert_eq!(m.pod_sets().len(), 2);
3296 assert_eq!(m.fin_set_start(), m.bodies().len() + 5);
3297 assert_eq!(m.component_count(), m.fin_set_start() + 2);
3298 for mach in [0.3, 2.0] {
3299 let f = flow(mach, 0.1, 0.4);
3300 let parts = m.components(&f).unwrap();
3301 assert_eq!(parts.len(), m.component_count());
3302 for (index, part) in parts.iter().enumerate() {
3303 let one = m.component_normal_force(index, &f).unwrap();
3304 assert_eq!(one, part.normal_force, "{} at {index}", part.id);
3305 let station = m.component_station_m(index, mach).unwrap();
3306 if let Some(cp) = m.components(&flow(mach, 0.0, 0.4)).unwrap()[index]
3307 .normal_force
3308 .cp_station_m
3309 && index >= m.bodies().len()
3310 {
3311 close(station, cp, 1e-12, &part.id);
3312 }
3313 }
3314 let ids: Vec<&str> = parts.iter().map(|p| p.id.as_str()).collect();
3315 let at = |id: &str| ids.iter().position(|i| *i == id).unwrap();
3316 assert!(at("pod-tube") < at("aft-pod-nose") && at("aft-pod-tail") < at("winglets"));
3317 assert!(m.component_normal_force(m.component_count(), &f).is_err());
3318 }
3319 }
3320
3321 #[test]
3324 fn the_aero_page_s_pod_example() {
3325 let (bare, pods) = (model(&finned_rocket(4)), model(&podded_rocket(3)));
3326 let at = |m: &AeroModel| {
3327 let n = m.normal_force(&Flow::axial(0.3)).unwrap();
3328 let c = DragConditions::coasting(5e6);
3329 let d = m.drag(&Flow::axial(0.3), &c).unwrap();
3330 (
3331 n.slope_per_rad,
3332 n.cp_station_m.unwrap(),
3333 d.zero_lift_coefficient,
3334 m.roll(0.3).unwrap().damping,
3335 )
3336 };
3337 let round = |x: f64, digits: i32| (x * 10_f64.powi(digits)).round() / 10_f64.powi(digits);
3338 let (slope, cp, drag, damping) = at(&bare);
3339 assert_eq!(
3340 [
3341 round(slope, 3),
3342 round(cp, 4),
3343 round(drag, 4),
3344 round(damping, 3)
3345 ],
3346 [12.374, 1.0662, 0.5051, -35.215]
3347 );
3348 let (slope, cp, drag, damping) = at(&pods);
3349 assert_eq!(
3350 [
3351 round(slope, 3),
3352 round(cp, 4),
3353 round(drag, 4),
3354 round(damping, 3)
3355 ],
3356 [13.197, 1.0237, 0.6386, -36.118]
3357 );
3358 let parts = pods
3359 .buildup_components(&Flow::axial(0.3), &DragConditions::coasting(5e6))
3360 .unwrap();
3361 let find = |id: &str| parts.iter().find(|c| c.id == id).unwrap().drag;
3362 let (cone, tube) = (find("pod-nose"), find("pod-tube"));
3363 assert_eq!(
3364 [
3365 round(cone.friction, 4),
3366 round(cone.pressure, 4),
3367 round(tube.friction, 4),
3368 round(tube.base, 4)
3369 ],
3370 [0.0074, 0.0127, 0.0592, 0.0542]
3371 );
3372 let single = model(&podded_rocket(1))
3373 .buildup_components(&Flow::axial(0.3), &DragConditions::coasting(5e6))
3374 .unwrap();
3375 let one_pod: f64 = single
3376 .iter()
3377 .filter(|c| c.id.starts_with("pod"))
3378 .map(|c| c.drag.zero_lift_coefficient)
3379 .sum();
3380 assert_eq!(round(one_pod, 4), 0.0445);
3381 let one_slope = model(&podded_rocket(1))
3382 .normal_force(&Flow::axial(0.3))
3383 .unwrap()
3384 .slope_per_rad;
3385 assert_eq!(round(one_slope, 2), 12.65);
3386 }
3387
3388 #[test]
3393 fn a_pod_s_bodies_damp_the_roll() {
3394 let bare = model(&finned_rocket(4));
3395 let m = model(&podded_rocket(3));
3396 let slope = 2.0 * (0.01_f64 / 0.027).powi(2);
3397 let want = -2.0 * slope * 3.0 * 0.04 * 0.04 / (0.054 * 0.054);
3398 close(want, -0.9032_f64, 1e-4, "worked number");
3399 for mach in [0.3, 1.5] {
3400 let got = m.roll(mach).unwrap().damping - bare.roll(mach).unwrap().damping;
3401 close(got, want, 1e-12, "pod roll damping");
3402 assert_eq!(
3403 m.roll(mach).unwrap().forcing,
3404 bare.roll(mach).unwrap().forcing
3405 );
3406 }
3407 }
3408
3409 #[test]
3413 fn a_pod_s_fins_are_the_pod_s_turned_with_it() {
3414 let bare = model(&finned_rocket(4));
3415 let a_ref = bare.reference_area_m2();
3416 let d = 0.054;
3417 let fin = FinAero::new(&pod_fin_planform(), a_ref).unwrap();
3418 let m = model(&winglet_rocket(1, 0.0));
3421 assert!(m.pod_sets().is_empty(), "a phantom body has no body terms");
3422 for mach in [0.3, 1.8] {
3423 let slope = fin.loading(mach).unwrap().slope_per_rad;
3424 let across = m.normal_force(&flow(mach, 0.0, FRAC_PI_2)).unwrap();
3425 let bare_across = bare.normal_force(&flow(mach, 0.0, FRAC_PI_2)).unwrap();
3426 close(
3427 across.slope_per_rad - bare_across.slope_per_rad,
3428 2.0 * slope,
3429 1e-12,
3430 "across",
3431 );
3432 let along = m.normal_force(&flow(mach, 0.0, 0.0)).unwrap();
3433 let bare_along = bare.normal_force(&flow(mach, 0.0, 0.0)).unwrap();
3434 close(
3435 along.slope_per_rad - bare_along.slope_per_rad,
3436 0.0,
3437 1e-12,
3438 "along",
3439 );
3440 let damping = m.roll(mach).unwrap().damping - bare.roll(mach).unwrap().damping;
3443 close(
3444 damping,
3445 2.0 * fin.roll(mach, 0.05, d).unwrap().damping,
3446 1e-12,
3447 "damping",
3448 );
3449 }
3450 let m = model(&winglet_rocket(2, 0.0));
3453 let mach = 0.5;
3454 let (c_r, c_t, s) = (0.06, 0.03, 0.04);
3455 let area = 0.5 * s * (c_r + c_t);
3456 let per_area = fin.loading(mach).unwrap().slope_per_rad * a_ref / area;
3457 let strips = 20_000;
3458 let mut second = 0.0;
3459 for i in 0..strips {
3460 let y = s * (f64::from(i) + 0.5) / f64::from(strips);
3461 let chord = c_r + (c_t - c_r) * y / s;
3462 for root in [0.05, -0.05] {
3463 second += 2.0 * (root + y) * (root + y) * chord * s / f64::from(strips);
3464 }
3465 }
3466 let want = -2.0 * per_area * second / (a_ref * d * d);
3467 let got = m.roll(mach).unwrap().damping - bare.roll(mach).unwrap().damping;
3468 close(got, want, 1e-7, "strips");
3469 let set = m.fin_sets().iter().find(|s| s.id == "winglets").unwrap();
3470 let pods = set.pods.as_ref().unwrap();
3471 assert_eq!(pods.fins, 4);
3472 assert_eq!(pods.roll_rad.len(), 2);
3473 }
3474
3475 fn flow(mach: f64, alpha_rad: f64, roll_rad: f64) -> Flow {
3476 Flow::new(mach, alpha_rad, roll_rad)
3477 }
3478
3479 #[test]
3484 fn angle_of_attack_terms() {
3485 let (l_n, l_t, r) = (0.2, 0.8, 0.03);
3486 let rocket = one_stage(
3487 vec![
3488 component("nose", nose(NoseShape::Conical {}, l_n, r), None),
3489 component("tube", body_part(l_t, r, r), None),
3490 ],
3491 ReferenceDiameter::Maximum {},
3492 );
3493 let m = model(&rocket);
3494 let a_ref = PI * r * r;
3495 let k = |alpha: f64| crate::crossflow::crossflow_factor(1.0 / 0.06, 0.3 * alpha.sin());
3497 close(m.fineness(), 1.0 / 0.06, 1e-15, "fineness");
3498 let lift_nose = k(FRAC_PI_2) * r * l_n / a_ref;
3499 let lift_tube = k(FRAC_PI_2) * 2.0 * r * l_t / a_ref;
3500
3501 let broadside = m.normal_force(&flow(0.3, FRAC_PI_2, 0.0)).unwrap();
3502 let want = 2.0 + lift_nose + lift_tube;
3503 close(broadside.coefficient, want, 1e-10, "C_N at 90°");
3504 let moment =
3505 2.0 * (2.0 * l_n / 3.0) + lift_nose * (2.0 * l_n / 3.0) + lift_tube * (l_n + 0.5 * l_t);
3506 close(
3507 broadside.cp_station_m.unwrap(),
3508 moment / want,
3509 1e-10,
3510 "CP at 90°",
3511 );
3512
3513 let zero = m.normal_force(&Flow::axial(0.3)).unwrap();
3514 assert_eq!(zero.coefficient, 0.0);
3515 close(zero.slope_per_rad, 2.0, 1e-15, "slope at 0");
3516 for alpha in [1e-6, 1e-3, 0.05] {
3518 let f = m.normal_force(&flow(0.3, alpha, 0.0)).unwrap();
3519 let lift = (lift_nose + lift_tube) * k(alpha) / k(FRAC_PI_2);
3520 let want = 2.0 * alpha.sin() + lift * alpha.sin().powi(2);
3521 close(f.coefficient, want, 1e-12, "C_N");
3522 close(f.slope_per_rad, want / alpha, 1e-12, "C_N/α");
3523 }
3524 let cp = |alpha| {
3526 m.normal_force(&flow(0.3, alpha, 0.0))
3527 .unwrap()
3528 .cp_station_m
3529 .unwrap()
3530 };
3531 assert!(cp(0.02) > cp(0.0) && cp(0.2) > cp(0.02));
3532 let parts = m.components(&flow(0.3, 0.2, 0.0)).unwrap();
3534 let total = m.normal_force(&flow(0.3, 0.2, 0.0)).unwrap();
3535 let sum: f64 = parts.iter().map(|c| c.normal_force.coefficient).sum();
3536 close(sum, total.coefficient, 1e-14, "component sum");
3537 }
3538
3539 #[test]
3542 fn body_lift_models() {
3543 let rocket = one_stage(
3544 vec![
3545 component("nose", nose(NoseShape::Conical {}, 0.2, 0.03), None),
3546 component("tube", body_part(0.8, 0.03, 0.03), None),
3547 ],
3548 ReferenceDiameter::Maximum {},
3549 );
3550 let layout = rocket.layout().unwrap();
3551 let old = AeroModel::with_body_model(&layout, BodyModel::BEFORE_M1_8E6).unwrap();
3552 let new = AeroModel::new(&layout).unwrap();
3553 assert_eq!(new.body_model(), BodyModel::default());
3554 for (mach, alpha) in [(0.3, 0.1), (2.0, 0.5), (4.0, 1.2)] {
3555 assert_eq!(old.body_lift_factor(&flow(mach, alpha, 0.0)), BODY_LIFT_K);
3556 }
3557 let slow = new.body_lift_factor(&flow(0.3, 0.1, 0.0));
3558 assert!(slow > 0.85 && slow < 0.9, "{slow}");
3559 let near_one = new.body_lift_factor(&flow(2.0, 0.5, 0.0));
3560 assert!(near_one > 1.4, "{near_one}");
3561 let k = BodyModel::BEFORE_M1_8E6.with_body_lift(BodyLift::Galejs { k: -0.5 });
3562 assert!(AeroModel::with_body_model(&layout, k).is_err());
3563 for (mach, alpha) in [(0.3, 0.1), (2.0, 0.5), (4.0, 1.2)] {
3565 let f = flow(mach, alpha, 0.0);
3566 assert_eq!(
3567 new.body_lift_factor(&f),
3568 crate::crossflow::crossflow_factor(new.fineness(), mach * alpha.sin())
3569 );
3570 }
3571 }
3572
3573 #[test]
3576 fn body_model_in_json() {
3577 assert_eq!(BodyModel::default(), BodyModel::CURRENT);
3578 let old = serde_json::to_string(&BodyModel::BEFORE_M1_8E6).unwrap();
3579 assert_eq!(
3580 old,
3581 r#"{"body_lift":{"kind":"galejs","k":1.1},"supersonic_boattail":"footnote8","supersonic_flare":"slender_body"}"#
3582 );
3583 let before_m1_8e17 =
3587 r#"{"body_lift":{"kind":"galejs","k":1.1},"supersonic_boattail":"footnote8"}"#;
3588 assert_eq!(
3589 serde_json::from_str::<BodyModel>(before_m1_8e17)
3590 .unwrap()
3591 .supersonic_flare,
3592 SupersonicFlare::Marched
3593 );
3594 assert_eq!(
3595 serde_json::from_str::<BodyModel>(&old).unwrap(),
3596 BodyModel::BEFORE_M1_8E6
3597 );
3598 assert_eq!(
3599 serde_json::from_str::<BodyModel>(r#"{"supersonic_boattail":"washington_pettis"}"#)
3600 .unwrap(),
3601 BodyModel::CURRENT
3602 );
3603 assert_eq!(
3604 serde_json::from_str::<BodyModel>("{}").unwrap(),
3605 BodyModel::CURRENT
3606 );
3607 for bad in [
3608 r#"{"boattail":"footnote8"}"#,
3609 r#"{"supersonic_boattail":"slender_body"}"#,
3610 ] {
3611 assert!(serde_json::from_str::<BodyModel>(bad).is_err(), "{bad}");
3612 }
3613 }
3614
3615 #[test]
3618 fn mach_changes_only_the_fins() {
3619 let m = model(&finned_rocket(4));
3620 let at = |mach| m.components(&Flow::axial(mach)).unwrap();
3621 let (slow, fast) = (at(0.0), at(0.8));
3622 for (a, b) in slow.iter().zip(&fast) {
3623 assert_eq!(
3624 a.normal_force.cp_station_m, b.normal_force.cp_station_m,
3625 "{}",
3626 a.id
3627 );
3628 if a.id == "fins" {
3629 let set = &m.fin_sets()[0];
3630 let ratio = set
3631 .fin
3632 .geometry()
3633 .single_fin_slope(m.reference_area_m2(), 0.8)
3634 .unwrap()
3635 / set
3636 .fin
3637 .geometry()
3638 .single_fin_slope(m.reference_area_m2(), 0.0)
3639 .unwrap();
3640 assert!(ratio > 1.05, "{ratio}");
3641 close(
3642 b.normal_force.slope_per_rad / a.normal_force.slope_per_rad,
3643 ratio,
3644 1e-14,
3645 "fin ratio",
3646 );
3647 } else {
3648 assert_eq!(
3649 a.normal_force.slope_per_rad, b.normal_force.slope_per_rad,
3650 "{}",
3651 a.id
3652 );
3653 }
3654 }
3655 }
3656
3657 #[test]
3659 fn two_fin_sets_depend_on_roll() {
3660 let four = model(&finned_rocket(4));
3661 let slope =
3662 |m: &AeroModel, roll| m.normal_force(&flow(0.2, 0.0, roll)).unwrap().slope_per_rad;
3663 close(slope(&four, 0.0), slope(&four, 0.4), 1e-15, "four fins");
3664
3665 let two = model(&finned_rocket(2));
3666 let bodies: f64 = two.bodies().iter().map(|b| b.slope_per_rad).sum();
3667 let set = &two.fin_sets()[0];
3668 let one_fin = set
3669 .fin
3670 .geometry()
3671 .single_fin_slope(two.reference_area_m2(), 0.2)
3672 .unwrap()
3673 * set.interference;
3674 close(slope(&two, 0.0), bodies, 1e-13, "along the fins");
3676 close(
3677 slope(&two, FRAC_PI_2),
3678 bodies + 2.0 * one_fin,
3679 1e-13,
3680 "across the fins",
3681 );
3682 close(
3684 slope(&two, FRAC_PI_2),
3685 slope(&four, 0.0),
3686 1e-13,
3687 "two across = four",
3688 );
3689 }
3690
3691 fn tube_finned_rocket(
3694 count: u32,
3695 length_m: f64,
3696 outer_radius_m: f64,
3697 thickness_m: f64,
3698 ) -> hpr_design::Rocket {
3699 let mut rocket = crate::testing::finned_rocket(4);
3700 rocket.stages[0].components[3].children.push(component(
3701 "tube-fins",
3702 Part::TubeFinSet(TubeFinSet {
3703 count,
3704 length_m,
3705 outer_radius_m,
3706 thickness_m,
3707 base_angle_rad: 0.3,
3708 material: material(),
3709 }),
3710 Some(Position::Bottom { aft_offset_m: 0.0 }),
3711 ));
3712 rocket
3713 }
3714
3715 #[test]
3716 fn tube_fins_add_their_rings_slopes_at_fletcher_s_center() {
3717 let plain = model(&crate::testing::finned_rocket(4));
3718 let (count, length, outer, wall) = (6_u32, 0.1, 0.022, 0.0005);
3719 let m = model(&tube_finned_rocket(count, length, outer, wall));
3720 let a_ref = m.reference_area_m2();
3721 let d = 2.0 * outer - wall;
3722 let fore = 1.3 - length;
3724 assert_eq!(m.component_count(), plain.component_count() + 1);
3725 assert_eq!(m.tube_fin_set_start(), plain.component_count());
3726 for (mach, roll) in [(0.0_f64, 0.0), (0.3, 0.4), (0.7, 1.1)] {
3727 let beta = (1.0 - mach * mach).sqrt();
3728 let lambda = length / d;
3729 let slope = f64::from(count)
3731 * (PI * PI
3732 / (1.0
3733 + 0.5 * PI * lambda / beta
3734 + lambda / beta * (1.2 * lambda / beta).atan()))
3735 / beta
3736 * d
3737 * length
3738 / a_ref;
3739 let a = beta * d / length;
3741 assert!(a < 2.0 / 3.0);
3742 let center = 0.143 * a / (2.0 / 3.0);
3743 let station = fore + center * length;
3744 let flow = Flow::new(mach, 0.05, roll);
3745 let (with, without) = (
3746 m.normal_force(&flow).unwrap(),
3747 plain.normal_force(&flow).unwrap(),
3748 );
3749 close(
3750 with.slope_per_rad - without.slope_per_rad,
3751 slope,
3752 1e-12,
3753 "slope",
3754 );
3755 close(
3756 with.moment_slope_m - without.moment_slope_m,
3757 slope * station,
3758 1e-12,
3759 "moment",
3760 );
3761 close(
3763 with.side_coefficient,
3764 without.side_coefficient,
3765 1e-12,
3766 "side",
3767 );
3768 let index = m.tube_fin_set_start();
3769 let own = m.component_normal_force(index, &flow).unwrap();
3770 close(own.slope_per_rad, slope, 1e-12, "component slope");
3771 close(
3772 m.component_station_m(index, mach).unwrap(),
3773 station,
3774 1e-12,
3775 "station",
3776 );
3777 let listed = m.components(&flow).unwrap();
3778 assert_eq!(listed[index].id, "tube-fins");
3779 close(
3780 listed[index].normal_force.slope_per_rad,
3781 slope,
3782 1e-12,
3783 "listed",
3784 );
3785 let rho = 0.022 + outer;
3787 let damping = -2.0 * slope * rho * rho / (0.054 * 0.054);
3788 close(
3789 m.roll(mach).unwrap().damping - plain.roll(mach).unwrap().damping,
3790 damping,
3791 1e-12,
3792 "roll damping",
3793 );
3794 assert_eq!(m.roll(mach).unwrap().forcing, 0.0);
3795 }
3796 }
3797
3798 #[test]
3799 fn short_tube_fins_take_fletcher_s_measured_center() {
3800 let (length, outer, wall) = (0.02, 0.022, 0.0005);
3805 let m = model(&tube_finned_rocket(6, length, outer, wall));
3806 let d = 2.0 * outer - wall;
3807 let index = m.tube_fin_set_start();
3808 for (mach, a, by_hand) in [(0.0, d / length, 5.0510), (0.6, 0.8 * d / length, 5.4837)] {
3809 let fraction = 0.253 + (0.355 - 0.253) * (a - 1.5) / 1.5;
3810 let station = 1.3 - length + fraction * length;
3811 close(
3812 m.component_station_m(index, mach).unwrap(),
3813 station,
3814 1e-12,
3815 "station",
3816 );
3817 let slope = m
3818 .component_normal_force(index, &Flow::new(mach, 0.05, 0.0))
3819 .unwrap()
3820 .slope_per_rad;
3821 let per_ring = slope / 6.0 * m.reference_area_m2() / (d * length);
3822 close(per_ring, by_hand, 1e-4, "slope on d L");
3823 }
3824 }
3825
3826 #[test]
3830 fn the_guide_s_tube_fin_example() {
3831 let r = 0.012_395_2;
3832 let mut tail = component("tail", body_part(0.4572, r, r), None);
3833 tail.children = vec![component(
3834 "tube-fins",
3835 Part::TubeFinSet(TubeFinSet {
3836 count: 6,
3837 length_m: 0.0762,
3838 outer_radius_m: r,
3839 thickness_m: 0.000_330_2,
3840 base_angle_rad: 0.0,
3841 material: material(),
3842 }),
3843 Some(Position::Bottom { aft_offset_m: 0.0 }),
3844 )];
3845 let rocket = one_stage(
3846 vec![
3847 component(
3848 "nose",
3849 nose(NoseShape::Ogive { radius_ratio: 1.0 }, 0.119_888, r),
3850 None,
3851 ),
3852 tail,
3853 ],
3854 hpr_design::ReferenceDiameter::Maximum {},
3855 );
3856 let m = model(&rocket);
3857 close(m.reference_area_m2(), 4.827e-4, 1e-4, "reference area");
3858 let set = &m.tube_fin_sets()[0];
3859 close(set.mean_diameter_m, 0.024_460_2, 1e-12, "d");
3860 close(set.length_m / set.mean_diameter_m, 3.115, 1e-3, "λ");
3861 let (slope, station) = set.loading(0.0).unwrap();
3862 close(slope, 22.93, 1e-3, "slope at Mach 0");
3863 close(
3864 set.loading(0.35).unwrap().0,
3865 22.96,
3866 1e-3,
3867 "slope at Mach 0.35",
3868 );
3869 close(
3870 (station - set.fore_station_m) / set.length_m,
3871 0.0689,
3872 1e-3,
3873 "center",
3874 );
3875 let terms = m.drag_terms.iter().find(|t| t.id == "tube-fins").unwrap();
3876 close(terms.friction_area_ratio, 145.6, 1e-3, "friction area");
3877 close(
3878 terms.fins.as_ref().unwrap().frontal_area_ratio,
3879 0.3154,
3880 1e-3,
3881 "wall area",
3882 );
3883 }
3884
3885 #[test]
3886 fn tube_fins_drag_inside_and_out_and_on_their_walls() {
3887 let (count, length, outer, wall) = (6_u32, 0.1, 0.022, 0.0005);
3888 let m = model(&tube_finned_rocket(count, length, outer, wall));
3889 let a_ref = m.reference_area_m2();
3890 let n = f64::from(count);
3891 let inner = outer - wall;
3892 let (mach, reynolds_per_m) = (0.5, 4e6);
3893 let conditions = DragConditions::coasting(reynolds_per_m);
3894 let parts = m
3895 .buildup_components(&Flow::axial(mach), &conditions)
3896 .unwrap();
3897 let tubes = parts.iter().find(|c| c.id == "tube-fins").unwrap().drag;
3898 let terms = m.drag_terms.iter().find(|t| t.id == "tube-fins").unwrap();
3900 let cf = crate::drag::skin_friction_coefficient(
3901 reynolds_per_m * m.length_m,
3902 terms.relative_roughness,
3903 mach,
3904 )
3905 .unwrap();
3906 let wetted = n * 2.0 * PI * length * (outer + inner);
3907 close(tubes.friction, cf * wetted / a_ref, 1e-12, "friction");
3908 let annulus = n * PI * (outer * outer - inner * inner);
3910 let square = crate::drag::stagnation_drag_coefficient(mach).unwrap()
3911 + base_drag_coefficient(mach).unwrap();
3912 close(tubes.pressure, square * annulus / a_ref, 1e-12, "pressure");
3913 assert_eq!((tubes.base, tubes.parasitic), (0.0, 0.0));
3914 }
3915
3916 #[test]
3917 fn tube_fins_refuse_mach_0_8() {
3918 let m = model(&tube_finned_rocket(6, 0.1, 0.022, 0.0005));
3919 let refused = |e: AeroError| {
3922 matches!(e, AeroError::Mach { mach, limit, model }
3923 if mach == 0.8 && limit == 0.8 && model == "the tube-fin model")
3924 };
3925 let flow = Flow::new(0.8, 0.05, 0.0);
3926 let index = m.tube_fin_set_start();
3927 assert!(refused(m.normal_force(&flow).unwrap_err()));
3928 assert!(refused(m.components(&flow).unwrap_err()));
3929 assert!(refused(m.component_normal_force(index, &flow).unwrap_err()));
3930 assert!(refused(m.component_station_m(index, 0.8).unwrap_err()));
3931 assert!(refused(m.roll(0.8).unwrap_err()));
3932 let coasting = DragConditions::coasting(4e6);
3933 assert!(refused(m.drag(&Flow::axial(0.8), &coasting).unwrap_err()));
3934 assert!(refused(
3935 m.buildup_components(&Flow::axial(0.8), &coasting)
3936 .unwrap_err()
3937 ));
3938 let terms = m.drag_terms().iter().find(|t| t.id == "tube-fins").unwrap();
3940 let own = terms.evaluate(4e6, 0.8, &coasting, m.reference_area_m2());
3941 assert!(matches!(own, Err(AeroError::InComponent { id, source })
3942 if id == "tube-fins" && matches!(*source, AeroError::Mach { limit, .. } if limit == 0.8)));
3943 assert!(m.component_normal_force(0, &flow).is_ok());
3945 let below = Flow::new(0.8 - 1e-9, 0.05, 0.0);
3947 assert!(m.normal_force(&below).is_ok());
3948 assert!(m.drag(&below, &coasting).is_ok());
3949 }
3950
3951 #[test]
3954 fn unsupported_inputs_are_refused() {
3955 for (count, length_m, outer_radius_m, thickness_m, why) in [
3959 (2, 0.1, 0.01, 0.001, "2 tube fins"),
3960 (6, 0.1, 0.01, 0.01, "solid tube fins"),
3961 (6, 0.1, 0.023, 0.001, "tube fins that overlap"),
3962 (6, 0.005, 0.022, 0.001, "tube fins shorter than a third"),
3963 ] {
3964 let rocket = tube_finned_rocket(count, length_m, outer_radius_m, thickness_m);
3965 let err = AeroModel::new(&rocket.layout().unwrap()).unwrap_err();
3966 assert!(
3967 matches!(&err, AeroError::InComponent { id, source } if id == "tube-fins"
3968 && matches!(&**source, AeroError::Unsupported(what) if what.starts_with(why))),
3969 "{err}"
3970 );
3971 }
3972
3973 let mut rocket = podded_rocket(2);
3976 let pods = rocket.stages[0].components[1].children.last_mut().unwrap();
3977 pods.children[1].children.push(component(
3978 "pod-fins",
3979 fin_set(3, pod_fin_planform()),
3980 Some(Position::Bottom { aft_offset_m: 0.0 }),
3981 ));
3982 if let Part::FinSet(set) = &mut pods.children[1].children[0].part {
3983 set.cant_rad = 0.01;
3984 }
3985 let err = AeroModel::new(&rocket.layout().unwrap()).unwrap_err();
3986 assert!(
3987 matches!(&err, AeroError::InComponent { id, source } if id == "pod-fins"
3988 && matches!(&**source, AeroError::Unsupported(what) if what.starts_with("cant on a pod"))),
3989 "{err}"
3990 );
3991 let mut rocket = podded_rocket(2);
3992 let pods = rocket.stages[0].components[1].children.last_mut().unwrap();
3993 pods.children
3994 .push(component("pod-disc", body_part(0.0, 0.01, 0.01), None));
3995 let err = AeroModel::new(&rocket.layout().unwrap()).unwrap_err();
3996 assert!(
3997 matches!(&err, AeroError::InComponent { id, source } if id == "pod-disc"
3998 && matches!(&**source, AeroError::Unsupported(what) if what.contains("flat disc"))),
3999 "{err}"
4000 );
4001 let bare = AeroModel::new(&crate::testing::finned_rocket(4).layout().unwrap()).unwrap();
4003 let mut rocket = podded_rocket(2);
4004 let pods = rocket.stages[0].components[1].children.last_mut().unwrap();
4005 pods.children.clear();
4006 let empty = AeroModel::new(&rocket.layout().unwrap()).unwrap();
4007 assert_eq!(format!("{empty:?}"), format!("{bare:?}"));
4008
4009 let err = AeroModel::new(&finned_rocket(9).layout().unwrap()).unwrap_err();
4010 assert!(
4011 matches!(&err, AeroError::InComponent { id, .. } if id == "fins"),
4012 "{err}"
4013 );
4014
4015 let m = model(&finned_rocket(4));
4016 for bad in [
4017 flow(NORMAL_FORCE_MACH_LIMIT, 0.0, 0.0),
4018 flow(-0.01, 0.0, 0.0),
4019 flow(f64::NAN, 0.0, 0.0),
4020 ] {
4021 assert!(matches!(
4022 m.normal_force(&bad),
4023 Err(AeroError::Mach { limit, .. }) if limit == NORMAL_FORCE_MACH_LIMIT
4024 ));
4025 }
4026 assert!(m.normal_force(&flow(1.0, 0.1, 0.0)).is_ok());
4028 let coasting = DragConditions::coasting(1e7);
4029 assert!(m.drag(&flow(1.0, 0.0, 0.0), &coasting).is_ok());
4030 assert!(matches!(
4031 m.drag(&flow(5.0, 0.0, 0.0), &coasting),
4032 Err(AeroError::Mach { limit, model, .. })
4033 if limit == BUILDUP_MACH_LIMIT && model == "the drag buildup"
4034 ));
4035 assert!(
4036 m.buildup_components(&flow(5.0, 0.0, 0.0), &coasting)
4037 .is_err()
4038 );
4039 for bad in [
4040 flow(0.3, -1e-9, 0.0),
4041 flow(0.3, PI + 1e-9, 0.0),
4042 flow(0.3, f64::NAN, 0.0),
4043 flow(0.3, 0.1, f64::INFINITY),
4044 ] {
4045 assert!(matches!(
4046 m.normal_force(&bad),
4047 Err(AeroError::Domain { .. })
4048 ));
4049 assert!(m.components(&bad).is_err());
4050 }
4051
4052 let mut lugged = finned_rocket(4);
4053 lugged.stages[0].components[1].children.push(component(
4054 "lug",
4055 Part::LaunchLug(LaunchLug {
4056 length_m: 0.05,
4057 outer_radius_m: 0.004,
4058 thickness_m: 0.0005,
4059 angle_rad: 0.0,
4060 count: 1,
4061 spacing_m: 0.0,
4062 material: material(),
4063 }),
4064 Some(Position::Middle { aft_offset_m: 0.0 }),
4065 ));
4066 let with = model(&lugged).normal_force(&flow(0.3, 0.1, 0.0)).unwrap();
4067 assert_eq!(with, m.normal_force(&flow(0.3, 0.1, 0.0)).unwrap());
4068 }
4069
4070 #[test]
4073 fn a_bare_tube_has_no_cp_at_zero_incidence() {
4074 let tube = one_stage(
4075 vec![component("tube", body_part(1.0, 0.05, 0.05), None)],
4076 ReferenceDiameter::Maximum {},
4077 );
4078 let m = model(&tube);
4079 assert_eq!(
4080 m.normal_force(&Flow::axial(0.5)).unwrap().cp_station_m,
4081 None
4082 );
4083 let f = m.normal_force(&flow(0.5, 0.1, 0.0)).unwrap();
4084 close(f.cp_station_m.unwrap(), 0.5, 1e-15, "tube lift CP");
4085 }
4086
4087 proptest! {
4088 #[test]
4092 fn components_sum_to_the_total(
4093 count in 1u32..=8,
4094 mach in 0.0f64..0.99,
4095 alpha in 0.0f64..PI,
4096 roll in -4.0f64..4.0,
4097 ) {
4098 let m = model(&finned_rocket(count));
4099 let f = flow(mach, alpha, roll);
4100 let total = m.normal_force(&f).unwrap();
4101 let parts = m.components(&f).unwrap();
4102 let (c, moment) = parts.iter().fold((0.0, 0.0), |(c, x), p| {
4103 (c + p.normal_force.coefficient, x + p.normal_force.moment_m)
4104 });
4105 let tol = 1e-12 * (1.0 + total.coefficient.abs());
4106 prop_assert!((c - total.coefficient).abs() <= tol);
4107 prop_assert!((moment - total.moment_m).abs() <= 1e-12 * (1.0 + total.moment_m.abs()));
4108 let side: f64 = parts.iter().map(|p| p.normal_force.side_coefficient).sum();
4109 prop_assert!((side - total.side_coefficient).abs() <= 1e-12 * (1.0 + total.side_coefficient.abs()));
4110 let slope: f64 = parts.iter().map(|p| p.normal_force.slope_per_rad).sum();
4111 prop_assert!((slope - total.slope_per_rad).abs() <= 1e-12 * total.slope_per_rad.abs());
4112 prop_assert_eq!(m.component_count(), parts.len());
4115 let small = flow(mach, 1e-6, roll);
4116 for (index, part) in parts.iter().enumerate() {
4117 prop_assert_eq!(&m.component_normal_force(index, &f).unwrap(), &part.normal_force);
4118 let station = m.component_station_m(index, mach).unwrap();
4119 if let Some(cp) = m.component_normal_force(index, &small).unwrap().cp_station_m {
4120 prop_assert!((station - cp).abs() <= 1e-6 * (1.0 + cp.abs()));
4121 }
4122 }
4123 prop_assert!(m.component_normal_force(parts.len(), &f).is_err());
4124 prop_assert!(m.component_station_m(parts.len(), mach).is_err());
4125 }
4126
4127 #[test]
4128 fn scaling_leaves_slopes_and_scales_the_cp(
4129 k in 0.1f64..10.0,
4130 nose_fineness in 1.5f64..8.0,
4131 radius in 0.01f64..0.1,
4132 root in 0.02f64..0.3,
4133 tip_ratio in 0.0f64..1.0,
4134 span in 0.01f64..0.3,
4135 sweep in -0.1f64..0.3,
4136 count in 1u32..=8,
4137 alpha in 0.0f64..0.5,
4138 reference in 0.5f64..2.0,
4139 ) {
4140 let build = |s: f64, custom: Option<f64>| {
4141 let (r, l) = (s * radius, s * (root + 0.5));
4142 let mut tube = component("tube", body_part(l, r, r), None);
4143 let planform = FinPlanform::Trapezoidal {
4144 root_chord_m: s * root,
4145 tip_chord_m: s * root * tip_ratio,
4146 span_m: s * span,
4147 sweep_m: s * sweep,
4148 };
4149 tube.children = vec![component(
4150 "fins",
4151 fin_set(count, planform),
4152 Some(Position::Bottom { aft_offset_m: 0.0 }),
4153 )];
4154 let ogive = NoseShape::Ogive { radius_ratio: 1.0 };
4155 let nose = component("nose", nose(ogive, 2.0 * nose_fineness * r, r), None);
4156 let mut rocket = one_stage(vec![nose, tube], ReferenceDiameter::Maximum {});
4157 if let Some(d) = custom {
4158 rocket.reference_diameter = ReferenceDiameter::Custom { diameter_m: d };
4159 }
4160 model(&rocket).normal_force(&flow(0.4, alpha, 0.3)).unwrap()
4161 };
4162 let rel = |a: f64, b: f64| (a / b - 1.0).abs();
4163 let base = build(1.0, None);
4164 let base_cp = base.cp_station_m.unwrap();
4165 let scaled = build(k, None);
4166 prop_assert!(rel(scaled.slope_per_rad, base.slope_per_rad) < 1e-9);
4167 prop_assert!(rel(scaled.cp_station_m.unwrap(), k * base_cp) < 1e-9);
4168 let d = 2.0 * radius * reference;
4169 let custom = build(1.0, Some(d));
4170 let factor = (2.0 * radius / d).powi(2);
4171 prop_assert!(rel(custom.slope_per_rad, factor * base.slope_per_rad) < 1e-12);
4172 prop_assert!(rel(custom.cp_station_m.unwrap(), base_cp) < 1e-12);
4173 }
4174 }
4175
4176 #[test]
4180 fn radius_steps_count_at_the_joint() {
4181 let rocket = one_stage(
4182 vec![
4183 component("nose", nose(NoseShape::Conical {}, 0.2, 0.027), None),
4184 component("tube", body_part(0.5, 0.029, 0.029), None),
4185 component("tail", body_part(0.3, 0.025, 0.025), None),
4186 ],
4187 ReferenceDiameter::Maximum {},
4188 );
4189 let m = model(&rocket);
4190 let a_ref = PI * 0.029 * 0.029;
4191 let total = m.normal_force(&Flow::axial(0.3)).unwrap();
4192 close(
4193 total.slope_per_rad,
4194 2.0 * PI * 0.025 * 0.025 / a_ref,
4195 1e-14,
4196 "eq. 10",
4197 );
4198 let parts = m.components(&Flow::axial(0.3)).unwrap();
4199 let tube = parts[1].normal_force;
4200 close(
4201 tube.slope_per_rad,
4202 2.0 * PI * (0.029f64.powi(2) - 0.027f64.powi(2)) / a_ref,
4203 1e-14,
4204 "step up",
4205 );
4206 close(
4207 tube.cp_station_m.unwrap(),
4208 0.2,
4209 1e-14,
4210 "step up at the joint",
4211 );
4212 let tail = parts[2].normal_force;
4213 assert!(tail.slope_per_rad < 0.0);
4214 close(
4215 tail.cp_station_m.unwrap(),
4216 0.7,
4217 1e-14,
4218 "step down at the joint",
4219 );
4220 assert_eq!(m.bodies()[0].step_area_m2, 0.0);
4221 }
4222
4223 #[test]
4225 fn freeform_fins_through_the_model() {
4226 let trapezoid = finned_rocket(3);
4227 let mut freeform = trapezoid.clone();
4228 if let Part::FinSet(set) = &mut freeform.stages[0].components[3].children[0].part {
4229 set.planform = FinPlanform::Freeform {
4230 points_m: vec![[0.0, 0.0], [0.07, 0.06], [0.12, 0.06], [0.12, 0.0]],
4231 root_m: Vec::new(),
4232 };
4233 }
4234 let (a, b) = (model(&trapezoid), model(&freeform));
4235 let f = flow(0.7, 0.1, 0.0);
4236 let (fa, fb) = (a.normal_force(&f).unwrap(), b.normal_force(&f).unwrap());
4237 close(fb.coefficient, fa.coefficient, 1e-13, "C_N");
4238 close(
4239 fb.cp_station_m.unwrap(),
4240 fa.cp_station_m.unwrap(),
4241 1e-13,
4242 "CP",
4243 );
4244 }
4245
4246 #[test]
4248 fn flow_and_results_round_trip() {
4249 let f = flow(0.3, 0.1, -0.2);
4250 let back: Flow = serde_json::from_str(&serde_json::to_string(&f).unwrap()).unwrap();
4251 assert_eq!(back, f);
4252 assert!(
4253 serde_json::from_str::<Flow>(
4254 r#"{"mach":0.3,"alpha_rad":0.1,"roll_rad":0,"aoa_deg":5}"#
4255 )
4256 .is_err()
4257 );
4258 let n = model(&finned_rocket(4)).normal_force(&f).unwrap();
4259 let back: NormalForce = serde_json::from_str(&serde_json::to_string(&n).unwrap()).unwrap();
4260 assert_eq!(back, n);
4261 }
4262
4263 #[test]
4265 fn inconsistent_layouts_are_refused() {
4266 let layout = finned_rocket(4).layout().unwrap();
4267 let (fins, _) = layout.find("fins").unwrap();
4268 let mut no_radius = layout.clone();
4269 no_radius.components[fins].body_radius_m = None;
4270 let err = AeroModel::new(&no_radius).unwrap_err();
4271 assert!(
4272 matches!(&err, AeroError::InComponent { id, source } if id == "fins"
4273 && matches!(**source, AeroError::Layout(_))),
4274 "{err}"
4275 );
4276 assert_eq!(err, err.clone());
4277 let mut nan = layout;
4278 nan.components[0].fore_station_m = f64::NAN;
4279 assert!(matches!(
4280 AeroModel::new(&nan),
4281 Err(AeroError::InComponent { .. })
4282 ));
4283 }
4284
4285 #[test]
4288 fn two_fin_sets_push_across_the_flow() {
4289 let two = model(&finned_rocket(2));
4290 let set = &two.fin_sets()[0];
4291 let alpha = 0.05;
4292 let f = two.normal_force(&flow(0.4, alpha, FRAC_PI_4)).unwrap();
4293 let one_fin = set
4294 .fin
4295 .geometry()
4296 .single_fin_slope(two.reference_area_m2(), 0.4)
4297 .unwrap()
4298 * set.interference;
4299 close(f.side_coefficient, one_fin * alpha, 1e-13, "side");
4300 close(
4301 f.side_moment_m,
4302 one_fin * alpha * set.cp_station_m(0.4).unwrap(),
4303 1e-13,
4304 "side moment",
4305 );
4306 let fins = &two.components(&flow(0.4, alpha, FRAC_PI_4)).unwrap()[4];
4307 close(
4308 fins.normal_force.coefficient,
4309 one_fin * alpha,
4310 1e-13,
4311 "in plane",
4312 );
4313 let four = model(&finned_rocket(4))
4314 .normal_force(&flow(0.4, alpha, 0.3))
4315 .unwrap();
4316 assert_eq!((four.side_coefficient, four.side_moment_m), (0.0, 0.0));
4317 }
4318
4319 #[test]
4322 fn a_cancelled_slope_has_no_cp() {
4323 let rocket = one_stage(
4324 vec![
4325 component("nose", nose(NoseShape::Conical {}, 0.2, 0.0254), None),
4326 component("tube", body_part(0.5, 0.0254, 0.0254), None),
4327 component("tail", body_part(0.3, 0.0254, 1e-9), None),
4328 ],
4329 ReferenceDiameter::Maximum {},
4330 );
4331 let m = model(&rocket);
4332 let f = m.normal_force(&Flow::axial(0.3)).unwrap();
4333 assert!(f.slope_per_rad.abs() < 1e-12, "{}", f.slope_per_rad);
4334 assert_eq!(f.cp_station_m, None);
4335 let moving = m.normal_force(&flow(0.3, 0.01, 0.0)).unwrap();
4336 assert!(moving.moment_m.is_finite() && moving.cp_station_m.is_some());
4337 }
4338
4339 #[test]
4342 fn a_normal_force_table_replaces_the_sum() {
4343 use crate::table::{NormalForceColumn, NormalForceTable};
4344 use hpr_core::interp::{Extrapolation, Interpolation, Table1D};
4345
4346 let m = model(&finned_rocket(4));
4347 let flat = |value| {
4348 Table1D::new(
4349 vec![0.0, 2.0],
4350 vec![value, value],
4351 Interpolation::Linear,
4352 Extrapolation::Clamp,
4353 )
4354 .unwrap()
4355 };
4356 let table = NormalForceTable::new(vec![NormalForceColumn::new(0.0, flat(10.0), flat(0.9))])
4357 .unwrap();
4358 let with = m.clone().with_normal_force_table(table.clone()).unwrap();
4359 let at = flow(0.5, 0.02, 0.3);
4360 let replaced = with.normal_force(&at).unwrap();
4361 close(replaced.coefficient, 10.0 * 0.02_f64.sin(), 1e-15, "C_N");
4362 close(
4363 replaced.moment_m,
4364 replaced.coefficient * 0.9,
4365 1e-15,
4366 "moment",
4367 );
4368 assert_eq!(replaced.cp_station_m, Some(0.9));
4369 assert_eq!(
4370 (replaced.side_coefficient, replaced.side_moment_m),
4371 (0.0, 0.0)
4372 );
4373 assert!(replaced.table.is_some_and(|lookup| lookup.beyond_alpha));
4375 assert_eq!(m.normal_force(&at).unwrap().table, None);
4376 assert_eq!(with.components(&at).unwrap(), m.components(&at).unwrap());
4377 assert_eq!(
4378 with.component_normal_force(3, &at).unwrap(),
4379 m.component_normal_force(3, &at).unwrap()
4380 );
4381 assert_ne!(m.normal_force(&at).unwrap(), replaced);
4382 let wider = m
4384 .clone()
4385 .with_normal_force_table(table.clone().with_reference_diameter_m(0.108).unwrap())
4386 .unwrap();
4387 let scaled = wider.normal_force(&at).unwrap();
4388 close(
4389 scaled.coefficient,
4390 4.0 * replaced.coefficient,
4391 1e-14,
4392 "rescaled",
4393 );
4394 assert_eq!(scaled.cp_station_m, Some(0.9));
4395 let largest = m
4397 .clone()
4398 .with_normal_force_table(
4399 table
4400 .clone()
4401 .with_reference(TableReference::LargestBody)
4402 .unwrap(),
4403 )
4404 .unwrap();
4405 close(
4406 largest.normal_force(&at).unwrap().coefficient,
4407 replaced.coefficient,
4408 1e-14,
4409 "largest body",
4410 );
4411 let mut half = finned_rocket(4);
4414 half.reference_diameter = ReferenceDiameter::Custom { diameter_m: 0.027 };
4415 let half = model(&half)
4416 .with_normal_force_table(
4417 table
4418 .clone()
4419 .with_reference(TableReference::LargestBody)
4420 .unwrap(),
4421 )
4422 .unwrap();
4423 close(
4424 half.normal_force(&at).unwrap().coefficient,
4425 4.0 * replaced.coefficient,
4426 1e-14,
4427 "largest body on a half-size reference",
4428 );
4429 let hypersonic = |cp_at_6: f64| {
4431 let cps = Table1D::new(
4432 vec![0.0, 5.0, 6.0, 25.0],
4433 vec![0.9, 0.9, cp_at_6, -3.0],
4434 Interpolation::Linear,
4435 Extrapolation::Clamp,
4436 )
4437 .unwrap();
4438 NormalForceTable::new(vec![NormalForceColumn::new(0.0, flat(10.0), cps)]).unwrap()
4439 };
4440 assert!(m.clone().with_normal_force_table(hypersonic(0.8)).is_ok());
4441 assert!(m.clone().with_normal_force_table(hypersonic(-0.1)).is_err());
4442 for cp in [-0.01, 1.4] {
4444 let outside =
4445 NormalForceTable::new(vec![NormalForceColumn::new(0.0, flat(10.0), flat(cp))])
4446 .unwrap();
4447 assert!(matches!(
4448 m.clone().with_normal_force_table(outside),
4449 Err(AeroError::Domain { .. })
4450 ));
4451 }
4452 assert!(with.normal_force(&flow(6.0, 0.02, 0.0)).is_ok());
4454 assert!(matches!(
4455 m.normal_force(&flow(6.0, 0.02, 0.0)),
4456 Err(AeroError::Mach { .. })
4457 ));
4458 }
4459
4460 fn straight_rocket() -> hpr_design::Rocket {
4463 let mut rocket = crate::testing::finned_rocket(4);
4464 rocket.stages[0].components[2].part = body_part(0.05, 0.027, 0.027);
4465 rocket.stages[0].components[3].part = body_part(0.3, 0.027, 0.027);
4466 rocket
4467 }
4468
4469 fn straight_rocket_body() -> ShockExpansionBody {
4471 let nose =
4472 hpr_design::Profile::nose(NoseShape::Ogive { radius_ratio: 1.0 }, 0.25, 0.027).unwrap();
4473 let cylinder = |length_m| BodySegment::Cylinder {
4474 length_m,
4475 radius_m: 0.027,
4476 };
4477 ShockExpansionBody::new(
4478 &[
4479 BodySegment::Profile { profile: nose },
4480 cylinder(0.7),
4481 cylinder(0.05),
4482 cylinder(0.3),
4483 ],
4484 DEFAULT_ELEMENTS_PER_CURVE,
4485 )
4486 .unwrap()
4487 }
4488
4489 fn body_values(model: &AeroModel, mach: f64) -> Vec<[f64; 3]> {
4491 (0..model.bodies().len())
4492 .map(|index| {
4493 let force = model
4494 .component_normal_force(index, &flow(mach, 0.0, 0.0))
4495 .unwrap();
4496 let moment = force
4498 .cp_station_m
4499 .map_or(0.0, |cp| cp * force.slope_per_rad);
4500 [
4501 force.slope_per_rad,
4502 moment,
4503 model.component_station_m(index, mach).unwrap(),
4504 ]
4505 })
4506 .collect()
4507 }
4508
4509 fn no_jump(model: &AeroModel, mach: f64) {
4511 let below = body_values(model, mach - 1e-9);
4512 let above = body_values(model, mach + 1e-9);
4513 for (b, a) in below.iter().zip(&above) {
4514 for k in 0..3 {
4515 let scale = b[k].abs().max(a[k].abs()).max(1.0);
4516 assert!(
4517 (a[k] - b[k]).abs() <= 1e-7 * scale,
4518 "a jump at Mach {mach}: {b:?} to {a:?}"
4519 );
4520 }
4521 }
4522 }
4523
4524 #[test]
4525 fn the_supersonic_join_has_no_jump() {
4526 let model = model(&straight_rocket());
4527 let join = model.supersonic_body().unwrap();
4528 assert_eq!(join.covered, 4);
4529 assert_eq!(join.join_start_mach, SUPERSONIC_JOIN_START_MACH);
4530 let start = join.join_start_mach;
4532 for mach in [
4533 start,
4534 start + SUPERSONIC_JOIN_WIDTH_MACH,
4535 1.35,
4536 2.0,
4537 2.05,
4538 3.0,
4539 4.63,
4540 4.95,
4541 4.999,
4542 ] {
4543 no_jump(&model, mach);
4544 }
4545 let (low, high) = (body_values(&model, 1.0), body_values(&model, 2.0));
4547 assert_eq!(low[1][0], 0.0);
4548 assert!(high[1][0] > 0.1, "{:?}", high[1]);
4549 }
4550
4551 #[test]
4556 fn crossflow_and_the_boattail_fly_without_a_jump() {
4557 use crate::crossflow::{CROSSFLOW_DRAG_MACHS, ETA_MACHS};
4558 let layout = finned_rocket(4).layout().unwrap();
4559 for boattail in [
4560 SupersonicBoattail::WashingtonPettis,
4561 SupersonicBoattail::Footnote8,
4562 ] {
4563 let model = AeroModel::with_body_model(
4564 &layout,
4565 BodyModel::CURRENT.with_supersonic_boattail(boattail),
4566 )
4567 .unwrap();
4568 let start = model.supersonic_body().unwrap().join_start_mach;
4569 for alpha_deg in [10.0_f64, 30.0] {
4570 let s = alpha_deg.to_radians().sin();
4571 let mut machs = vec![start, start + SUPERSONIC_JOIN_WIDTH_MACH, 2.0, 2.05, 4.999];
4572 machs.extend(
4573 CROSSFLOW_DRAG_MACHS
4574 .iter()
4575 .chain(&ETA_MACHS)
4576 .map(|m| m / s)
4577 .filter(|&m| m > 1e-3 && m < 4.999),
4578 );
4579 for mach in machs {
4580 let at = |m: f64| {
4581 let f = model
4582 .normal_force(&flow(m, alpha_deg.to_radians(), 0.0))
4583 .unwrap();
4584 [f.coefficient, f.cp_station_m.unwrap()]
4585 };
4586 let (below, above) = (at(mach - 1e-9), at(mach + 1e-9));
4587 for k in 0..2 {
4588 let scale = below[k].abs().max(1.0);
4589 assert!(
4590 (above[k] - below[k]).abs() <= 1e-7 * scale,
4591 "{boattail:?} at {alpha_deg}° and Mach {mach}: {below:?} to {above:?}"
4592 );
4593 }
4594 }
4595 }
4596 }
4597 }
4598
4599 #[test]
4605 fn vertical_tips_fly_the_method_without_a_jump() {
4606 for shape in [
4607 NoseShape::PowerSeries { exponent: 0.6369 },
4608 NoseShape::PowerSeries { exponent: 0.5 },
4609 NoseShape::VON_KARMAN,
4610 NoseShape::LV_HAACK,
4611 NoseShape::Elliptical {},
4612 ] {
4613 for boattail in [true, false] {
4614 let mut rocket = if boattail {
4615 finned_rocket(4)
4616 } else {
4617 straight_rocket()
4618 };
4619 rocket.stages[0].components[0].part = nose(shape, 0.25, 0.027);
4620 let model = model(&rocket);
4621 let table = model
4622 .supersonic_body()
4623 .unwrap_or_else(|| panic!("{shape:?}: no table"));
4624 assert_eq!(table.covered, 4, "{shape:?}");
4625 if !boattail {
4627 let Part::NoseCone(cone) = &rocket.stages[0].components[0].part else {
4628 unreachable!("`nose` builds a nose cone")
4629 };
4630 let cylinder = |length_m| BodySegment::Cylinder {
4631 length_m,
4632 radius_m: 0.027,
4633 };
4634 let body = ShockExpansionBody::new(
4635 &[
4636 BodySegment::Profile {
4637 profile: cone.profile().unwrap(),
4638 },
4639 cylinder(0.7),
4640 cylinder(0.05),
4641 cylinder(0.3),
4642 ],
4643 DEFAULT_ELEMENTS_PER_CURVE,
4644 )
4645 .unwrap();
4646 let method = body.segment_slopes(3.0, model.reference_area_m2()).unwrap();
4647 let (nose_slope, _) = table.share(0, 3.0).unwrap();
4648 assert!(
4649 (nose_slope - method[0].slope_per_rad).abs() <= 1e-12,
4650 "{shape:?}: {nose_slope} against {:?}",
4651 method[0]
4652 );
4653 }
4654 let start = table.join_start_mach;
4655 let mut machs = vec![start, start + SUPERSONIC_JOIN_WIDTH_MACH, 2.1, 4.999];
4656 machs.extend(
4660 (SUPERSONIC_FIRST_STEP..SUPERSONIC_LAST_STEP).flat_map(|step| {
4661 let row = step as f64 / SUPERSONIC_STEPS_PER_MACH;
4662 [row, row + 0.013, row + 0.027, row + 0.041]
4663 }),
4664 );
4665 for alpha_deg in [1.0_f64, 10.0] {
4666 for &mach in &machs {
4667 let at = |m: f64| {
4668 let f = model
4669 .normal_force(&flow(m, alpha_deg.to_radians(), 0.0))
4670 .unwrap();
4671 [f.coefficient, f.cp_station_m.unwrap()]
4672 };
4673 let (below, above) = (at(mach - 1e-9), at(mach + 1e-9));
4674 for k in 0..2 {
4675 let scale = below[k].abs().max(1.0);
4676 assert!(
4677 (above[k] - below[k]).abs() <= 1e-7 * scale,
4678 "{shape:?} at {alpha_deg}° and Mach {mach}: {below:?} to {above:?}"
4679 );
4680 }
4681 }
4682 }
4683 }
4684 }
4685 }
4686
4687 #[test]
4693 fn a_blunt_tips_join_starts_where_its_cap_first_ends_on_the_nose() {
4694 let mut rocket = crate::testing::committed_design("wind-tunnel-arcas-robin-short.json");
4695 rocket.stages[0].components.truncate(3);
4697 let model = model(&rocket);
4698 let table = model.supersonic_body().unwrap();
4699 let Part::NoseCone(cone) = &rocket.stages[0].components[0].part else {
4700 unreachable!("the committed design starts with its nose")
4701 };
4702 let profile = cone.profile().unwrap();
4703 let base_angle = profile.radius_and_slope(profile.length_m()).1.atan();
4704 let (mut low, mut high) = (1.0 + 1e-9, 2.0);
4705 for _ in 0..200 {
4706 let mid = 0.5 * (low + high);
4707 if crate::blunt_tip::handover_angle_rad(mid).unwrap() < base_angle {
4708 low = mid;
4709 } else {
4710 high = mid;
4711 }
4712 }
4713 assert!(
4714 (table.join_start_mach - high).abs() <= 1e-12,
4715 "{} against {high}",
4716 table.join_start_mach
4717 );
4718 let run = model.supersonic_run.as_ref().unwrap();
4719 let body = ShockExpansionBody::new(&run.segments, DEFAULT_ELEMENTS_PER_CURVE).unwrap();
4720 let in_its_place: Vec<Option<ShockExpansionBody>> = run
4721 .boattails
4722 .iter()
4723 .map(|b| {
4724 b.as_ref().map(|b| {
4725 ShockExpansionBody::new(&b.in_its_place, DEFAULT_ELEMENTS_PER_CURVE).unwrap()
4726 })
4727 })
4728 .collect();
4729 let holds = |m: f64| {
4730 run.shares(&body, &in_its_place, None, m, model.reference_area_m2())
4731 .is_some()
4732 };
4733 for i in 1..=2000 {
4734 let d = f64::from(i) * 1e-10;
4735 assert!(
4736 !holds(table.join_start_mach - d),
4737 "holds {d} below the start"
4738 );
4739 assert!(
4740 holds(table.join_start_mach + d),
4741 "fails {d} above the start"
4742 );
4743 }
4744 }
4745
4746 #[test]
4747 fn a_long_lip_is_reported_as_the_runs_fallback() {
4748 let rocket = crate::testing::committed_design("wind-tunnel-arcas-robin-short.json");
4751 let components = &rocket.stages[0].components;
4752 let drop_m = components
4753 .iter()
4754 .find_map(|c| match &c.part {
4755 Part::Transition(t) if t.aft_radius_m < t.fore_radius_m => {
4756 Some(2.0 * (t.fore_radius_m - t.aft_radius_m))
4757 }
4758 _ => None,
4759 })
4760 .unwrap();
4761 let lip_index = components.iter().position(|c| c.id == "lip").unwrap();
4762 let with_lip = |length_m: f64| {
4763 let mut rocket = rocket.clone();
4764 let lip = &mut rocket.stages[0].components[lip_index];
4765 match &mut lip.part {
4766 Part::Transition(t) => t.length_m = length_m,
4767 Part::BodyTube(t) => t.length_m = length_m,
4768 other => panic!("the lip is a {}", other.kind_name()),
4769 }
4770 model(&rocket).supersonic_fallback()
4771 };
4772 assert_eq!(with_lip(drop_m * (1.0 - 1e-9)), None);
4773 assert_eq!(
4774 with_lip(drop_m * (1.0 + 1e-9)),
4775 Some(SupersonicFallback::LongLip {
4776 component: "lip".to_owned()
4777 })
4778 );
4779 }
4780
4781 #[test]
4786 fn a_lip_in_a_boattails_wake_carries_nothing() {
4787 for name in [
4788 "wind-tunnel-arcas-robin-short.json",
4789 "wind-tunnel-arcas-robin-long.json",
4790 ] {
4791 let rocket = crate::testing::committed_design(name);
4792 let model = model(&rocket);
4793 let table = model
4794 .supersonic_body()
4795 .unwrap_or_else(|| panic!("{name}: no table"));
4796 assert_eq!(table.covered, 4, "{name}");
4798 let lip = model.bodies().last().unwrap();
4799 assert!(lip.slope_per_rad > 0.1, "{name}: the lip is a flare");
4800 for mach in [1.3, 2.0, 3.0, 5.0] {
4801 let (slope, moment) = table.share(3, mach).unwrap();
4802 assert_eq!((slope, moment), (0.0, 0.0), "{name} at Mach {mach}");
4803 }
4804 let start = table.join_start_mach;
4806 let lip_slope = |mach: f64| {
4808 model
4809 .components(&Flow::axial(mach))
4810 .unwrap()
4811 .into_iter()
4812 .find(|c| c.id == "lip")
4813 .unwrap()
4814 .normal_force
4815 .slope_per_rad
4816 };
4817 let below = lip_slope(1.0);
4818 assert!(
4819 (below - lip.slope_per_rad).abs() <= 0.01 * lip.slope_per_rad,
4820 "{name}: {below} against slender-body theory's {}",
4821 lip.slope_per_rad
4822 );
4823 let above = lip_slope(start + SUPERSONIC_JOIN_WIDTH_MACH);
4824 assert!(above.abs() <= 1e-12, "{name}: {above} above the join");
4825 for alpha_deg in [1.0_f64, 10.0] {
4826 let mut machs = vec![start, start + SUPERSONIC_JOIN_WIDTH_MACH, 4.999];
4827 machs.extend(
4828 (SUPERSONIC_FIRST_STEP..SUPERSONIC_LAST_STEP).flat_map(|step| {
4829 let row = step as f64 / SUPERSONIC_STEPS_PER_MACH;
4830 [row, row + 0.017, row + 0.033]
4831 }),
4832 );
4833 for mach in machs {
4834 let at = |m: f64| {
4835 let f = model
4836 .normal_force(&flow(m, alpha_deg.to_radians(), 0.0))
4837 .unwrap();
4838 [f.coefficient, f.cp_station_m.unwrap()]
4839 };
4840 let (below, above) = (at(mach - 1e-9), at(mach + 1e-9));
4841 for k in 0..2 {
4842 let scale = below[k].abs().max(1.0);
4843 assert!(
4844 (above[k] - below[k]).abs() <= 1e-7 * scale,
4845 "{name} at {alpha_deg}° and Mach {mach}: {below:?} to {above:?}"
4846 );
4847 }
4848 }
4849 }
4850 }
4851 let lipped = |rise: f64| {
4856 let mut rocket = crate::testing::finned_rocket(4);
4857 rocket.stages[0].components.truncate(3);
4861 rocket.stages[0].components.push(component(
4862 "lip",
4863 body_part(0.01, 0.022, 0.022 + 0.5 * rise * 0.010),
4864 None,
4865 ));
4866 model(&rocket)
4867 };
4868 assert_eq!(lipped(0.2).supersonic_body().map(|t| t.covered), Some(4));
4869 assert!(lipped(0.6).supersonic_body().is_none());
4870 let at = |rise: f64| {
4871 let model = lipped(rise);
4872 let f = model
4873 .normal_force(&flow(3.0, 4f64.to_radians(), 0.0))
4874 .unwrap();
4875 (f.coefficient, f.cp_station_m.unwrap())
4876 };
4877 let (below, above) = (at(0.2499), at(0.2501));
4882 assert!(
4883 (above.0 - below.0).abs() <= 2e-4 * below.0.abs()
4884 && (above.1 - below.1).abs() <= 2e-4 * below.1.abs(),
4885 "{below:?} to {above:?} across the wake's full-shelter rise"
4886 );
4887 let out_of_wake = {
4891 let mut rocket = crate::testing::finned_rocket(4);
4892 rocket.stages[0].components.truncate(3);
4893 rocket.stages[0].components.push(component(
4894 "lip",
4895 body_part(0.0101, 0.022, 0.022 + 0.5 * 0.25 * 0.010),
4896 None,
4897 ));
4898 let m = model(&rocket);
4899 assert!(
4900 m.supersonic_body().is_none(),
4901 "out of the wake by its length"
4902 );
4903 let f = m.normal_force(&flow(3.0, 4f64.to_radians(), 0.0)).unwrap();
4904 (
4905 f.coefficient,
4906 f.cp_station_m.unwrap() / m.reference_diameter_m(),
4907 )
4908 };
4909 let in_wake = at(0.25);
4910 let diameter_m = model(&crate::testing::finned_rocket(4)).reference_diameter_m();
4911 let gap_force = out_of_wake.0 / in_wake.0 - 1.0;
4912 let gap_calibers = out_of_wake.1 - in_wake.1 / diameter_m;
4913 assert!(
4914 (gap_force + 0.3295).abs() < 5e-4 && (gap_calibers + 1.774).abs() < 5e-3,
4915 "at one shape the two models differ by {gap_force} in force and {gap_calibers} \
4916 calibres in center of pressure"
4917 );
4918 let (full, nearly_none) = (at(0.25), at(0.4999));
4921 let calibers = (nearly_none.1 - full.1) / diameter_m;
4922 assert!(
4923 (nearly_none.0 / full.0 - 1.0 + 0.291).abs() < 5e-4 && (calibers + 0.932).abs() < 5e-3,
4924 "across the band: {full:?} to {nearly_none:?}, {calibers} calibres"
4925 );
4926 assert!((lipped(0.25).supersonic_body().unwrap().shape_weight - 1.0).abs() < 1e-12);
4929 for (rise, want) in [(0.3_f64, 0.8_f64), (0.375, 0.5), (0.45, 0.2)] {
4932 let weight = lipped(rise).supersonic_body().unwrap().shape_weight;
4933 assert!(
4934 (weight - want).abs() < 1e-9,
4935 "at a rise of {rise} the wake covers {weight}, not {want}"
4936 );
4937 }
4938 assert!(
4941 lipped(0.49999).supersonic_body().unwrap().shape_weight < 1e-4,
4942 "a hair inside the wake's far edge the method has almost no weight left"
4943 );
4944 assert!(
4945 lipped(0.5).supersonic_body().is_none(),
4946 "at the wake's far edge there is no run"
4947 );
4948 let gapped = |gap_m: f64| {
4951 let mut rocket = crate::testing::finned_rocket(4);
4952 rocket.stages[0].components.truncate(3);
4953 rocket.stages[0].components.push(component(
4954 "gap",
4955 body_part(gap_m, 0.022, 0.022),
4956 None,
4957 ));
4958 rocket.stages[0].components.push(component(
4959 "lip",
4960 body_part(0.01, 0.022, 0.022 + 0.5 * 0.17 * 0.010),
4961 None,
4962 ));
4963 let m = model(&rocket);
4964 let weight = m.supersonic_body().map_or(0.0, |t| t.shape_weight);
4965 let f = m.normal_force(&flow(3.0, 4f64.to_radians(), 0.0)).unwrap();
4966 (
4967 weight,
4968 f.coefficient,
4969 f.cp_station_m.unwrap() / m.reference_diameter_m(),
4970 )
4971 };
4972 let (near, far) = (gapped(1e-6), gapped(0.010));
4973 assert!(
4974 near.0 > 0.999 && far.0 == 0.0,
4975 "a tube of the boattail's own drop in diameter carries the lip out of the wake: \
4976 {near:?} to {far:?}"
4977 );
4978 assert!(
4979 (far.2 - near.2 + 1.973).abs() < 0.01 && (far.1 / near.1 - 1.0 + 0.3381).abs() < 5e-4,
4980 "over that tube the force moves {} and the center of pressure {} calibres",
4981 far.1 / near.1 - 1.0,
4982 far.2 - near.2
4983 );
4984 let half = lipped(0.375);
4989 let tube = half
4990 .component_normal_force(1, &flow(3.0, 0.0, 0.0))
4991 .unwrap();
4992 let station = half.component_station_m(1, 3.0).unwrap();
4993 assert!(
4994 (station - tube.cp_station_m.unwrap() - 0.1142).abs() < 5e-4,
4995 "the tube's station {station} against its center of pressure {:?}",
4996 tube.cp_station_m
4997 );
4998 assert!(
4999 lipped(0.55).supersonic_body().is_none(),
5000 "past the wake there is no run at all"
5001 );
5002 for rise in [0.2499_f64, 0.25, 0.3, 0.375, 0.45, 0.4999] {
5004 let (low, high) = (at(rise - 1e-9), at(rise + 1e-9));
5005 assert!(
5006 (high.0 - low.0).abs() <= 1e-7 * low.0.abs().max(1.0)
5007 && (high.1 - low.1).abs() <= 1e-7 * low.1.abs().max(1.0),
5008 "at a rise of {rise}: {low:?} to {high:?}"
5009 );
5010 }
5011 for shape in [NoseShape::Conical {}, NoseShape::Haack { parameter: 0.5 }] {
5015 let mut rocket = crate::testing::committed_design("wind-tunnel-arcas-robin-short.json");
5016 let components = &mut rocket.stages[0].components;
5017 let last = components.len() - 1;
5018 let Part::Transition(lip) = &mut components[last].part else {
5019 unreachable!("the committed design ends in its lip")
5020 };
5021 lip.aft_radius_m = lip.fore_radius_m + 0.60 * (0.028575 - 0.0166116);
5023 lip.shape = shape;
5024 assert!(
5025 model(&rocket).supersonic_body().is_none(),
5026 "{shape:?} out of the wake"
5027 );
5028 }
5029 let mut rocket = finned_rocket(4);
5032 rocket.stages[0].components.truncate(3);
5033 rocket.stages[0].components.push(component(
5034 "long flare",
5035 body_part(1.0, 0.022, 0.0229),
5036 None,
5037 ));
5038 assert!(model(&rocket).supersonic_body().is_none());
5039 let mut rocket = finned_rocket(4);
5042 rocket.stages[0].components.truncate(3);
5043 rocket.stages[0].components.push(component(
5044 "second boattail",
5045 body_part(0.03, 0.0231, 0.021),
5046 None,
5047 ));
5048 let narrowing = model(&rocket);
5049 assert!(narrowing.bodies().last().unwrap().slope_per_rad < 0.0);
5050 assert!(narrowing.supersonic_body().is_none());
5051 }
5052
5053 #[test]
5059 fn a_separating_boattail_reads_the_correlation_at_its_steepest_measured_angle() {
5060 let mach = 1.5;
5065 let ogive = nose(NoseShape::Ogive { radius_ratio: 1.0 }, 0.25, 0.027);
5066 let Part::NoseCone(ogive) = ogive else {
5067 unreachable!("`nose` builds a nose cone")
5068 };
5069 let kept_pair = |at: f64, half_angle_deg: f64| {
5072 let length_m = (0.027 - 0.022) / half_angle_deg.to_radians().tan();
5073 let mut rocket = finned_rocket(4);
5074 rocket.stages[0].components[2].part = body_part(length_m, 0.027, 0.022);
5075 let model = model(&rocket);
5076 let a_ref = model.reference_area_m2();
5077 let table = model
5078 .supersonic_body()
5079 .unwrap_or_else(|| panic!("a table at {half_angle_deg}°"));
5080 let share = table.share(2, at).expect("the boattail's share").0;
5081 let in_its_place = ShockExpansionBody::new(
5082 &[
5083 BodySegment::Profile {
5084 profile: ogive.profile().unwrap(),
5085 },
5086 BodySegment::Cylinder {
5087 length_m: 0.7,
5088 radius_m: 0.027,
5089 },
5090 BodySegment::Cylinder {
5091 length_m,
5092 radius_m: 0.027,
5093 },
5094 ],
5095 DEFAULT_ELEMENTS_PER_CURVE,
5096 )
5097 .unwrap();
5098 let cylinder = in_its_place.segment_slopes(at, a_ref).unwrap()[2].slope_per_rad;
5099 let raw = |l: f64| {
5100 crate::supersonic_boattail::wp_slope(at, 0.027, 0.022, l).unwrap()
5101 * PI
5102 * 0.027
5103 * 0.027
5104 / a_ref
5105 };
5106 (share - cylinder, raw(length_m))
5108 };
5109 let kept = |half_angle_deg: f64| kept_pair(mach, half_angle_deg);
5110 let kept_at = |at: f64, half_angle_deg: f64| kept_pair(at, half_angle_deg).0;
5111 let at_16 = (0.027 - 0.022) / 16.0_f64.to_radians().tan();
5112 let close = |got: f64, want: f64, what: &str| {
5113 assert!(
5114 (got - want).abs() < 0.02 * want.abs(),
5115 "{what}: {got} against {want}"
5116 );
5117 };
5118 for angle in [8.0, 15.0, 16.0] {
5120 let (flown, raw) = kept(angle);
5121 close(flown, raw, "read at the true angle");
5122 }
5123 let a_ref = model(&finned_rocket(4)).reference_area_m2();
5129 let at_16_read = crate::supersonic_boattail::wp_slope(mach, 0.027, 0.022, at_16).unwrap()
5130 * PI
5131 * 0.027
5132 * 0.027
5133 / a_ref;
5134 for angle in [17.0, 23.0, 30.0, 40.0] {
5135 let (flown, raw) = kept(angle);
5136 let (held, _) = kept(16.0);
5137 close(flown, held, "held at 16°");
5138 close(flown, at_16_read, "the correlation at the 16° geometry");
5140 assert!(flown < 0.0, "{angle}°: the boattail still takes lift off");
5141 assert!(flown < raw, "{angle}°: {flown} against {raw} read raw");
5144 }
5145 let length_m = (0.027 - 0.022) / 30.0_f64.to_radians().tan();
5150 let mut steep = finned_rocket(4);
5151 steep.stages[0].components[2].part = body_part(length_m, 0.027, 0.022);
5152 steep.stages[0].components.truncate(3);
5153 let steep = model(&steep);
5154 let body_cp_calibers = |at: f64, increment_kept: bool| {
5155 let parts = steep.components(&Flow::axial(at)).unwrap();
5156 let (mut slope, mut moment) = (0.0, 0.0);
5157 for (index, part) in parts.iter().take(steep.bodies().len()).enumerate() {
5158 let mut share = part.normal_force.slope_per_rad;
5159 let station = part.normal_force.cp_station_m.unwrap_or(0.0);
5160 if index == 2 && !increment_kept {
5161 share -= kept_at(at, 30.0);
5162 }
5163 slope += share;
5164 moment += share * station;
5165 }
5166 moment / slope / steep.reference_diameter_m()
5167 };
5168 let gaps: Vec<(f64, f64)> = [1.5_f64, 2.0, 3.0, 4.63]
5170 .iter()
5171 .map(|at| {
5172 (
5173 *at,
5174 body_cp_calibers(*at, false) - body_cp_calibers(*at, true),
5175 )
5176 })
5177 .collect();
5178 for (at, gap) in &gaps {
5179 assert!(
5180 *gap > 0.6,
5181 "Mach {at}: the fading rule sits {gap} calibres aft"
5182 );
5183 }
5184 for (at, want) in [(1.5, 1.35), (2.0, 0.91), (3.0, 0.75), (4.63, 0.67)] {
5186 let got = gaps
5187 .iter()
5188 .find(|(m, _)| (m - at).abs() < 1e-9)
5189 .expect("a measured gap")
5190 .1;
5191 assert!(
5192 (got - want).abs() < 0.02,
5193 "Mach {at}: the gap is {got} calibres, the guide says {want}"
5194 );
5195 }
5196 for angle in [15.9_f64, 16.0, 16.1] {
5201 let (below, above) = (kept(angle - 1e-6).0, kept(angle + 1e-6).0);
5202 assert!((above - below).abs() < 1e-6, "{angle}°: {below} to {above}");
5203 }
5204 assert!((at_16 - (0.027 - 0.022) / 16.0_f64.to_radians().tan()).abs() < 1e-15);
5205 }
5206
5207 #[test]
5212 fn holding_the_correlation_stops_at_potential_flow() {
5213 let (fore_radius_m, aft_radius_m, mach) = (0.027_f64, 0.05 * 0.027_f64, 1.42_f64);
5218 let slender = 2.0 * ((aft_radius_m / fore_radius_m).powi(2) - 1.0);
5219 let length_m = (fore_radius_m - aft_radius_m) / 30.0_f64.to_radians().tan();
5220 let mut rocket = finned_rocket(4);
5221 rocket.stages[0].components[2].part = body_part(length_m, fore_radius_m, aft_radius_m);
5222 rocket.stages[0].components[3].part = body_part(0.3, aft_radius_m, aft_radius_m);
5223 let steep = model(&rocket);
5224 let a_ref = steep.reference_area_m2();
5225 let per_boattail_area = PI * fore_radius_m * fore_radius_m / a_ref;
5226 let table = steep.supersonic_body().expect("the method covers it");
5227 assert!(
5230 mach > table.join_start_mach,
5231 "Mach {mach} is below the table's start, {}",
5232 table.join_start_mach
5233 );
5234 let cylinder = ShockExpansionBody::new(
5235 &[
5236 BodySegment::Profile {
5237 profile: match nose(NoseShape::Ogive { radius_ratio: 1.0 }, 0.25, 0.027) {
5238 Part::NoseCone(ogive) => ogive.profile().unwrap(),
5239 _ => unreachable!("`nose` builds a nose cone"),
5240 },
5241 },
5242 BodySegment::Cylinder {
5243 length_m: 0.7,
5244 radius_m: fore_radius_m,
5245 },
5246 BodySegment::Cylinder {
5247 length_m,
5248 radius_m: fore_radius_m,
5249 },
5250 ],
5251 DEFAULT_ELEMENTS_PER_CURVE,
5252 )
5253 .unwrap()
5254 .segment_slopes(mach, a_ref)
5255 .unwrap()[2]
5256 .slope_per_rad;
5257 let flown = table.share(2, mach).expect("the boattail's share").0 - cylinder;
5258 assert!(
5260 (flown / per_boattail_area - slender).abs() < 2e-3,
5261 "{} against potential flow's {slender}",
5262 flown / per_boattail_area
5263 );
5264 let held = crate::supersonic_boattail::wp_slope(
5265 mach,
5266 fore_radius_m,
5267 aft_radius_m,
5268 (fore_radius_m - aft_radius_m) / SEPARATION_ONSET_RAD.tan(),
5269 )
5270 .unwrap();
5271 assert!(
5272 held < slender - 0.05,
5273 "the held read {held} must pass the bound's {slender} to pin it"
5274 );
5275 let gentle_m = (fore_radius_m - 0.6 * fore_radius_m) / 4.0_f64.to_radians().tan();
5278 let gentle =
5279 crate::supersonic_boattail::wp_slope(1.5, fore_radius_m, 0.6 * fore_radius_m, gentle_m)
5280 .unwrap();
5281 let gentle_slender = 2.0 * (0.6_f64.powi(2) - 1.0);
5282 assert!(
5283 gentle < gentle_slender,
5284 "the 4° read {gentle} should pass {gentle_slender}"
5285 );
5286 let mut gentle_rocket = finned_rocket(4);
5287 gentle_rocket.stages[0].components[2].part =
5288 body_part(gentle_m, fore_radius_m, 0.6 * fore_radius_m);
5289 gentle_rocket.stages[0].components[3].part =
5290 body_part(0.3, 0.6 * fore_radius_m, 0.6 * fore_radius_m);
5291 let gentle_model = model(&gentle_rocket);
5292 let gentle_table = gentle_model.supersonic_body().expect("a table");
5293 let gentle_cylinder = ShockExpansionBody::new(
5294 &[
5295 BodySegment::Profile {
5296 profile: match nose(NoseShape::Ogive { radius_ratio: 1.0 }, 0.25, 0.027) {
5297 Part::NoseCone(ogive) => ogive.profile().unwrap(),
5298 _ => unreachable!("`nose` builds a nose cone"),
5299 },
5300 },
5301 BodySegment::Cylinder {
5302 length_m: 0.7,
5303 radius_m: fore_radius_m,
5304 },
5305 BodySegment::Cylinder {
5306 length_m: gentle_m,
5307 radius_m: fore_radius_m,
5308 },
5309 ],
5310 DEFAULT_ELEMENTS_PER_CURVE,
5311 )
5312 .unwrap()
5313 .segment_slopes(1.5, gentle_model.reference_area_m2())
5314 .unwrap()[2]
5315 .slope_per_rad;
5316 let gentle_flown = gentle_table.share(2, 1.5).expect("the share").0 - gentle_cylinder;
5317 let gentle_area = PI * fore_radius_m * fore_radius_m / gentle_model.reference_area_m2();
5318 assert!(
5319 (gentle_flown / gentle_area - gentle).abs() < 2e-3,
5320 "the 4° boattail flies {} and its correlation reads {gentle}",
5321 gentle_flown / gentle_area
5322 );
5323 }
5324
5325 #[test]
5333 fn what_the_potential_flow_bound_reaches() {
5334 let fore_radius_m = 0.027_f64;
5335 let ogive = nose(NoseShape::Ogive { radius_ratio: 1.0 }, 0.25, 0.027);
5336 let Part::NoseCone(ogive) = ogive else {
5337 unreachable!("`nose` builds a nose cone")
5338 };
5339 let mut worst = (0.0_f64, 0.0_f64, 0.0_f64, 0.0_f64);
5340 let mut steepest_tabled = 0.0_f64;
5341 let mut floor_angles: Vec<f64> = Vec::new();
5342 for angle_deg in [16.0_f64, 16.5, 17.0, 17.25, 17.5, 30.0, 53.0, 53.5, 53.6] {
5343 for ratio in [0.001_f64, 0.02, 0.25, 0.3] {
5344 let aft_radius_m = ratio * fore_radius_m;
5345 let drop_m = fore_radius_m - aft_radius_m;
5346 let length_m = drop_m / angle_deg.to_radians().tan();
5347 let held_length_m = drop_m / SEPARATION_ONSET_RAD.tan();
5348 let ceiling = 2.0 * (ratio * ratio - 1.0);
5349 let read = |mach: f64, length: f64| {
5350 crate::supersonic_boattail::wp_slope(mach, fore_radius_m, aft_radius_m, length)
5351 .unwrap()
5352 };
5353 if read(1.2, held_length_m) >= ceiling {
5356 continue;
5357 }
5358 let mut rocket = finned_rocket(4);
5359 rocket.stages[0].components[2].part =
5360 body_part(length_m, fore_radius_m, aft_radius_m);
5361 rocket.stages[0].components[3].part = body_part(0.3, aft_radius_m, aft_radius_m);
5362 let flown = model(&rocket);
5363 let Some(table) = flown.supersonic_body() else {
5364 continue;
5365 };
5366 steepest_tabled = steepest_tabled.max(angle_deg);
5367 let a_ref = flown.reference_area_m2();
5368 let per_area = PI * fore_radius_m * fore_radius_m / a_ref;
5369 let cylinder_body = ShockExpansionBody::new(
5372 &[
5373 BodySegment::Profile {
5374 profile: ogive.profile().unwrap(),
5375 },
5376 BodySegment::Cylinder {
5377 length_m: 0.7,
5378 radius_m: fore_radius_m,
5379 },
5380 BodySegment::Cylinder {
5381 length_m,
5382 radius_m: fore_radius_m,
5383 },
5384 ],
5385 DEFAULT_ELEMENTS_PER_CURVE,
5386 )
5387 .unwrap();
5388 let mut rows = vec![table.join_start_mach];
5393 let mut step = (table.join_start_mach * SUPERSONIC_STEPS_PER_MACH).ceil();
5394 while step / SUPERSONIC_STEPS_PER_MACH <= 1.55 {
5395 rows.push(step / SUPERSONIC_STEPS_PER_MACH);
5396 step += 1.0;
5397 }
5398 for mach in rows {
5399 let Some((share, _)) = table.share(2, mach) else {
5400 break;
5401 };
5402 let cylinder =
5403 cylinder_body.segment_slopes(mach, a_ref).unwrap()[2].slope_per_rad;
5404 let increment = (share - cylinder) / per_area;
5405 if read(mach, length_m) >= ceiling {
5406 assert!(
5410 increment >= ceiling - 1e-12,
5411 "{angle_deg}° to {ratio} of the radius at Mach {mach}: the boattail \
5412 takes {increment} off, past potential flow's {ceiling}"
5413 );
5414 } else {
5415 assert!(
5418 increment <= ceiling,
5419 "{angle_deg}° to {ratio} at Mach {mach}: {increment} against {ceiling}"
5420 );
5421 let published = read(mach, length_m);
5424 assert!(
5425 (increment - published).abs() < 1e-9,
5426 "{angle_deg}° to {ratio} at Mach {mach}: {increment} against the \
5427 correlation's own {published}"
5428 );
5429 if !floor_angles.contains(&angle_deg) {
5430 floor_angles.push(angle_deg);
5431 }
5432 }
5433 let moved = (increment - read(mach, held_length_m)).abs()
5435 * per_area
5436 * table.weight(mach);
5437 if moved > worst.0 {
5438 worst = (moved, angle_deg, ratio, mach);
5439 }
5440 }
5441 }
5442 }
5443 assert_eq!(
5447 floor_angles,
5448 [16.0, 16.5, 17.0, 17.25],
5449 "the swept angles where a boattail's own read already passes potential flow"
5450 );
5451 assert!(
5454 (steepest_tabled - 53.5).abs() < 1e-12,
5455 "the steepest boattail swept that the method tables is {steepest_tabled}°"
5456 );
5457 assert!(
5458 (worst.0 - 0.060).abs() < 5e-4,
5459 "the bound moves a printed coefficient by at most {:.4} per rad, at {}° to {} of the \
5460 radius at Mach {:.3}",
5461 worst.0,
5462 worst.1,
5463 worst.2,
5464 worst.3
5465 );
5466 }
5467
5468 #[test]
5480 fn footnote_eights_boattail_share_by_hand() {
5481 let a_ref = PI * 0.027 * 0.027;
5482 let ogive = nose(NoseShape::Ogive { radius_ratio: 1.0 }, 0.25, 0.027);
5483 let Part::NoseCone(ogive) = ogive else {
5484 unreachable!("`nose` builds a nose cone")
5485 };
5486 let (length_m, fore_radius_m, aft_radius_m) = (0.05, 0.027, 0.022);
5487 const TUBE_LENGTH_M: f64 = 0.2;
5488 let body = ShockExpansionBody::new(
5489 &[
5490 BodySegment::Profile {
5491 profile: ogive.profile().unwrap(),
5492 },
5493 BodySegment::Cylinder {
5494 length_m: 0.7,
5495 radius_m: 0.027,
5496 },
5497 BodySegment::Profile {
5498 profile: hpr_design::Profile::transition(
5499 NoseShape::Conical {},
5500 length_m,
5501 fore_radius_m,
5502 aft_radius_m,
5503 false,
5504 )
5505 .unwrap(),
5506 },
5507 BodySegment::Cylinder {
5508 length_m: TUBE_LENGTH_M,
5509 radius_m: aft_radius_m,
5510 },
5511 ],
5512 DEFAULT_ELEMENTS_PER_CURVE,
5513 )
5514 .unwrap();
5515 let mach = 2.0;
5516 let shares = body.segment_slopes(mach, a_ref).unwrap();
5517 let flows = body.element_flows(mach).unwrap();
5520 let boattail = flows
5521 .iter()
5522 .rev()
5523 .nth(1)
5524 .copied()
5525 .expect("the boattail's element");
5526 let (load, decay) = (boattail.loading_per_rad, boattail.decay_per_m);
5527 assert!((boattail.corner_radius_m - fore_radius_m).abs() < 1e-12);
5529 let delta = ((aft_radius_m - fore_radius_m) / length_m).atan();
5532 let cone_load = delta.tan() * 2.0;
5533 assert!((boattail.tangent_cone_loading_per_rad - cone_load).abs() < 1e-12);
5534 assert!((boattail.tangent_cone_pressure_ratio - 1.0).abs() < 1e-12);
5535 let integrate = |span_m: f64, fore_r: f64, aft_r: f64, load: f64, cone: f64, decay: f64| {
5537 let steps = 4000;
5538 let mut sum = 0.0;
5539 for i in 0..=steps {
5540 let t = f64::from(i) / f64::from(steps);
5541 let x = t * span_m;
5542 let r = fore_r + (aft_r - fore_r) * t;
5543 let e = (-decay * x).exp();
5544 let lambda = (1.0 - e) * cone + e * load;
5545 let weight = if i == 0 || i == steps {
5546 1.0
5547 } else if i % 2 == 1 {
5548 4.0
5549 } else {
5550 2.0
5551 };
5552 sum += weight * lambda * r;
5553 }
5554 2.0 * PI * (sum * span_m / (3.0 * f64::from(steps))) / a_ref
5555 };
5556 let by_hand = integrate(
5557 length_m,
5558 fore_radius_m,
5559 aft_radius_m,
5560 load,
5561 cone_load,
5562 decay,
5563 );
5564 assert!(
5565 (shares[2].slope_per_rad - by_hand).abs() < 1e-6,
5566 "{} against {by_hand}",
5567 shares[2].slope_per_rad
5568 );
5569 let tube = flows.last().expect("the tube's element");
5572 assert!(tube.angle_rad.abs() < 1e-12 && tube.tangent_cone_loading_per_rad.abs() < 1e-12);
5573 let tube_by_hand = integrate(
5574 TUBE_LENGTH_M,
5575 aft_radius_m,
5576 aft_radius_m,
5577 tube.loading_per_rad,
5578 0.0,
5579 tube.decay_per_m,
5580 );
5581 assert!(
5582 (shares[3].slope_per_rad - tube_by_hand).abs() < 1e-6,
5583 "the tube: {} against {tube_by_hand}",
5584 shares[3].slope_per_rad
5585 );
5586 assert!(
5588 tube_by_hand < by_hand && tube_by_hand < 0.0,
5589 "the tube's {tube_by_hand} against the boattail's {by_hand}"
5590 );
5591 let slender = 2.0 * (aft_radius_m.powi(2) - fore_radius_m.powi(2)) / (0.027 * 0.027);
5594 assert!(
5595 by_hand < 0.0 && by_hand / slender < 0.2,
5596 "{by_hand} against slender-body theory's {slender}"
5597 );
5598 }
5599
5600 #[test]
5604 fn the_boattail_takes_washington_and_pettis_increment() {
5605 let layout = finned_rocket(4).layout().unwrap();
5606 let current = AeroModel::new(&layout).unwrap();
5607 let before = AeroModel::with_body_model(&layout, BodyModel::BEFORE_M1_8E6).unwrap();
5608 let a_ref = current.reference_area_m2();
5609 let ogive = nose(NoseShape::Ogive { radius_ratio: 1.0 }, 0.25, 0.027);
5610 let Part::NoseCone(ogive) = ogive else {
5611 unreachable!("`nose` builds a nose cone")
5612 };
5613 let segments = |last: BodySegment| {
5614 [
5615 BodySegment::Profile {
5616 profile: ogive.profile().unwrap(),
5617 },
5618 BodySegment::Cylinder {
5619 length_m: 0.7,
5620 radius_m: 0.027,
5621 },
5622 last,
5623 ]
5624 };
5625 let in_its_place = ShockExpansionBody::new(
5626 &segments(BodySegment::Cylinder {
5627 length_m: 0.05,
5628 radius_m: 0.027,
5629 }),
5630 DEFAULT_ELEMENTS_PER_CURVE,
5631 )
5632 .unwrap();
5633 let rocket = finned_rocket(4);
5634 let Part::Transition(tail) = &rocket.stages[0].components[2].part else {
5635 unreachable!("the test rocket's third part is its boattail")
5636 };
5637 let with_boattail = ShockExpansionBody::new(
5638 &segments(BodySegment::Profile {
5639 profile: tail.profile().unwrap(),
5640 }),
5641 DEFAULT_ELEMENTS_PER_CURVE,
5642 )
5643 .unwrap();
5644 for mach in [2.0, 3.0, 4.5] {
5645 let cylinder = in_its_place.segment_slopes(mach, a_ref).unwrap()[2];
5646 let increment = crate::supersonic_boattail::wp_slope(mach, 0.027, 0.022, 0.05).unwrap()
5647 * PI
5648 * 0.027
5649 * 0.027
5650 / a_ref;
5651 let center = 0.25 + 0.7 + wp_center_fraction(mach) * 0.05;
5652 let (slope, moment) = current.supersonic_body().unwrap().share(2, mach).unwrap();
5653 close(
5654 slope,
5655 cylinder.slope_per_rad + increment,
5656 1e-12,
5657 "W&P slope",
5658 );
5659 close(
5660 moment,
5661 cylinder.moment_slope_m + increment * center,
5662 1e-12,
5663 "W&P moment",
5664 );
5665 assert!(increment < 0.0 && cylinder.slope_per_rad.abs() < 0.1 * increment.abs());
5667 let footnote_8 = with_boattail.segment_slopes(mach, a_ref).unwrap()[2];
5668 let (old, _) = before.supersonic_body().unwrap().share(2, mach).unwrap();
5669 close(old, footnote_8.slope_per_rad, 1e-12, "footnote 8");
5670 let slender = -2.0 * (1.0 - (0.022_f64 / 0.027).powi(2));
5673 assert!(slender < slope && slope < old, "Mach {mach}: {slope} {old}");
5674 }
5675 }
5676
5677 #[test]
5682 fn two_boattails_each_take_their_own_increment() {
5683 let (l_n, l_1, l_b1, l_2, l_b2) = (0.25, 0.4, 0.05, 0.3, 0.04);
5684 let (r_0, r_1, r_2) = (0.027, 0.024, 0.02);
5685 let rocket = one_stage(
5686 vec![
5687 component(
5688 "nose",
5689 nose(NoseShape::Ogive { radius_ratio: 1.0 }, l_n, r_0),
5690 None,
5691 ),
5692 component("tube", body_part(l_1, r_0, r_0), None),
5693 component("boattail", body_part(l_b1, r_0, r_1), None),
5694 component("waist", body_part(l_2, r_1, r_1), None),
5695 component("tail", body_part(l_b2, r_1, r_2), None),
5696 ],
5697 ReferenceDiameter::Maximum {},
5698 );
5699 let model = AeroModel::new(&rocket.layout().unwrap()).unwrap();
5700 let before =
5701 AeroModel::with_body_model(&rocket.layout().unwrap(), BodyModel::BEFORE_M1_8E6)
5702 .unwrap();
5703 let table = model.supersonic_body().unwrap();
5704 assert_eq!(table.covered, 5);
5705 let a_ref = model.reference_area_m2();
5706 let profile = |part: &Part| match part {
5707 Part::NoseCone(n) => n.profile().unwrap(),
5708 Part::Transition(t) => t.profile().unwrap(),
5709 _ => unreachable!("the test's profiled parts are a nose and transitions"),
5710 };
5711 let parts: Vec<Part> = rocket.stages[0]
5712 .components
5713 .iter()
5714 .map(|c| c.part.clone())
5715 .collect();
5716 let cylinder = |length_m, radius_m| BodySegment::Cylinder { length_m, radius_m };
5717 let profiled = |i: usize| BodySegment::Profile {
5718 profile: profile(&parts[i]),
5719 };
5720 let body = |segments: &[BodySegment]| {
5721 ShockExpansionBody::new(segments, DEFAULT_ELEMENTS_PER_CURVE).unwrap()
5722 };
5723 let whole = body(&[
5724 profiled(0),
5725 cylinder(l_1, r_0),
5726 profiled(2),
5727 cylinder(l_2, r_1),
5728 profiled(4),
5729 ]);
5730 let first_in_place = body(&[profiled(0), cylinder(l_1, r_0), cylinder(l_b1, r_0)]);
5731 let second_in_place = body(&[
5732 profiled(0),
5733 cylinder(l_1, r_0),
5734 profiled(2),
5735 cylinder(l_2, r_1),
5736 cylinder(l_b2, r_1),
5737 ]);
5738 let increment = |mach: f64, fore: f64, aft: f64, length: f64| {
5739 crate::supersonic_boattail::wp_slope(mach, fore, aft, length).unwrap()
5740 * PI
5741 * fore
5742 * fore
5743 / a_ref
5744 };
5745 for mach in [2.0, 3.5] {
5746 let method = whole.segment_slopes(mach, a_ref).unwrap();
5747 let share = |i| table.share(i, mach).unwrap().0;
5748 for i in [0, 1, 3] {
5749 close(
5750 share(i),
5751 method[i].slope_per_rad,
5752 1e-12,
5753 "the method's share",
5754 );
5755 }
5756 let first = first_in_place.segment_slopes(mach, a_ref).unwrap()[2].slope_per_rad
5757 + increment(mach, r_0, r_1, l_b1);
5758 let second = second_in_place.segment_slopes(mach, a_ref).unwrap()[4].slope_per_rad
5759 + increment(mach, r_1, r_2, l_b2);
5760 close(share(2), first, 1e-12, "first boattail");
5761 close(share(4), second, 1e-12, "second boattail");
5762 let old = before.supersonic_body().unwrap();
5764 close(
5765 old.share(2, mach).unwrap().0,
5766 method[2].slope_per_rad,
5767 1e-12,
5768 "fn 8",
5769 );
5770 close(
5771 old.share(4, mach).unwrap().0,
5772 method[4].slope_per_rad,
5773 1e-12,
5774 "fn 8",
5775 );
5776 assert!(share(2) < method[2].slope_per_rad && share(4) < method[4].slope_per_rad);
5777 }
5778 }
5779
5780 #[test]
5781 fn a_boattailed_body_flies_the_method_without_a_jump() {
5782 let model = model(&finned_rocket(4));
5785 let join = model.supersonic_body().unwrap();
5786 assert_eq!(join.covered, 4);
5787 let start = join.join_start_mach;
5788 for mach in [
5789 start,
5790 start + SUPERSONIC_JOIN_WIDTH_MACH,
5791 1.35,
5792 2.0,
5793 2.05,
5794 3.0,
5795 4.63,
5796 4.95,
5797 4.999,
5798 ] {
5799 no_jump(&model, mach);
5800 }
5801 let (low, high) = (body_values(&model, 1.0), body_values(&model, 2.0));
5802 let (boattail, share) = (&model.bodies()[2], join.share(2, 2.0).unwrap());
5803 assert!(boattail.slope_per_rad < 0.0 && share.0 < 0.0, "{share:?}");
5804 assert!(
5806 (high[2][0] - share.0).abs() <= 1e-12 * share.0.abs(),
5807 "{share:?}"
5808 );
5809 assert_eq!(low[2][2], high[2][2]);
5812 assert!((0.95..=1.0).contains(&high[2][2]), "{:?}", high[2]);
5813 assert!(high[3][0] < 0.0, "{:?}", high[3]);
5814 for mach in [1.2, 1.3, 1.5, 3.0, 4.999] {
5815 let at = body_values(&model, mach);
5816 assert_eq!((at[2][2], at[3][2]), (low[2][2], low[3][2]), "Mach {mach}");
5817 }
5818 assert!((1.0..=1.3).contains(&low[3][2]), "{:?}", low[3]);
5819 assert_eq!(low[1][0], 0.0);
5821 assert!(high[1][0] > 0.1, "{:?}", high[1]);
5822 }
5823
5824 #[test]
5825 fn a_blunter_cone_joins_where_the_method_starts_to_hold() {
5826 let mut rocket = straight_rocket();
5829 let length = 0.027 / 20.0_f64.to_radians().tan();
5830 rocket.stages[0].components[0].part = nose(NoseShape::Conical {}, length, 0.027);
5831 let model = model(&rocket);
5832 let join = model.supersonic_body().unwrap();
5833 let start = join.join_start_mach;
5834 assert!(start > SUPERSONIC_JOIN_START_MACH, "{start}");
5835 for mach in [start, start + SUPERSONIC_JOIN_WIDTH_MACH, start + 0.5] {
5836 no_jump(&model, mach);
5837 }
5838 assert_eq!(body_values(&model, start), body_values(&model, 0.5));
5840 }
5841
5842 #[test]
5843 fn the_joins_start_moves_with_the_nose_not_in_steps() {
5844 let start_at = |degrees: f64| {
5847 let mut rocket = straight_rocket();
5848 let length = 0.027 / degrees.to_radians().tan();
5849 rocket.stages[0].components[0].part = nose(NoseShape::Conical {}, length, 0.027);
5850 let model = model(&rocket);
5851 let start = model.supersonic_body().unwrap().join_start_mach;
5852 (model, start)
5853 };
5854 let (model, start) = start_at(20.0);
5855 let grid = start * SUPERSONIC_STEPS_PER_MACH;
5856 assert!((grid - grid.round()).abs() > 1e-3, "on the grid: {start}");
5857 let first_row = grid.ceil() / SUPERSONIC_STEPS_PER_MACH;
5859 for mach in [start, first_row, start + SUPERSONIC_JOIN_WIDTH_MACH] {
5860 no_jump(&model, mach);
5861 }
5862 assert_eq!(body_values(&model, start), body_values(&model, 0.5));
5863 assert!((start - 1.341910).abs() < 1e-6, "{start}");
5865 let join = model.supersonic_body().unwrap();
5870 let (lead, row) = (
5871 join.share(1, start).unwrap(),
5872 join.share(1, first_row).unwrap(),
5873 );
5874 assert!(lead.0 < 1e-5 && row.0 > 0.3, "{lead:?} against {row:?}");
5875 let (_, nudged) = start_at(20.0 + 1e-6);
5877 assert!((nudged - start).abs() < 1e-7, "{start} to {nudged}");
5878 let starts: Vec<f64> = [20.0, 20.1, 20.2, 20.3, 20.4, 20.5]
5880 .into_iter()
5881 .map(|degrees| start_at(degrees).1)
5882 .collect();
5883 assert!((starts[5] - 1.355500).abs() < 1e-6, "{starts:?}");
5884 for pair in starts.windows(2) {
5885 assert!(pair[1] > pair[0] && pair[1] - pair[0] < 0.02, "{starts:?}");
5886 }
5887 }
5888
5889 #[test]
5890 fn below_the_join_the_bodies_keep_slender_body_terms() {
5891 let model = model(&straight_rocket());
5892 let slender: Vec<[f64; 3]> = model
5893 .bodies()
5894 .iter()
5895 .enumerate()
5896 .map(|(index, body)| {
5897 [
5898 body.slope_per_rad,
5899 body.moment_slope_m,
5900 model.component_station_m(index, 0.0).unwrap(),
5901 ]
5902 })
5903 .collect();
5904 for mach in [0.0, 0.5, 0.99, 1.1, SUPERSONIC_JOIN_START_MACH] {
5905 let got = body_values(&model, mach);
5906 for (index, (g, want)) in got.iter().zip(&slender).enumerate() {
5907 assert_eq!(g[0], want[0], "body {index} at Mach {mach}");
5908 assert_eq!(g[2], want[2], "body {index} at Mach {mach}");
5909 if want[0] != 0.0 {
5910 close(g[1], want[1], 1e-14, "moment");
5911 }
5912 }
5913 }
5914 }
5915
5916 #[test]
5917 fn past_the_join_the_covered_bodies_take_the_method() {
5918 let model = model(&straight_rocket());
5919 let area = model.reference_area_m2();
5920 let body = straight_rocket_body();
5921 let join_end = SUPERSONIC_JOIN_START_MACH + SUPERSONIC_JOIN_WIDTH_MACH;
5922 for (mach, rel) in [
5924 (1.5, 1e-12),
5925 (2.0, 1e-12),
5926 (3.0, 1e-12),
5927 (4.95, 1e-12),
5928 (2.96, 1e-3),
5929 ] {
5930 assert!(mach >= join_end);
5931 let want = body.slope(mach, area).unwrap();
5932 let values = body_values(&model, mach);
5933 let slope: f64 = values.iter().map(|v| v[0]).sum();
5934 let moment: f64 = values.iter().map(|v| v[1]).sum();
5935 close(slope, want.slope_per_rad, rel, "slope");
5936 close(
5937 moment / slope,
5938 want.center_of_pressure_m,
5939 rel,
5940 "center of pressure",
5941 );
5942 let bounds = [(0.0, 0.25), (0.25, 0.95), (0.95, 1.0), (1.0, 1.3)];
5944 for (v, (fore, aft)) in values.iter().zip(bounds) {
5945 assert!(v[2] > fore && v[2] < aft, "{v:?} not in {fore} to {aft}");
5946 }
5947 }
5948 }
5949
5950 fn flared_rocket(flare_deg: f64) -> hpr_design::Rocket {
5953 let (fore_r, flare_l) = (0.027, 0.3);
5954 let aft_r = fore_r + flare_l * flare_deg.to_radians().tan();
5955 let mut tail = component("tail", body_part(0.2, aft_r, aft_r), None);
5956 tail.children = vec![component(
5957 "fins",
5958 fin_set(
5959 4,
5960 FinPlanform::Trapezoidal {
5961 root_chord_m: 0.12,
5962 tip_chord_m: 0.05,
5963 span_m: 0.06,
5964 sweep_m: 0.07,
5965 },
5966 ),
5967 Some(Position::Bottom { aft_offset_m: 0.0 }),
5968 )];
5969 one_stage(
5970 vec![
5971 component(
5972 "nose",
5973 nose(NoseShape::Ogive { radius_ratio: 1.0 }, 0.25, fore_r),
5974 None,
5975 ),
5976 component("body", body_part(0.7, fore_r, fore_r), None),
5977 component("flare", body_part(flare_l, fore_r, aft_r), None),
5978 tail,
5979 ],
5980 ReferenceDiameter::Maximum {},
5981 )
5982 }
5983
5984 fn flared_at(rocket: &hpr_design::Rocket, mach: f64) -> (f64, f64) {
5987 let model = model(rocket);
5988 let force = model.normal_force(&flow(mach, 1e-4, 0.0)).unwrap();
5989 (force.coefficient / 1e-4, force.cp_station_m.unwrap())
5990 }
5991
5992 fn flared_run(rocket: &hpr_design::Rocket) -> ShockExpansionBody {
5995 let model = model(rocket);
5996 let run = model.supersonic_run.as_ref().expect("a run with a flare");
5997 let flare = run.flare.as_ref().expect("the flare in the run");
5998 ShockExpansionBody::new(&flare.ahead, DEFAULT_ELEMENTS_PER_CURVE).unwrap()
5999 }
6000
6001 fn corner_limit_deg(ahead: &ShockExpansionBody, mach: f64) -> f64 {
6005 let aft = ahead.aft_flow(mach).unwrap();
6006 (crate::shock_expansion::flare_corner_limit_rad(aft.surface_mach).unwrap() + aft.angle_rad)
6007 .min(crate::blunt_tip::CONE_TABLE_CAP_RAD)
6008 .to_degrees()
6009 }
6010
6011 #[test]
6016 fn a_flared_body_flies_the_method() {
6017 let rocket = flared_rocket(10.0);
6018 let model = model(&rocket);
6019 let body = model
6020 .supersonic_body()
6021 .expect("the method covers the flare");
6022 assert_eq!(body.covered, 3);
6025 let old = AeroModel::with_body_model(
6026 &rocket.layout().unwrap(),
6027 BodyModel::CURRENT.with_supersonic_flare(SupersonicFlare::SlenderBody),
6028 )
6029 .unwrap();
6030 assert!(old.supersonic_body().is_none());
6031 for mach in [2.0, 3.0, 4.95] {
6033 let (slope, moment) = body.share(2, mach).unwrap();
6034 assert!(slope > 0.0, "the flare's share at Mach {mach} is {slope}");
6035 let station = moment / slope;
6036 assert!(
6037 (0.95..=1.25).contains(&station),
6038 "the flare's share acts at {station} m at Mach {mach}"
6039 );
6040 }
6041 let diameter_m = model.reference_diameter_m();
6045 assert!((diameter_m - 0.159_796).abs() < 5e-7, "{diameter_m} m");
6046 for (mach, want_slope, want_calibers, want_old_slope, want_forward) in [
6047 (2.0, 3.537_087_3, 7.365_003, 3.715_259_2, 0.089_048),
6048 (3.0, 2.889_886_8, 7.014_772, 3.081_544_6, 0.168_593),
6049 (4.95, 2.488_028_5, 6.680_843, 2.643_091_4, 0.237_654),
6050 ] {
6051 let (slope, station_m) = flared_at(&rocket, mach);
6052 let old = old.normal_force(&flow(mach, 1e-4, 0.0)).unwrap();
6053 let forward = (old.cp_station_m.unwrap() - station_m) / diameter_m;
6054 assert!(
6055 (slope - want_slope).abs() < 5e-7
6056 && (station_m / diameter_m - want_calibers).abs() < 5e-6,
6057 "Mach {mach}: {slope} per rad at {} calibres",
6058 station_m / diameter_m
6059 );
6060 assert!(
6061 (old.coefficient / 1e-4 - want_old_slope).abs() < 5e-7
6062 && (forward - want_forward).abs() < 5e-6,
6063 "Mach {mach} on slender-body theory: {} per rad, the method {forward} calibres \
6064 forward",
6065 old.coefficient / 1e-4
6066 );
6067 }
6068 }
6069
6070 #[test]
6078 fn a_flare_is_read_no_steeper_than_its_corners_shock_holds() {
6079 let ahead = flared_run(&flared_rocket(18.5));
6080 for (mach, want_surface_mach, want_limit_deg) in [
6081 (1.5, 1.499_968_751_529, 12.111_850_220_062),
6082 (2.0, 1.999_780_928_628, 22.969_761_173_077),
6083 (2.5, 2.498_954_764_087, 29.786_310_648_004),
6084 (3.0, 2.996_526_706_276, 30.0),
6086 (4.95, 4.892_298_698_955, 30.0),
6087 ] {
6088 let aft = ahead.aft_flow(mach).unwrap();
6089 assert!(
6090 (aft.surface_mach - want_surface_mach).abs() < 5e-10 && aft.angle_rad == 0.0,
6091 "Mach {mach}: the corner turns Mach {} at {} rad",
6092 aft.surface_mach,
6093 aft.angle_rad
6094 );
6095 let limit = corner_limit_deg(&ahead, mach);
6096 assert!(
6097 (limit - want_limit_deg).abs() < 5e-10,
6098 "Mach {mach}: the corner is read to {limit}°"
6099 );
6100 }
6101 let steep = model(&flared_rocket(30.0));
6103 let run = steep.supersonic_run.as_ref().unwrap();
6104 let flare = run.flare.as_ref().unwrap();
6105 let rise_m = flare.aft_radius_m - flare.fore_radius_m;
6106 let drawn_length_m = rise_m / corner_limit_deg(&ahead, 2.0).to_radians().tan();
6107 assert!(drawn_length_m > flare.length_m);
6108 let mut segments = flare.ahead.clone();
6109 segments.push(BodySegment::Profile {
6110 profile: hpr_design::Profile::transition(
6111 NoseShape::Conical {},
6112 drawn_length_m,
6113 flare.fore_radius_m,
6114 flare.aft_radius_m,
6115 false,
6116 )
6117 .unwrap(),
6118 });
6119 let drawn = ShockExpansionBody::new(&segments, DEFAULT_ELEMENTS_PER_CURVE).unwrap();
6120 let area = steep.reference_area_m2();
6121 let want = *drawn.segment_slopes(2.0, area).unwrap().last().unwrap();
6122 let body = steep.supersonic_body().unwrap();
6123 let (slope, moment) = body.share(2, 2.0).unwrap();
6124 assert!(
6125 (slope - want.slope_per_rad).abs() < 1e-12 * want.slope_per_rad,
6126 "the held share is {slope}, the drawn-out flare's {}",
6127 want.slope_per_rad
6128 );
6129 let (fore_m, aft_m) = (run.fore_m[2], run.fore_m[2] + flare.length_m);
6133 let drawn_station_m = want.moment_slope_m / want.slope_per_rad;
6134 let station_m = moment / slope;
6135 assert!(
6136 station_m > fore_m && station_m < aft_m,
6137 "the held share acts at {station_m}, off the flare's {fore_m} to {aft_m}"
6138 );
6139 assert!(
6142 (drawn_station_m - 1.207_758).abs() < 5e-6 && station_m < drawn_station_m - 0.05,
6143 "the drawn-out share acts at {drawn_station_m}, the mapped one at {station_m}"
6144 );
6145 let at = |deg: f64| {
6150 let rocket = flared_rocket(deg);
6151 let new = model(&rocket).normal_force(&flow(2.0, 1e-4, 0.0)).unwrap();
6152 let old = AeroModel::with_body_model(
6153 &rocket.layout().unwrap(),
6154 BodyModel::CURRENT.with_supersonic_flare(SupersonicFlare::SlenderBody),
6155 )
6156 .unwrap()
6157 .normal_force(&flow(2.0, 1e-4, 0.0))
6158 .unwrap();
6159 (new.coefficient / 1e-4, old.coefficient / 1e-4)
6160 };
6161 for (deg, want_method, want_slender) in
6162 [(30.0, 1.946_513, 2.307_693), (75.0, 1.637_856, 2.010_350)]
6163 {
6164 let (method, slender) = at(deg);
6165 assert!(
6166 (method - want_method).abs() < 5e-6 && (slender - want_slender).abs() < 5e-6,
6167 "{deg}°: the method reads {method}, slender-body theory {slender}"
6168 );
6169 assert!(
6170 method < slender,
6171 "{deg}° should read below slender-body theory"
6172 );
6173 }
6174 let at_limit = model(&flared_rocket(corner_limit_deg(&ahead, 2.0)));
6177 let limit_run = at_limit.supersonic_run.as_ref().unwrap();
6178 let limit_body =
6179 ShockExpansionBody::new(&limit_run.segments, DEFAULT_ELEMENTS_PER_CURVE).unwrap();
6180 let marched = *limit_body
6181 .segment_slopes(2.0, at_limit.reference_area_m2())
6182 .unwrap()
6183 .last()
6184 .unwrap();
6185 let (limit_slope, limit_moment) =
6186 at_limit.supersonic_body().unwrap().share(2, 2.0).unwrap();
6187 assert!(
6188 (limit_slope - marched.slope_per_rad).abs() < 1e-12 * marched.slope_per_rad
6189 && (limit_moment - marched.moment_slope_m).abs() < 1e-12 * marched.moment_slope_m,
6190 "at the limit the reading is {limit_slope} at {}, the march's {} at {}",
6191 limit_moment / limit_slope,
6192 marched.slope_per_rad,
6193 marched.moment_slope_m / marched.slope_per_rad
6194 );
6195 }
6196
6197 #[test]
6206 fn nothing_jumps_where_the_flares_shock_detaches() {
6207 let ahead = flared_run(&flared_rocket(18.5));
6210 let boundary_deg = corner_limit_deg(&ahead, 2.0);
6211 assert!((boundary_deg - 22.969_761_173_077).abs() < 5e-12);
6212 for (epsilon, want) in [(1e-9, 4.527e-11), (1e-7, 4.527e-9), (1e-5, 4.527e-7)] {
6213 let below = flared_at(&flared_rocket(boundary_deg - epsilon), 2.0);
6214 let above = flared_at(&flared_rocket(boundary_deg + epsilon), 2.0);
6215 let gap = (above.0 / below.0 - 1.0).abs();
6216 assert!(
6217 (gap - want).abs() < 0.01 * want,
6218 "±{epsilon}° across the boundary moves the slope by {gap}, not {want}"
6219 );
6220 assert!(
6221 (above.1 - below.1).abs() < 2e-3 * epsilon,
6222 "the station moves"
6223 );
6224 }
6225 let rocket = flared_rocket(18.5);
6234 let model = model(&rocket);
6235 let run = model.supersonic_run.as_ref().expect("a run with a flare");
6236 let body = ShockExpansionBody::new(&run.segments, DEFAULT_ELEMENTS_PER_CURVE).unwrap();
6237 let in_its_place: Vec<Option<ShockExpansionBody>> =
6238 run.segments.iter().map(|_| None).collect();
6239 let area = model.reference_area_m2();
6240 let (mut low, mut high) = (1.2_f64, 5.0_f64);
6241 for _ in 0..80 {
6242 let mid = 0.5 * (low + high);
6243 if corner_limit_deg(&ahead, mid) < 18.5 {
6244 low = mid;
6245 } else {
6246 high = mid;
6247 }
6248 }
6249 assert!(
6250 (high - 1.767_666_917_849).abs() < 5e-12,
6251 "18.5° attaches from Mach {high}"
6252 );
6253 assert!(corner_limit_deg(&ahead, high - 1e-9) < 18.5);
6256 assert!(corner_limit_deg(&ahead, high + 1e-9) > 18.5);
6257 let flare_share = |mach: f64| {
6258 let shares = run
6259 .shares(&body, &in_its_place, Some(&ahead), mach, area)
6260 .expect("the row holds either side of the crossing");
6261 shares[2]
6262 };
6263 for (epsilon, want) in [(1e-9, 3.622e-10), (1e-7, 3.622e-8), (1e-5, 3.622e-6)] {
6264 let (below, above) = (flare_share(high - epsilon), flare_share(high + epsilon));
6265 let gap = (above.slope_per_rad / below.slope_per_rad - 1.0).abs();
6266 assert!(
6267 (gap - want).abs() < 0.01 * want,
6268 "±{epsilon} in Mach across the boundary moves the flare's share by {gap}, not \
6269 {want}"
6270 );
6271 let station = |s: SegmentSlope| s.moment_slope_m / s.slope_per_rad;
6273 assert!(
6274 (station(above) - station(below)).abs() < 0.1 * epsilon
6275 && (0.95..=1.25).contains(&station(below)),
6276 "the flare's share acts at {} then {}",
6277 station(below),
6278 station(above)
6279 );
6280 }
6281 for epsilon in [1e-9, 1e-7, 1e-5] {
6283 let below = flared_at(&rocket, high - epsilon);
6284 let above = flared_at(&rocket, high + epsilon);
6285 assert!(
6286 (above.0 / below.0 - 1.0).abs() < 2.0 * epsilon,
6287 "±{epsilon} in Mach moves the rocket's slope by {}",
6288 above.0 / below.0 - 1.0
6289 );
6290 }
6291 }
6292
6293 #[test]
6314 fn a_near_flat_flare_marches_every_row_and_the_fallback_is_still_measured() {
6315 let ahead = flared_run(&flared_rocket(1.0));
6318 let at_mach_3 =
6319 crate::shock_expansion::flare_reduction_turns_rad(&ahead.aft_flow(3.0).unwrap())
6320 .unwrap();
6321 let reduced_at_mach_3 =
6322 at_mach_3.crossing_rad.to_degrees()..at_mach_3.balance_rad.to_degrees();
6323 assert!(
6324 reduced_at_mach_3.contains(&0.007) && reduced_at_mach_3.contains(&0.008),
6325 "Mach 3 reduces {reduced_at_mach_3:?}, which was meant to hold 0.007° and 0.008°"
6326 );
6327 for deg in [0.001, 0.007, 0.008, 0.01, 0.03, 0.038_17, 0.045, 0.0589] {
6328 let model = model(&flared_rocket(deg));
6329 let run = model.supersonic_run.as_ref().expect("a run with a flare");
6330 let body = ShockExpansionBody::new(&run.segments, DEFAULT_ELEMENTS_PER_CURVE).unwrap();
6331 let area = model.reference_area_m2();
6332 let refused = (SUPERSONIC_FIRST_STEP..=SUPERSONIC_LAST_STEP).find(|step| {
6333 body.slope(*step as f64 / SUPERSONIC_STEPS_PER_MACH, area)
6334 .is_err()
6335 });
6336 assert_eq!(
6337 refused,
6338 None,
6339 "a {deg}° flare loses Mach {:?}",
6340 refused.map(|step| step as f64 / SUPERSONIC_STEPS_PER_MACH)
6341 );
6342 let flows = body.element_flows(3.0).unwrap();
6345 assert_eq!(
6346 flows[flows.len() - 1].decay_per_m == 0.0,
6347 reduced_at_mach_3.contains(°),
6348 "a {deg}° flare at Mach 3, where the region is {reduced_at_mach_3:?}"
6349 );
6350 }
6351 let dropped = |deg: f64, mach: f64| {
6354 let rocket = flared_rocket(deg);
6355 let layout = rocket.layout().unwrap();
6356 let read = |model: AeroModel| {
6357 let force = model
6358 .normal_force(&flow(mach, 4f64.to_radians(), 0.0))
6359 .unwrap();
6360 (
6361 force.coefficient,
6362 force.cp_station_m.unwrap() / model.reference_diameter_m(),
6363 )
6364 };
6365 let method = read(AeroModel::new(&layout).unwrap());
6366 let bare = read(
6367 AeroModel::with_body_model(
6368 &layout,
6369 BodyModel::CURRENT.with_supersonic_flare(SupersonicFlare::SlenderBody),
6370 )
6371 .unwrap(),
6372 );
6373 (bare.0 / method.0 - 1.0, bare.1 - method.1)
6374 };
6375 let (force, calibers) = dropped(0.058_820_517_4, 3.0);
6376 assert!(
6377 (force + 0.083_0).abs() < 5e-5 && (calibers - 1.157_4).abs() < 5e-5,
6378 "at the band's steep edge the fallback moves the force {force:.4} and the center of \
6379 pressure {calibers:.4} calibres"
6380 );
6381 let (force, calibers) = dropped(9.018_246_55e-4, 2.0);
6382 assert!(
6383 (force + 0.046_2).abs() < 5e-4 && (calibers - 0.752_2).abs() < 5e-4,
6384 "at the join's shallowest step the fallback moves the force {force:.4} and the center \
6385 of pressure {calibers:.4} calibres"
6386 );
6387 }
6388
6389 #[test]
6393 fn where_the_corners_turn_runs_out_the_join_carries_the_reading() {
6394 let model = model(&flared_rocket(18.5));
6395 let body = model.supersonic_body().unwrap();
6396 assert!(
6397 (body.join_start_mach - 1.555_220_046_128_9).abs() < 5e-13,
6398 "the table starts at Mach {}",
6399 body.join_start_mach
6400 );
6401 assert_eq!(body.shape_weight, 1.0);
6404 assert_eq!(body.weight(body.join_start_mach), 0.0);
6405 assert_eq!(
6406 body.weight(body.join_start_mach + SUPERSONIC_JOIN_WIDTH_MACH),
6407 1.0
6408 );
6409 let old = AeroModel::with_body_model(
6412 &flared_rocket(18.5).layout().unwrap(),
6413 BodyModel::CURRENT.with_supersonic_flare(SupersonicFlare::SlenderBody),
6414 )
6415 .unwrap();
6416 let at = |model: &AeroModel, mach| {
6417 model
6418 .normal_force(&flow(mach, 1e-4, 0.0))
6419 .unwrap()
6420 .coefficient
6421 / 1e-4
6422 };
6423 assert!((at(&model, body.join_start_mach) - at(&old, body.join_start_mach)).abs() < 1e-12);
6424 for (epsilon, want) in [
6425 (1e-9, 8.547_36e-10),
6426 (1e-7, 8.547_36e-8),
6427 (1e-5, 8.547_38e-6),
6428 ] {
6429 let below = at(&model, body.join_start_mach - epsilon);
6430 let above = at(&model, body.join_start_mach + epsilon);
6431 let gap = (above / below - 1.0).abs();
6432 assert!(
6433 (gap - want).abs() < 1e-5 * want,
6434 "±{epsilon} across the join's start moves the slope by {gap}, not {want}"
6435 );
6436 }
6437 }
6438
6439 #[test]
6450 fn the_cap_makes_a_kink_in_the_slope_even_though_the_reading_holds() {
6451 let ahead = flared_run(&flared_rocket(18.5));
6452 let boundary_deg = corner_limit_deg(&ahead, 2.0);
6453 let whole = |deg: f64| flared_at(&flared_rocket(deg), 2.0).0;
6454 let share = |deg: f64| {
6455 model(&flared_rocket(deg))
6456 .supersonic_body()
6457 .unwrap()
6458 .share(2, 2.0)
6459 .unwrap()
6460 .0
6461 };
6462 let kink = |f: &dyn Fn(f64) -> f64, at: f64| {
6463 let epsilon = 1e-6;
6464 let (down, up) = (
6465 (f(at) - f(at - epsilon)) / epsilon,
6466 (f(at + epsilon) - f(at)) / epsilon,
6467 );
6468 up / down - 1.0
6469 };
6470 for (what, f, at_boundary, at_20, at_26) in [
6471 (
6472 "the whole rocket",
6473 &whole as &dyn Fn(f64) -> f64,
6474 -0.313_9,
6475 0.035_0,
6476 0.0,
6477 ),
6478 ("the flare's share", &share, -1.415_6, 0.410_6, 0.0),
6479 ] {
6480 assert!(
6481 (kink(f, boundary_deg) - at_boundary).abs() < 5e-4,
6482 "{what} kinks {} at the boundary",
6483 kink(f, boundary_deg)
6484 );
6485 assert!(
6486 (kink(f, 20.0) - at_20).abs() < 5e-4,
6487 "{what} kinks {} at 20°",
6488 kink(f, 20.0)
6489 );
6490 assert!(
6492 (kink(f, 26.0) - at_26).abs() < 1e-6,
6493 "{what} kinks {} at 26°",
6494 kink(f, 26.0)
6495 );
6496 }
6497 }
6498
6499 #[test]
6523 fn a_near_flat_flare_reads_through_and_leaves_only_the_corners_crossing() {
6524 let join = |deg: f64| {
6525 model(&flared_rocket(deg))
6526 .supersonic_body()
6527 .map(|body| body.join_start_mach)
6528 };
6529 for deg in [
6532 0.0,
6533 2.4e-4,
6534 2.5e-4,
6535 3e-4,
6536 9.018_246_55e-4,
6537 0.001,
6538 0.01,
6539 0.03,
6540 0.038_161_270_2,
6541 0.045,
6542 0.05,
6543 0.058_820_517_4,
6544 0.06,
6545 0.1,
6546 1.0,
6547 ] {
6548 assert_eq!(
6549 join(deg),
6550 Some(SUPERSONIC_JOIN_START_MACH),
6551 "a {deg}° flare starts its table at {:?}",
6552 join(deg)
6553 );
6554 }
6555 let at = |deg: f64, mach: f64| {
6558 let model = model(&flared_rocket(deg));
6559 let force = model
6560 .normal_force(&flow(mach, 4f64.to_radians(), 0.0))
6561 .unwrap();
6562 (
6563 force.coefficient,
6564 force.cp_station_m.unwrap() / model.reference_diameter_m(),
6565 )
6566 };
6567 let across = |deg: f64, mach: f64, epsilon: f64| {
6568 let (below, above) = (at(deg - epsilon, mach), at(deg + epsilon, mach));
6569 (above.0 / below.0 - 1.0, above.1 - below.1)
6570 };
6571 for (deg, mach) in [(0.058_820_517_4, 3.0), (9.018_246_55e-4, 2.0)] {
6573 let (force, calibers) = across(deg, mach, 1e-9);
6574 assert!(
6575 force.abs() < 1e-7 && calibers.abs() < 1e-7,
6576 "crossing {deg}° at Mach {mach} moves the force {force:.3e} and the center of \
6577 pressure {calibers:.3e} calibres"
6578 );
6579 }
6580 let ahead = flared_run(&flared_rocket(1.0));
6584 let turns = |mach: f64| {
6585 crate::shock_expansion::flare_reduction_turns_rad(&ahead.aft_flow(mach).unwrap())
6586 .unwrap()
6587 };
6588 for (mach, force_want, caliber_want) in [
6589 (2.0_f64, 3.234e-6, -1.589e-6),
6590 (3.0, 1.1369e-4, 1.8258e-4),
6591 (4.0, 5.5377e-4, 1.6709e-3),
6592 (4.95, 1.2874e-3, 5.1095e-3),
6593 ] {
6594 let crossing = turns(mach).crossing_rad.to_degrees();
6595 let probes: &[f64] = if mach > 2.0 { &[1e-9, 1e-7] } else { &[1e-9] };
6598 for epsilon in probes {
6599 let (force, calibers) = across(crossing, mach, *epsilon);
6600 assert!(
6601 (force - force_want).abs() < 0.02 * force_want.abs()
6602 && (calibers - caliber_want).abs() < 0.02 * caliber_want.abs(),
6603 "Mach {mach}: ±{epsilon}° across the crossing at {crossing}° moves the force \
6604 {force:.4e} and the center of pressure {calibers:.4e} calibres"
6605 );
6606 }
6607 let balance = turns(mach).balance_rad.to_degrees();
6610 let (force, calibers) = across(balance, mach, 1e-9);
6611 assert!(
6612 force.abs() < 1e-7 && calibers.abs() < 1e-7,
6613 "Mach {mach}: the balance at {balance}° moves the force {force:.3e} and the \
6614 center of pressure {calibers:.3e} calibres"
6615 );
6616 }
6617 let machs = [3.0, 4.0];
6622 let step = 0.08 / 16.0;
6623 let swept: Vec<[f64; 2]> = (0..=16)
6624 .map(|index| {
6625 let model = model(&flared_rocket(index as f64 * step));
6626 machs.map(|mach| {
6627 model
6628 .normal_force(&flow(mach, 4f64.to_radians(), 0.0))
6629 .unwrap()
6630 .coefficient
6631 })
6632 })
6633 .collect();
6634 for (column, mach) in machs.iter().enumerate() {
6635 let crossing = turns(*mach).crossing_rad.to_degrees();
6636 let mut rises = Vec::new();
6637 for index in 1..swept.len() {
6638 let moved = swept[index][column] / swept[index - 1][column] - 1.0;
6639 let deg = index as f64 * step;
6640 assert!(
6641 moved.abs() < 2e-3,
6642 "Mach {mach}: from {}° to {deg}° the force moves {moved:.4}",
6643 deg - step
6644 );
6645 if moved > 0.0 {
6646 rises.push(deg);
6647 }
6648 }
6649 assert!(
6650 rises.len() <= 1
6651 && rises
6652 .first()
6653 .is_none_or(|deg| (deg - step..*deg).contains(&crossing)),
6654 "Mach {mach}: the force rises over {rises:?}, against a crossing at {crossing}°"
6655 );
6656 }
6657 }
6658
6659 #[test]
6666 fn issue_87s_switches_are_this_big() {
6667 let at = |rocket: &hpr_design::Rocket| {
6668 let model = model(rocket);
6669 let force = model
6670 .normal_force(&flow(3.0, 4f64.to_radians(), 0.0))
6671 .unwrap();
6672 (
6673 force.coefficient,
6674 force.cp_station_m.unwrap() / model.reference_diameter_m(),
6675 model.supersonic_body().is_some(),
6676 )
6677 };
6678 let stepped = |drop_m: f64| {
6680 let mut rocket = straight_rocket();
6681 rocket.stages[0].components[3].part = body_part(0.3, 0.027 - drop_m, 0.027 - drop_m);
6684 rocket
6685 };
6686 let flared = |rise_m: f64| {
6690 let mut rocket = crate::testing::finned_rocket(4);
6691 rocket.stages[0].components.truncate(3);
6692 rocket.stages[0].components.push(component(
6693 "flare",
6694 body_part(1.0, 0.022, 0.022 + rise_m),
6695 None,
6696 ));
6697 rocket
6698 };
6699 let coned = |half_angle_deg: f64| {
6701 let mut rocket = straight_rocket();
6702 rocket.stages[0].components[0].part = nose(
6703 NoseShape::Conical {},
6704 0.027 / half_angle_deg.to_radians().tan(),
6705 0.027,
6706 );
6707 rocket
6708 };
6709 let blunt = |length_m: f64| {
6711 let mut rocket = straight_rocket();
6712 rocket.stages[0].components[0].part =
6713 nose(NoseShape::PowerSeries { exponent: 0.5 }, length_m, 0.027);
6714 rocket
6715 };
6716 let at_16 = 0.027 / 24.0_f64.to_radians().tan();
6717 assert!(model(&stepped(2.6e-11)).supersonic_body().is_some());
6721 assert!(model(&stepped(2.8e-11)).supersonic_body().is_none());
6722 type Side = (f64, f64, bool);
6726 let switches: [(&str, Side, Side, (f64, f64)); 4] = [
6727 (
6728 "a step in radius",
6729 at(&stepped(1e-12)),
6730 at(&stepped(1e-9)),
6731 (-0.0865, 1.0285),
6732 ),
6733 (
6734 "a flare behind a boattail",
6735 at(&flared(1e-12)),
6736 at(&flared(1e-9)),
6737 (-0.2749, 0.2872),
6738 ),
6739 (
6740 "a pointed tip past the cone tables' 30°",
6741 at(&coned(29.999)),
6742 at(&coned(30.002)),
6743 (-0.0770, 0.8107),
6744 ),
6745 (
6746 "a vertical tip steeper than the handover",
6747 at(&blunt(0.5 * at_16 * 1.0002)),
6748 at(&blunt(0.5 * at_16 * 0.9998)),
6749 (-0.0699, 0.6383),
6750 ),
6751 ];
6752 let (below, above) = (at(&coned(23.999)), at(&coned(24.002)));
6755 assert!(
6756 below.2 && above.2,
6757 "a 24° tip flies the method on both sides now"
6758 );
6759 assert!(
6760 (above.0 / below.0 - 1.0).abs() < 1e-4 && (above.1 - below.1).abs() < 1e-4,
6761 "across Fig. 2's old edge: {below:?} to {above:?}"
6762 );
6763 for (what, covered, bare, (want_force, want_calibers)) in switches {
6764 assert!(
6765 covered.2 && !bare.2,
6766 "{what}: the method should cover one side only ({covered:?}, {bare:?})"
6767 );
6768 let force = bare.0 / covered.0 - 1.0;
6769 let calibers = bare.1 - covered.1;
6770 assert!(
6771 (force - want_force).abs() < 5e-4 && (calibers - want_calibers).abs() < 5e-4,
6772 "{what}: the force moves {force:.4} and the center of pressure {calibers:.4} \
6773 calibres, against {want_force} and {want_calibers}"
6774 );
6775 }
6776 }
6777
6778 #[test]
6779 fn a_steep_tip_is_reported_as_the_runs_fallback() {
6780 let with_tip = |half_angle_deg: f64| {
6783 let mut rocket = straight_rocket();
6784 let Part::NoseCone(nose) = &mut rocket.stages[0].components[0].part else {
6785 panic!("the first part is the nose");
6786 };
6787 nose.shape = NoseShape::Conical {};
6788 nose.length_m = nose.base_radius_m / half_angle_deg.to_radians().tan();
6789 model(&rocket).supersonic_fallback()
6790 };
6791 assert_eq!(with_tip(29.99), None);
6792 assert_eq!(with_tip(30.01), Some(SupersonicFallback::SteepTip));
6793 }
6794
6795 #[test]
6796 fn a_step_is_reported_as_the_runs_fallback() {
6797 let stepped = |drop_m: f64| {
6800 let mut rocket = straight_rocket();
6801 for (i, length_m) in [0.0, 0.7, 0.05, 0.3].iter().enumerate().skip(2) {
6802 rocket.stages[0].components[i].part =
6803 body_part(*length_m, 0.027 - drop_m, 0.027 - drop_m);
6804 }
6805 let id = rocket.stages[0].components[2].id.clone();
6806 (model(&rocket).supersonic_fallback(), id)
6807 };
6808 assert_eq!(stepped(0.0).0, None);
6809 let (fallback, id) = stepped(1e-3);
6810 assert_eq!(
6811 fallback,
6812 Some(SupersonicFallback::RadiusStep { component: id })
6813 );
6814 let (fallback, id) = stepped(-1e-3);
6815 assert_eq!(
6816 fallback,
6817 Some(SupersonicFallback::RadiusStep { component: id })
6818 );
6819 assert_eq!(stepped(2.6e-11).0, None);
6823 for drop_m in [2.8e-11, 1e-9, -1e-9, 1e-8] {
6824 let (fallback, id) = stepped(drop_m);
6825 assert_eq!(
6826 fallback,
6827 Some(SupersonicFallback::RadiusStep { component: id }),
6828 "a step of {drop_m} m"
6829 );
6830 }
6831 }
6832
6833 #[test]
6849 fn a_step_takes_the_whole_body_off_the_method() {
6850 use hpr_design::ReferenceDiameter;
6851 let at_joint = |joint: usize, drop_m: f64| {
6858 let mut rocket = straight_rocket();
6859 rocket.reference_diameter = ReferenceDiameter::Custom { diameter_m: 0.054 };
6860 let lengths = [0.0, 0.7, 0.05, 0.3];
6861 for (i, length_m) in lengths.iter().enumerate().skip(joint) {
6862 rocket.stages[0].components[i].part =
6863 body_part(*length_m, 0.027 - drop_m, 0.027 - drop_m);
6864 }
6865 rocket
6866 };
6867 let read = |joint: usize, drop_m: f64| {
6868 let rocket = at_joint(joint, drop_m);
6869 let model = model(&rocket);
6870 let force = model
6871 .normal_force(&flow(3.0, 4f64.to_radians(), 0.0))
6872 .unwrap();
6873 (
6874 force.coefficient,
6875 force.cp_station_m.unwrap() / model.reference_diameter_m(),
6876 model.supersonic_body().is_some(),
6877 )
6878 };
6879 let (flush, flush_cp, marched) = read(3, 0.0);
6880 assert!(marched, "with no step the method covers the body");
6881 assert!(
6882 (flush - 0.899_591_695).abs() < 5e-7 && (flush_cp - 16.949_177_4).abs() < 5e-5,
6883 "the flush rocket reads {flush} at {flush_cp} calibres"
6884 );
6885
6886 for (joint, drop_m, want_force, want_calibers) in [
6892 (1, 2.8e-11, -0.086_518_790, 1.028_480_4),
6893 (2, 2.8e-11, -0.086_518_790, 1.028_480_4),
6894 (3, 2.8e-11, -0.086_518_789, 1.028_480_4),
6895 (1, 1e-3, -0.106_197_868, 1.193_751_6),
6896 (3, 1e-3, -0.102_873_512, 0.994_855_8),
6897 (1, 2e-3, -0.125_539_709, 1.359_310_7),
6898 (3, 2e-3, -0.118_890_998, 0.959_745_0),
6899 (1, -2.8e-11, -0.086_518_788, 1.028_480_4),
6900 (3, -2.8e-11, -0.086_518_789, 1.028_480_4),
6901 (1, -1e-3, -0.067_194_034, 0.867_387_3),
6902 (3, -1e-3, -0.070_503_028, 1.064_443_7),
6903 (1, -2e-3, -0.047_519_789, 0.706_432_7),
6904 (3, -2e-3, -0.054_109_172, 1.098_652_0),
6905 ] {
6906 let (force, cp, marched) = read(joint, drop_m);
6907 assert!(
6908 !marched,
6909 "a {drop_m} m step at joint {joint} kept the method"
6910 );
6911 let (force, calibers) = (force / flush - 1.0, cp - flush_cp);
6912 assert!(
6913 (force - want_force).abs() < 5e-6 && (calibers - want_calibers).abs() < 5e-5,
6914 "a {drop_m} m step at joint {joint} moves the force {force:.9} and the center of \
6915 pressure {calibers:.7} calibres, against {want_force} and {want_calibers}"
6916 );
6917 }
6918
6919 let (mut low, mut high) = (1e-13, 1e-8);
6925 assert!(
6926 read(3, low).2 && !read(3, high).2,
6927 "the bisection has to start with the march covering one end and refusing the other"
6928 );
6929 while high - low > 1e-6 * high {
6930 let mid = 0.5 * (low + high);
6931 if read(3, mid).2 {
6932 low = mid
6933 } else {
6934 high = mid
6935 }
6936 }
6937 assert!(
6938 (high / 2.7e-11 - 1.0).abs() < 2e-6,
6939 "the march refuses a step of {high} m, not a billionth of the 0.027 m radius"
6940 );
6941 for joint in [1, 2, 3] {
6942 for sign in [1.0, -1.0] {
6943 assert!(
6944 read(joint, sign * 2.6e-11).2 && !read(joint, sign * 2.8e-11).2,
6945 "at joint {joint} with sign {sign} the switch is not at 2.7e-11 m"
6946 );
6947 }
6948 }
6949
6950 let run_segments = |drop_m: f64| {
6955 let rocket = at_joint(3, drop_m);
6956 model(&rocket)
6957 .supersonic_run
6958 .as_ref()
6959 .map(|run| run.segments.len())
6960 };
6961 for drop_m in [1e-8, -1e-8] {
6962 assert_eq!(
6963 run_segments(drop_m),
6964 Some(4),
6965 "{drop_m} m is inside the gate"
6966 );
6967 }
6968 for drop_m in [2e-8, -2e-8] {
6969 assert_eq!(run_segments(drop_m), None, "{drop_m} m is outside the gate");
6970 }
6971
6972 let boattail = |drop_m: f64| {
6980 let mut rocket = crate::testing::finned_rocket(4);
6981 rocket.reference_diameter = ReferenceDiameter::Custom { diameter_m: 0.054 };
6982 rocket.stages[0].components[2].part = body_part(0.05, 0.027 - drop_m, 0.022);
6983 let model = model(&rocket);
6984 let force = model
6985 .normal_force(&flow(3.0, 4f64.to_radians(), 0.0))
6986 .unwrap();
6987 (
6988 force.coefficient,
6989 force.cp_station_m.unwrap() / model.reference_diameter_m(),
6990 model.supersonic_body().is_some(),
6991 )
6992 };
6993 let (boattail_flush, boattail_flush_cp, marched) = boattail(0.0);
6994 assert!(marched, "flush, the method covers the boattailed body too");
6995 let (mut low, mut high) = (1e-18, 1e-8);
6996 assert!(boattail(-low).2 && !boattail(-high).2);
6997 while high - low > 1e-6 * high {
6998 let mid = 0.5 * (low + high);
6999 if boattail(-mid).2 {
7000 low = mid
7001 } else {
7002 high = mid
7003 }
7004 }
7005 assert!(
7006 (high / (1e-12 * 1.3 * 0.1) - 1.0).abs() < 2e-4,
7007 "a step up at the boattail's joint is refused at {high} m, not 1e-12 × 1.3 m × 0.1"
7008 );
7009 assert!(
7010 boattail(2.6e-11).2 && !boattail(2.8e-11).2,
7011 "stepping down at that joint, the merge tolerance still binds"
7012 );
7013 for drop_m in [-1.4e-13, -2.8e-11, 2.8e-11] {
7014 let (force, cp, marched) = boattail(drop_m);
7015 assert!(
7016 !marched,
7017 "a {drop_m} m step at the boattail kept the method"
7018 );
7019 let (force, calibers) = (force / boattail_flush - 1.0, cp - boattail_flush_cp);
7020 assert!(
7021 (force + 0.113_409_121).abs() < 5e-6 && (calibers - 1.095_116_4).abs() < 5e-5,
7022 "a {drop_m} m step at the boattail moves the force {force:.9} and the center of \
7023 pressure {calibers:.7} calibres"
7024 );
7025 }
7026 }
7027
7028 #[test]
7029 fn a_body_the_method_cannot_finish_keeps_slender_body_terms() {
7030 let at = |m: &AeroModel, mach| body_values(m, mach);
7031 let mut rocket = straight_rocket();
7036 rocket.stages[0].components[2].part = body_part(0.05, 0.027, 0.032);
7039 if let Part::Transition(transition) = &mut rocket.stages[0].components[2].part {
7040 transition.shape = NoseShape::Ogive { radius_ratio: 1.0 };
7041 }
7042 rocket.stages[0].components[3].part = body_part(0.3, 0.032, 0.032);
7043 let flared = model(&rocket);
7044 let mut rocket = crate::testing::finned_rocket(4);
7047 rocket.stages[0]
7048 .components
7049 .push(component("lip", body_part(0.01, 0.022, 0.025), None));
7050 let lipped = model(&rocket);
7051 let mut rocket = crate::testing::finned_rocket(4);
7052 rocket.stages[0].components[2].part = body_part(0.05, 0.026, 0.022);
7053 let stepped_boattail = model(&rocket);
7054 let mut rocket = straight_rocket();
7055 rocket.stages[0].components[0].part =
7056 nose(NoseShape::PowerSeries { exponent: 0.5 }, 0.027, 0.027);
7057 let blunt = model(&rocket);
7058 let mut rocket = straight_rocket();
7059 rocket.stages[0].components[1].part = body_part(0.7, 0.03, 0.03);
7060 let stepped = model(&rocket);
7061 for model in [&flared, &lipped, &stepped_boattail, &blunt, &stepped] {
7062 assert!(model.supersonic_body().is_none());
7063 assert_eq!(at(model, 3.0), at(model, 0.5));
7064 }
7065 }
7066
7067 #[test]
7068 fn clones_share_the_table_and_compare_equal() {
7069 let model = model(&straight_rocket());
7070 let clone = model.clone();
7071 assert_eq!(model, clone);
7072 let built = model.supersonic_body().unwrap() as *const SupersonicBody;
7073 assert_eq!(model, clone);
7074 assert!(std::ptr::eq(built, clone.supersonic_body().unwrap()));
7075 }
7076
7077 #[test]
7081 fn models_built_apart_share_the_table_only_on_the_same_body() {
7082 let rocket = straight_rocket();
7083 let nominal = model(&rocket);
7084 let alone = model(&rocket).with_drag_scale(1.1).unwrap();
7085 let built = nominal.supersonic_body().unwrap() as *const SupersonicBody;
7086 assert!(!std::ptr::eq(built, alone.supersonic_body().unwrap()));
7087 let mut scaled = model(&rocket).with_drag_scale(1.1).unwrap();
7088 assert!(!scaled.supersonic_table_built());
7089 assert!(scaled.share_supersonic_table(&nominal));
7090 assert!(scaled.supersonic_table_built());
7091 assert!(std::ptr::eq(built, scaled.supersonic_body().unwrap()));
7092 let first = model(&rocket);
7094 let mut second = model(&rocket);
7095 assert!(second.share_supersonic_table(&first));
7096 assert!(!first.supersonic_table_built());
7097 let built = second.supersonic_body().unwrap() as *const SupersonicBody;
7098 assert!(first.supersonic_table_built());
7099 assert!(std::ptr::eq(built, first.supersonic_body().unwrap()));
7100 let mut longer = rocket.clone();
7102 longer.stages[0].components[0].part =
7103 nose(NoseShape::Ogive { radius_ratio: 1.0 }, 0.4, 0.027);
7104 let mut longer = model(&longer);
7105 assert!(!longer.share_supersonic_table(&nominal));
7106 assert!(!std::ptr::eq(
7107 nominal.supersonic_body().unwrap(),
7108 longer.supersonic_body().unwrap()
7109 ));
7110 let mut wider = rocket.clone();
7112 wider.reference_diameter = ReferenceDiameter::Custom { diameter_m: 0.07 };
7113 let mut wider = model(&wider);
7114 assert!(!wider.share_supersonic_table(&nominal));
7115 assert!(!wider.supersonic_table_built());
7116 let mut blunt = rocket;
7117 blunt.stages[0].components[0].part =
7118 nose(NoseShape::PowerSeries { exponent: 0.5 }, 0.027, 0.027);
7119 let mut blunt = model(&blunt);
7120 assert!(!blunt.share_supersonic_table(&nominal));
7121 assert!(blunt.supersonic_body().is_none());
7122 }
7123
7124 fn shares(rocket: &hpr_design::Rocket) -> (Vec<(String, Drag)>, f64) {
7126 let m = model(rocket);
7127 let conditions = DragConditions::coasting(0.3 * 340.294 / 1.4607e-5);
7128 let parts: Vec<(String, Drag)> = m
7129 .buildup_components(&Flow::axial(0.3), &conditions)
7130 .unwrap()
7131 .into_iter()
7132 .map(|part| (part.id, part.drag))
7133 .collect();
7134 let total = m.drag(&Flow::axial(0.3), &conditions).unwrap();
7135 let sum: f64 = parts.iter().map(|(_, d)| d.zero_lift_coefficient).sum();
7136 close(total.zero_lift_coefficient, sum, 1e-12, "the parts add up");
7137 (parts, total.zero_lift_coefficient)
7138 }
7139
7140 fn share(parts: &[(String, Drag)], id: &str) -> Drag {
7141 parts
7142 .iter()
7143 .find(|(i, _)| i == id)
7144 .map(|(_, d)| *d)
7145 .unwrap()
7146 }
7147
7148 fn stated(coefficient: f64, include_children: bool) -> Option<hpr_design::DragOverride> {
7149 Some(hpr_design::DragOverride {
7150 coefficient,
7151 include_children,
7152 })
7153 }
7154
7155 fn stepped_rocket() -> hpr_design::Rocket {
7158 let mut rocket = finned_rocket(3);
7159 let components = &mut rocket.stages[0].components;
7160 components.remove(2);
7161 rocket
7162 }
7163
7164 #[test]
7168 fn a_stated_drag_coefficient_replaces_the_part_s_own() {
7169 let plain = stepped_rocket();
7170 let (before, total) = shares(&plain);
7171 let step = base_drag_coefficient(0.3).unwrap() * (1.0 - (0.022_f64 / 0.027).powi(2));
7172 assert!(
7174 share(&before, "tail").pressure >= step,
7175 "{:?}",
7176 share(&before, "tail")
7177 );
7178
7179 let mut fins = plain.clone();
7181 fins.stages[0].components[2].children[0].drag_override = stated(0.5, false);
7182 let (after, with) = shares(&fins);
7183 close(share(&after, "fins").stated, 1.5, 1e-15, "per fin");
7184 close(
7185 with,
7186 total - share(&before, "fins").zero_lift_coefficient + 1.5,
7187 1e-12,
7188 "fins",
7189 );
7190
7191 let mut body = plain.clone();
7193 body.stages[0].components[1].drag_override = stated(0.2, false);
7194 let (after, with) = shares(&body);
7195 close(
7196 share(&after, "body").zero_lift_coefficient,
7197 0.2,
7198 1e-15,
7199 "body",
7200 );
7201 close(
7202 share(&before, "tail").zero_lift_coefficient
7203 - share(&after, "tail").zero_lift_coefficient,
7204 step,
7205 1e-12,
7206 "the step goes with the body",
7207 );
7208 close(
7209 with,
7210 total - share(&before, "body").zero_lift_coefficient - step + 0.2,
7211 1e-12,
7212 "body",
7213 );
7214
7215 let mut tail = plain.clone();
7217 tail.stages[0].components[2].drag_override = stated(0.2, false);
7218 let (after, with_tail) = shares(&tail);
7219 let own = share(&after, "tail");
7220 close(own.stated, 0.2, 1e-15, "tail");
7221 close(own.pressure, step, 1e-12, "the step stays");
7222 assert_eq!((own.friction, own.base, own.parasitic), (0.0, 0.0, 0.0));
7223 assert_eq!(share(&after, "fins"), share(&before, "fins"));
7224 close(
7225 with_tail,
7226 total - share(&before, "tail").zero_lift_coefficient + step + 0.2,
7227 1e-12,
7228 "tail",
7229 );
7230
7231 tail.stages[0].components[2].drag_override = stated(0.2, true);
7233 let (after, _) = shares(&tail);
7234 assert_eq!(share(&after, "fins").zero_lift_coefficient, 0.0);
7235
7236 let mut rails = plain.clone();
7239 rails.stages[0].components[1].children = vec![
7240 component(
7241 "lug",
7242 Part::LaunchLug(LaunchLug {
7243 length_m: 0.05,
7244 outer_radius_m: 0.004,
7245 thickness_m: 0.0005,
7246 angle_rad: 0.0,
7247 count: 1,
7248 spacing_m: 0.0,
7249 material: material(),
7250 }),
7251 Some(Position::Middle { aft_offset_m: 0.0 }),
7252 ),
7253 component(
7254 "buttons",
7255 Part::RailButton(RailButton {
7256 outer_diameter_m: 0.01,
7257 inner_diameter_m: 0.006,
7258 height_m: 0.008,
7259 base_height_m: 0.002,
7260 flange_height_m: 0.002,
7261 screw_height_m: 0.0,
7262 angle_rad: 0.0,
7263 count: 2,
7264 spacing_m: 0.2,
7265 material: material(),
7266 }),
7267 Some(Position::Middle { aft_offset_m: 0.0 }),
7268 ),
7269 ];
7270 let (before_rails, total_rails) = shares(&rails);
7271 for id in ["lug", "buttons"] {
7272 assert!(share(&before_rails, id).zero_lift_coefficient > 0.0, "{id}");
7273 }
7274 rails.stages[0].components[1].children[1].drag_override = stated(0.5, false);
7275 let (after, _) = shares(&rails);
7276 close(share(&after, "buttons").stated, 1.0, 1e-15, "per button");
7277 rails.stages[0].components[1].children[1].drag_override = None;
7278 rails.stages[0].components[1].drag_override = stated(0.2, true);
7279 let (after, with) = shares(&rails);
7280 for id in ["lug", "buttons"] {
7281 assert_eq!(share(&after, id).zero_lift_coefficient, 0.0, "{id}");
7282 }
7283 close(
7284 with,
7285 total_rails
7286 - share(&before_rails, "body").zero_lift_coefficient
7287 - share(&before_rails, "lug").zero_lift_coefficient
7288 - share(&before_rails, "buttons").zero_lift_coefficient
7289 - step
7290 + 0.2,
7291 1e-12,
7292 "body covering its lug and buttons",
7293 );
7294
7295 let mut stage = plain.clone();
7297 stage.stages[0].drag_override = stated(0.1, false);
7298 let (after, with) = shares(&stage);
7299 close(share(&after, "stage").stated, 0.1, 1e-15, "stage");
7300 close(with, total + 0.1, 1e-12, "stage added");
7301 stage.stages[0].drag_override = stated(0.1, true);
7302 let (_, with) = shares(&stage);
7303 close(with, 0.1, 1e-15, "stage alone");
7304
7305 let m = model(&stage).with_drag_scale(2.0).unwrap();
7307 let conditions = DragConditions::coasting(0.3 * 340.294 / 1.4607e-5);
7308 let scaled = m.drag(&Flow::axial(0.3), &conditions).unwrap();
7309 close(scaled.zero_lift_coefficient, 0.2, 1e-15, "scaled");
7310 }
7311
7312 #[test]
7315 fn a_stated_drag_coefficient_out_of_range_or_in_a_pod_is_refused() {
7316 for bad in [-0.1, f64::NAN, f64::INFINITY] {
7317 let mut rocket = stepped_rocket();
7318 rocket.stages[0].components[1].drag_override = stated(bad, false);
7319 match AeroModel::new(&rocket.layout().unwrap()) {
7320 Err(AeroError::InComponent { id, source }) => {
7321 assert_eq!(id, "body");
7322 assert!(
7323 matches!(
7324 *source,
7325 AeroError::Domain {
7326 what: "stated drag coefficient",
7327 ..
7328 }
7329 ),
7330 "{source:?}"
7331 );
7332 }
7333 other => panic!("{bad}: {other:?}"),
7334 }
7335 }
7336 let mut rocket = stepped_rocket();
7337 rocket.stages[0].drag_override = stated(-1.0, true);
7338 assert!(matches!(
7339 AeroModel::new(&rocket.layout().unwrap()),
7340 Err(AeroError::InComponent { ref id, .. }) if id == "stage"
7341 ));
7342 let mut rocket = stepped_rocket();
7345 rocket.stages[0].components[2].drag_override = stated(0.1, true);
7346 rocket.stages[0].components[2].children[0].drag_override = stated(f64::NAN, false);
7347 assert!(matches!(
7348 AeroModel::new(&rocket.layout().unwrap()),
7349 Err(AeroError::InComponent { ref id, .. }) if id == "fins"
7350 ));
7351 let refused = |rocket: &hpr_design::Rocket, why: &str| match AeroModel::new(
7352 &rocket.layout().unwrap(),
7353 ) {
7354 Err(AeroError::InComponent { source, .. }) => {
7355 assert!(
7356 matches!(*source, AeroError::Unsupported(ref what) if what == why),
7357 "{source:?}"
7358 );
7359 }
7360 other => panic!("{why}: {other:?}"),
7361 };
7362 for (pod_set, part) in [(true, 0), (false, 1)] {
7363 let mut rocket = podded_rocket(2);
7364 let pods = rocket.stages[0].components[1].children.last_mut().unwrap();
7365 if pod_set {
7366 pods.drag_override = stated(0.1, true);
7367 } else {
7368 pods.children[part].drag_override = stated(0.1, false);
7369 }
7370 refused(&rocket, "a drag override on a pod set or in a pod");
7371 }
7372 let mut rocket = podded_rocket(2);
7374 rocket.stages[0].components[1].drag_override = stated(0.1, true);
7375 refused(&rocket, "a drag override covering a pod set");
7376 let mut rocket = podded_rocket(2);
7377 rocket.stages[0].drag_override = stated(0.1, true);
7378 refused(&rocket, "a drag override covering a pod set");
7379 let mut rocket = podded_rocket(2);
7381 rocket.stages[0].components[1].drag_override = stated(0.1, false);
7382 assert!(AeroModel::new(&rocket.layout().unwrap()).is_ok());
7383 let mut rocket = tube_finned_rocket(6, 0.0762, 0.011, 0.0005);
7385 rocket.stages[0].components[3].children[1].drag_override = stated(0.1, false);
7386 refused(&rocket, "a drag override on a tube fin set");
7387 }
7388
7389 #[test]
7392 fn a_looped_layout_does_not_hang_the_override_walk() {
7393 let mut lugged = stepped_rocket();
7394 lugged.stages[0].components[1].children.push(component(
7395 "lug",
7396 Part::LaunchLug(LaunchLug {
7397 length_m: 0.05,
7398 outer_radius_m: 0.004,
7399 thickness_m: 0.0005,
7400 angle_rad: 0.0,
7401 count: 1,
7402 spacing_m: 0.0,
7403 material: material(),
7404 }),
7405 Some(Position::Middle { aft_offset_m: 0.0 }),
7406 ));
7407 let mut layout = lugged.layout().unwrap();
7408 let (lug, _) = layout.find("lug").unwrap();
7409 layout.components[lug].parent = Some(lug);
7410 let _ = AeroModel::new(&layout);
7412 }
7413
7414 #[test]
7416 fn a_stated_drag_coefficient_flies_a_shape_the_buildup_refuses() {
7417 let mut rocket = stepped_rocket();
7418 rocket.stages[0].components[0].part =
7420 nose(NoseShape::Ogive { radius_ratio: 0.8 }, 0.25, 0.027);
7421 let conditions = DragConditions::coasting(0.3 * 340.294 / 1.4607e-5);
7422 assert!(matches!(
7423 model(&rocket).drag(&Flow::axial(0.3), &conditions),
7424 Err(AeroError::InComponent { .. })
7425 ));
7426 rocket.stages[0].components[0].drag_override = stated(0.3, false);
7427 let (parts, _) = shares(&rocket);
7428 close(
7429 share(&parts, "nose").zero_lift_coefficient,
7430 0.3,
7431 1e-15,
7432 "nose",
7433 );
7434 }
7435}