Reduce abstraction for lift calculation; remove RigidWings struct

common/src/comp/fluid_dynamics.rs

@@ -1,6 +1,6 @@
 use super::{
     body::{object, Body},
-    Density, Ori, Vel,
+    CharacterState, Density, Ori, Vel,
 };
 use crate::{
     consts::{AIR_DENSITY, WATER_DENSITY},
@@ -90,10 +90,9 @@ impl Default for Fluid {
 impl Body {
     pub fn aerodynamic_forces(
         &self,
-        ori: &Ori,
         rel_flow: &Vel,
         fluid_density: f32,
-        wings: Option<&RigidWings>,
+        character_state: Option<&CharacterState>,
     ) -> Vec3<f32> {
         let v_sq = rel_flow.0.magnitude_squared();
         if v_sq < 0.25 {
@@ -104,39 +103,41 @@ impl Body {
             // All the coefficients come pre-multiplied by their reference area
             0.5 * fluid_density
                 * v_sq
-                * wings
+                * character_state
-                    .map(|wings| {
+                    .and_then(|cs| match cs {
-                        // Since we have wings, we proceed to calculate the lift and drag
+                        CharacterState::Glide(data) => {
+                            Some((data.aspect_ratio, data.planform_area, data.ori))
+                        },
+                        _ => None,
+                    })
+                    .map(|(ar, area, ori)| {
+                        // We have an elliptical wing; proceed to calculate its lift and drag
 
-                        let ar = wings.aspect_ratio();
                         // aoa will be positive when we're pitched up and negative otherwise
-                        let aoa = angle_of_attack(ori, &rel_flow_dir);
+                        let aoa = angle_of_attack(&ori, &rel_flow_dir);
                         // c_l will be positive when aoa is positive (we have positive lift,
                         // producing an upward force) and negative otherwise
-                        let c_l = wings.lift_coefficient(aoa);
+                        let c_l = lift_coefficient(ar, area, aoa);
 
                         // lift dir will be orthogonal to the local relative flow vector.
-                        // Local relative flow is the resulting vector of (relative) freestream flow
+                        // Local relative flow is the resulting vector of (relative) freestream
-                        // + downwash (created by the vortices of the wing tips)
+                        // flow + downwash (created by the vortices
+                        // of the wing tips)
                         let lift_dir: Dir = {
                             // induced angle of attack
                             let aoa_i = c_l / (PI * ar);
-                            // effective angle of attack; the aoa as seen by aerofoil after downwash
+                            // effective angle of attack; the aoa as seen by aerofoil after
+                            // downwash
                             let aoa_eff = aoa - aoa_i;
-                            /*println!(
+                            // Angle between chord line and local relative wind is aoa_eff
-                                "CL={:.1}, α={:.1}°, αᵢ={:.1}°, αₑ={:.1}°, AR={:.1}",
+                            // radians. Direction of lift is
-                                c_l,
+                            // perpendicular to local relative wind.
-                                aoa.to_degrees(),
+                            // At positive lift, local relative wind will be below our cord line
-                                aoa_i.to_degrees(),
+                            // at an angle of aoa_eff. Thus if
-                                aoa_eff.to_degrees(),
+                            // we pitch down by aoa_eff radians then
-                                ar
+                            // our chord line will be colinear with local relative wind vector
-                            );*/
+                            // and our up will be the direction
-                            // Angle between chord line and local relative wind is aoa_eff radians.
+                            // of lift.
-                            // Direction of lift is perpendicular to local relative wind.
-                            // At positive lift, local relative wind will be below our cord line at
-                            // an angle of aoa_eff. Thus if we pitch down by aoa_eff radians then
-                            // our chord line will be colinear with local relative wind vector and
-                            // our up will be the direction of lift.
                             ori.pitched_down(aoa_eff).up()
                         };
 
@@ -146,20 +147,12 @@ impl Body {
                             // (this definition should not be used for aspect ratios > 25)
                             let e = 1.78 * (1.0 - 0.045 * ar.powf(0.68)) - 0.64;
 
-                            wings.zero_lift_drag_coefficient()
+                            zero_lift_drag_coefficient(area)
                                 + self.parasite_drag_coefficient()
                                 + c_l.powi(2) / (PI * e * ar)
                         };
                         debug_assert!(c_d.is_sign_positive());
                         debug_assert!(c_l.is_sign_positive() || aoa.is_sign_negative());
-                        /*println!(
-                            "L/D (at α={:.1}, AR={:.1}) = {:.1}/{:.1} = {:.1}",
-                            aoa.to_degrees(),
-                            ar,
-                            0.5 * fluid_density * v_sq * c_l,
-                            0.5 * fluid_density * v_sq * c_d,
-                            c_l / c_d
-                        );*/
 
                         c_l * *lift_dir + c_d * *rel_flow_dir
                     })
@@ -192,24 +185,28 @@ impl Body {
                 } else {
                     1.0
                 };
-                cd * std::f32::consts::PI * dim.x * dim.z
+                cd * PI * dim.x * dim.z
             },
 
             // Cross-section, zero-lift angle; exclude the wings (width * 0.2)
             Body::BirdMedium(_) | Body::BirdLarge(_) | Body::Dragon(_) => {
                 let dim = self.dimensions().map(|a| a * 0.5);
+                // "Field Estimates of Body Drag Coefficient on the Basis of Dives in Passerine
+                // Birds", Anders Hedenström and Felix Liechti, 2001
                 let cd = match self {
-                    Body::BirdMedium(_) => 0.2,
+                    Body::BirdLarge(_) | Body::BirdMedium(_) => 0.2,
-                    Body::BirdLarge(_) => 0.4,
+                    // arbitrary
                     _ => 0.7,
                 };
-                cd * std::f32::consts::PI * dim.x * 0.2 * dim.z
+                cd * PI * dim.x * 0.2 * dim.z
             },
 
             // Cross-section, zero-lift angle; exclude the fins (width * 0.2)
             Body::FishMedium(_) | Body::FishSmall(_) => {
                 let dim = self.dimensions().map(|a| a * 0.5);
-                0.031 * std::f32::consts::PI * dim.x * 0.2 * dim.z
+                // "A Simple Method to Determine Drag Coefficients in Aquatic Animals",
+                // D. Bilo and W. Nachtigall, 1980
+                0.031 * PI * dim.x * 0.2 * dim.z
             },
 
             Body::Object(object) => match object {
@@ -225,7 +222,7 @@ impl Body {
                 | object::Body::FireworkYellow
                 | object::Body::MultiArrow => {
                     let dim = self.dimensions().map(|a| a * 0.5);
-                    0.02 * std::f32::consts::PI * dim.x * dim.z
+                    0.02 * PI * dim.x * dim.z
                 },
 
                 // spherical-ish objects
@@ -243,12 +240,12 @@ impl Body {
                 | object::Body::Pumpkin4
                 | object::Body::Pumpkin5 => {
                     let dim = self.dimensions().map(|a| a * 0.5);
-                    0.5 * std::f32::consts::PI * dim.x * dim.z
+                    0.5 * PI * dim.x * dim.z
                 },
 
                 _ => {
                     let dim = self.dimensions();
-                    2.0 * (std::f32::consts::PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0)
+                    2.0 * (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0)
                 },
             },
 
@@ -256,7 +253,7 @@ impl Body {
                 // Airships tend to use the square of the cube root of its volume for
                 // reference area
                 let dim = self.dimensions();
-                (std::f32::consts::PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0)
+                (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0)
             },
         }
     }
@@ -267,122 +264,73 @@ fn angle_of_attack(ori: &Ori, rel_flow_dir: &Dir) -> f32 {
     PI / 2.0 - ori.up().angle_between(rel_flow_dir.to_vec())
 }
 
-/// An elliptical fixed rigid wing. Plurally named simply because it's a shape
+/// Total lift coefficient for a finite wing of symmetric aerofoil shape and
-/// typically composed of two wings forming an elliptical lift distribution.
+/// elliptical pressure distribution.
-//
+pub fn lift_coefficient(aspect_ratio: f32, planform_area: f32, aoa: f32) -> f32 {
-// Animal wings are technically flexible, not rigid, (difference being that the
+    let aoa_abs = aoa.abs();
-// former's shape is affected by the flow) and usually has the ability to
+    let stall_angle = PI * 0.1;
-// assume complex shapes with properties like curved camber line, span-wise
+    inline_tweak::tweak!(1.0)
-// twist, dihedral angle, sweep angle, and partitioned sections. However, we
+        * planform_area
-// could make do with this model for fully extended animal wings, enabling them
+        * if aoa_abs < stall_angle {
-// to glide.
+            lift_slope(aspect_ratio, None) * aoa
-#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
+        } else if inline_tweak::tweak!(true) {
-pub struct RigidWings {
+            // This is when flow separation and turbulence starts to kick in.
-    aspect_ratio: f32,
+            // Going to just make something up (based on some data), as the alternative is
-    planform_area: f32,
+            // to just throw your hands up and return 0
-    // sweep_angle: Option<f32>,
+            let aoa_s = aoa.signum();
-}
+            let c_l_max = lift_slope(aspect_ratio, None) * stall_angle;
-
+            let deg_45 = PI / 4.0;
-impl RigidWings {
+            if aoa_abs < deg_45 {
-    /// Wings from total span (wing-tip to wing-tip) and
+                // drop directly to 0.6 * max lift at stall angle
-    /// chord length (leading edge to trailing edge)
+                // then climb back to max at 45°
-    pub fn new(span_length: f32, chord_length: f32) -> Self {
+                Lerp::lerp(0.6 * c_l_max, c_l_max, aoa_abs / deg_45) * aoa_s
-        let planform_area = std::f32::consts::PI * chord_length * span_length * 0.25;
-        Self {
-            aspect_ratio: span_length.powi(2) / planform_area,
-            planform_area,
-        }
-    }
-
-    /// The aspect ratio is the ratio of the span squared to actual planform
-    /// area
-    pub fn aspect_ratio(&self) -> f32 { self.aspect_ratio }
-
-    pub fn planform_area(&self) -> f32 { self.planform_area }
-}
-
-impl RigidWings {
-    /// Total lift coefficient for a finite wing of symmetric aerofoil shape and
-    /// elliptical pressure distribution.
-    pub fn lift_coefficient(&self, aoa: f32) -> f32 {
-        let aoa_abs = aoa.abs();
-        let stall_angle = PI * 0.1;
-        inline_tweak::tweak!(1.0)
-            * self.planform_area()
-            * if aoa_abs < stall_angle {
-                self.lift_slope(None) * aoa
-            } else if inline_tweak::tweak!(true) {
-                // This is when flow separation and turbulence starts to kick in.
-                // Going to just make something up (based on some data), as the alternative is
-                // to just throw your hands up and return 0
-                let aoa_s = aoa.signum();
-                let c_l_max = self.lift_slope(None) * stall_angle;
-                let deg_45 = PI / 4.0;
-                if aoa_abs < deg_45 {
-                    // drop directly to 0.6 * max lift at stall angle
-                    // then climb back to max at 45°
-                    Lerp::lerp(0.6 * c_l_max, c_l_max, aoa_abs / deg_45) * aoa_s
-                } else {
-                    // let's just say lift goes down linearly again until we're at 90°
-                    Lerp::lerp(c_l_max, 0.0, (aoa_abs - deg_45) / deg_45) * aoa_s
-                }
             } else {
-                0.0
+                // let's just say lift goes down linearly again until we're at 90°
+                Lerp::lerp(c_l_max, 0.0, (aoa_abs - deg_45) / deg_45) * aoa_s
             }
-    }
-
-    /// The zero-lift profile drag coefficient is the parasite drag on the wings
-    /// at the angle of attack which generates no lift
-    pub fn zero_lift_drag_coefficient(&self) -> f32 {
-        // avg value for Harris' hawk (Parabuteo unicinctus) [1]
-        self.planform_area() * 0.02
-    }
-
-    /// The change in lift over change in angle of attack¹. Multiplying by angle
-    /// of attack gives the lift coefficient (for a finite wing, not aerofoil).
-    ///
-    /// Aspect ratio is the ratio of total wing span squared over planform area.
-    ///
-    /// # Notes
-    ///
-    /// Only valid for symmetric, elliptical wings at small² angles of attack³.
-    /// Does not apply to twisted, cambered or delta wings. (It still gives a
-    /// reasonably accurate approximation if the wing shape is not truly
-    /// elliptical.)
-    ///
-    /// 1. geometric angle of attack, i.e. the pitch angle relative to
-    /// freestream flow
-    /// 2. up to around ~18°, at which point maximum lift has been achieved and
-    /// thereafter falls precipitously, causing a stall (this is the stall
-    /// angle) 3. effective aoa, i.e. geometric aoa - induced aoa; assumes
-    /// no sideslip
-    fn lift_slope(&self, sweep_angle: Option<f32>) -> f32 {
-        // lift slope for a thin aerofoil, given by Thin Aerofoil Theory
-        let ar = self.aspect_ratio();
-        let a0 = 2.0 * PI;
-        if let Some(sweep) = sweep_angle {
-            // for swept wings we use Kuchemann's modification to Helmbold's
-            // equation
-            let a0_cos_sweep = a0 * sweep.cos();
-            let x = a0_cos_sweep / (PI * ar);
-            a0_cos_sweep / ((1.0 + x.powi(2)).sqrt() + x)
-        } else if ar < 4.0 {
-            // for low aspect ratio wings (AR < 4) we use Helmbold's equation
-            let x = a0 / (PI * ar);
-            a0 / ((1.0 + x.powi(2)).sqrt() + x)
         } else {
-            // for high aspect ratio wings (AR > 4) we use the equation given by
+            0.0
-            // Prandtl's lifting-line theory
-            a0 / (1.0 + (a0 / (PI * ar)))
         }
-    }
 }
 
-/*
+/// The zero-lift profile drag coefficient is the parasite drag on the wings
-## References:
+/// at the angle of attack which generates no lift
+pub fn zero_lift_drag_coefficient(planform_area: f32) -> f32 {
+    // avg value for Harris' hawk (Parabuteo unicinctus) [1]
+    planform_area * 0.02
+}
 
-1. "Field Estimates of Body Drag Coefficient on the Basis of Dives in Passerine Birds",
+/// The change in lift over change in angle of attack¹. Multiplying by angle
-    Anders Hedenström and Felix Liechti, 2001
+/// of attack gives the lift coefficient (for a finite wing, not aerofoil).
-2. "A Simple Method to Determine Drag Coefficients in Aquatic Animals",
+/// Aspect ratio is the ratio of total wing span squared over planform area.
-    D. Bilo and W. Nachtigall, 1980
+///
-*/
+/// # Notes
+/// Only valid for symmetric, elliptical wings at small² angles of attack³.
+/// Does not apply to twisted, cambered or delta wings. (It still gives a
+/// reasonably accurate approximation if the wing shape is not truly
+/// elliptical.)
+/// 1. geometric angle of attack, i.e. the pitch angle relative to
+/// freestream flow
+/// 2. up to around ~18°, at which point maximum lift has been achieved and
+/// thereafter falls precipitously, causing a stall (this is the stall
+/// angle) 3. effective aoa, i.e. geometric aoa - induced aoa; assumes
+/// no sideslip
+fn lift_slope(aspect_ratio: f32, sweep_angle: Option<f32>) -> f32 {
+    // lift slope for a thin aerofoil, given by Thin Aerofoil Theory
+    let a0 = 2.0 * PI;
+    if let Some(sweep) = sweep_angle {
+        // for swept wings we use Kuchemann's modification to Helmbold's
+        // equation
+        let a0_cos_sweep = a0 * sweep.cos();
+        let x = a0_cos_sweep / (PI * aspect_ratio);
+        a0_cos_sweep / ((1.0 + x.powi(2)).sqrt() + x)
+    } else if aspect_ratio < 4.0 {
+        // for low aspect ratio wings (AR < 4) we use Helmbold's equation
+        let x = a0 / (PI * aspect_ratio);
+        a0 / ((1.0 + x.powi(2)).sqrt() + x)
+    } else {
+        // for high aspect ratio wings (AR > 4) we use the equation given by
+        // Prandtl's lifting-line theory
+        a0 / (1.0 + (a0 / (PI * aspect_ratio)))
+    }
+}

common/src/comp/mod.rs

@@ -68,7 +68,7 @@ pub use self::{
         InputKind, InventoryAction, InventoryEvent, InventoryManip, MountState, Mounting,
     },
     energy::{Energy, EnergyChange, EnergySource},
-    fluid_dynamics::{Fluid, RigidWings},
+    fluid_dynamics::Fluid,
     group::Group,
     home_chunk::HomeChunk,
     inputs::CanBuild,

common/src/states/glide.rs

@@ -1,6 +1,6 @@
 use super::utils::handle_climb;
 use crate::{
-    comp::{inventory::slot::EquipSlot, CharacterState, Ori, RigidWings, StateUpdate},
+    comp::{inventory::slot::EquipSlot, CharacterState, Ori, StateUpdate},
     states::behavior::{CharacterBehavior, JoinData},
     util::Dir,
 };
@@ -9,14 +9,19 @@ use vek::*;
 
 #[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
 pub struct Data {
-    pub wings: RigidWings,
+    /// The aspect ratio is the ratio of the span squared to actual planform
+    /// area
+    pub aspect_ratio: f32,
+    pub planform_area: f32,
     pub ori: Ori,
 }
 
 impl Data {
     pub fn new(span_length: f32, chord_length: f32, ori: Ori) -> Self {
+        let planform_area = std::f32::consts::PI * chord_length * span_length * 0.25;
         Self {
-            wings: RigidWings::new(span_length, chord_length),
+            aspect_ratio: span_length.powi(2) / planform_area,
+            planform_area,
             ori,
         }
     }

common/systems/src/phys.rs

@@ -9,7 +9,6 @@ use common::{
     event::{EventBus, ServerEvent},
     outcome::Outcome,
     resources::DeltaTime,
-    states,
     terrain::{Block, TerrainGrid},
     uid::Uid,
     util::{Projection, SpatialGrid},
@@ -43,7 +42,6 @@ fn fluid_density(height: f32, fluid: &Fluid) -> Density {
 fn integrate_forces(
     dt: &DeltaTime,
     mut vel: Vel,
-    ori: &Ori,
     body: &Body,
     density: &Density,
     mass: &Mass,
@@ -61,17 +59,7 @@ fn integrate_forces(
     // Aerodynamic/hydrodynamic forces
     if !rel_flow.0.is_approx_zero() {
         debug_assert!(!rel_flow.0.map(|a| a.is_nan()).reduce_or());
-        let glider: Option<&states::glide::Data> = character_state.and_then(|cs| match cs {
+        let impulse = dt.0 * body.aerodynamic_forces(&rel_flow, fluid_density.0, character_state);
-            CharacterState::Glide(data) => Some(data),
-            _ => None,
-        });
-        let impulse = dt.0
-            * body.aerodynamic_forces(
-                glider.map(|g| &g.ori).unwrap_or(ori),
-                &rel_flow,
-                fluid_density.0,
-                glider.map(|g| g.wings).as_ref(),
-            );
         debug_assert!(!impulse.map(|a| a.is_nan()).reduce_or());
         if !impulse.is_approx_zero() {
             let new_v = vel.0 + impulse / mass.0;
@@ -576,7 +564,6 @@ impl<'a> PhysicsData<'a> {
         (
             positions,
             velocities,
-            &write.orientations,
             read.stickies.maybe(),
             &read.bodies,
             read.character_states.maybe(),
@@ -595,7 +582,6 @@ impl<'a> PhysicsData<'a> {
                  (
                     pos,
                     vel,
-                    ori,
                     sticky,
                     body,
                     character_state,
@@ -627,7 +613,6 @@ impl<'a> PhysicsData<'a> {
                                 vel.0 = integrate_forces(
                                     &dt,
                                     *vel,
-                                    ori,
                                     body,
                                     density,
                                     mass,