Refactor and abstract fluid dynamics

common/src/comp/body.rs

@@ -16,6 +16,7 @@ pub mod theropod;
 
 use crate::{
     assets::{self, Asset},
+    comp::{fluid_dynamics::*, Ori, Vel},
     consts::{HUMAN_DENSITY, WATER_DENSITY},
     make_case_elim,
     npc::NpcKind,
@@ -23,6 +24,7 @@ use crate::{
 use serde::{Deserialize, Serialize};
 use specs::{Component, DerefFlaggedStorage};
 use specs_idvs::IdvStorage;
+use std::f32::consts::PI;
 use vek::*;
 
 use super::{BuffKind, Density, Mass};
@@ -316,9 +318,9 @@ impl Body {
 
                 _ => Vec3::new(1.0, 0.75, 1.4),
             },
-            Body::BirdMedium(_) => Vec3::new(2.0, 1.0, 1.5),
-            Body::BirdLarge(_) => Vec3::new(2.0, 6.0, 3.5),
-            Body::Dragon(_) => Vec3::new(16.0, 10.0, 16.0),
+            Body::BirdMedium(body) => body.dimensions(),
+            Body::BirdLarge(body) => body.dimensions(),
+            Body::Dragon(body) => body.dimensions(),
             Body::FishMedium(_) => Vec3::new(0.5, 2.0, 0.8),
             Body::FishSmall(_) => Vec3::new(0.3, 1.2, 0.6),
             Body::Golem(_) => Vec3::new(5.0, 5.0, 7.5),
@@ -732,3 +734,121 @@ impl Body {
 impl Component for Body {
     type Storage = DerefFlaggedStorage<Self, IdvStorage<Self>>;
 }
+
+impl Body {
+    pub fn aerodynamic_forces(
+        &self,
+        ori: Option<&Ori>,
+        rel_flow: &Vel,
+        fluid_density: f32,
+    ) -> Vec3<f32> {
+        match (self, ori) {
+            (Body::BirdMedium(bird), Some(ori)) => bird_medium::FlyingBirdMedium::from((bird, ori))
+                .aerodynamic_forces(rel_flow, fluid_density),
+            (Body::BirdLarge(bird), Some(ori)) => bird_large::FlyingBirdLarge::from((bird, ori))
+                .aerodynamic_forces(rel_flow, fluid_density),
+            (Body::Dragon(not_bird), Some(ori)) => dragon::FlyingDragon::from((not_bird, ori))
+                .aerodynamic_forces(rel_flow, fluid_density),
+            _ => self.drag(rel_flow, fluid_density),
+        }
+    }
+}
+
+impl Drag for Body {
+    /// Parasite drag is the sum of pressure drag and skin friction.
+    /// Skin friction is the drag arising from the shear forces between a fluid
+    /// and a surface, while pressure drag is due to flow separation. Both are
+    /// viscous effects.
+    ///
+    /// This coefficient includes the reference area.
+    fn parasite_drag_coefficient(&self) -> f32 {
+        // Reference area and drag coefficient assumes best-case scenario of the
+        // orientation producing least amount of drag
+        match self {
+            // Cross-section, head/feet first
+            Body::BipedLarge(_) | Body::BipedSmall(_) | Body::Golem(_) | Body::Humanoid(_) => {
+                const SCALE: f32 = 0.7;
+                let radius = self.dimensions().xy().map(|a| SCALE * a * 0.5);
+                const CD: f32 = 0.7;
+                CD * PI * radius.x * radius.y
+            },
+
+            // Cross-section, nose/tail first
+            Body::Theropod(_)
+            | Body::QuadrupedMedium(_)
+            | Body::QuadrupedSmall(_)
+            | Body::QuadrupedLow(_) => {
+                let radius = self.dimensions().map(|a| a * 0.5);
+                let cd: f32 = if matches!(self, Body::QuadrupedLow(_)) {
+                    0.7
+                } else {
+                    1.0
+                };
+                cd * PI * radius.x * radius.z
+            },
+
+            Body::BirdMedium(bird) => bird.parasite_drag_coefficient(),
+            Body::BirdLarge(bird) => bird.parasite_drag_coefficient(),
+            Body::Dragon(not_bird) => not_bird.parasite_drag_coefficient(),
+
+            // Cross-section, zero-lift angle; exclude the fins
+            Body::FishMedium(_) | Body::FishSmall(_) => {
+                let radius = self.dimensions().map(|a| a * 0.5);
+                // "A Simple Method to Determine Drag Coefficients in Aquatic Animals",
+                // D. Bilo and W. Nachtigall, 1980
+                const CD: f32 = 0.031;
+                CD * PI * radius.x * radius.z
+            },
+
+            Body::Object(object) => match object {
+                // very streamlined objects
+                object::Body::Arrow
+                | object::Body::ArrowSnake
+                | object::Body::ArrowTurret
+                | object::Body::FireworkBlue
+                | object::Body::FireworkGreen
+                | object::Body::FireworkPurple
+                | object::Body::FireworkRed
+                | object::Body::FireworkWhite
+                | object::Body::FireworkYellow
+                | object::Body::MultiArrow => {
+                    let radius = self.dimensions().map(|a| a * 0.5);
+                    const CD: f32 = 0.02;
+                    CD * PI * radius.x * radius.z
+                },
+
+                // spherical-ish objects
+                object::Body::BoltFire
+                | object::Body::BoltFireBig
+                | object::Body::BoltNature
+                | object::Body::Bomb
+                | object::Body::PotionBlue
+                | object::Body::PotionGreen
+                | object::Body::PotionRed
+                | object::Body::Pouch
+                | object::Body::Pumpkin
+                | object::Body::Pumpkin2
+                | object::Body::Pumpkin3
+                | object::Body::Pumpkin4
+                | object::Body::Pumpkin5 => {
+                    let radius = self.dimensions().xy().map(|a| a * 0.5);
+                    const CD: f32 = 0.5;
+                    CD * PI * radius.product()
+                },
+
+                _ => {
+                    let dim = self.dimensions();
+                    const CD: f32 = 2.0;
+                    CD * (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0)
+                },
+            },
+
+            Body::Ship(_) => {
+                // Airships tend to use the square of the cube root of its volume for
+                // reference area
+                let dim = self.dimensions();
+                (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0)
+            },
+        }
+    }
+}

common/src/comp/body/bird_large.rs

@@ -1,6 +1,13 @@
-use crate::{make_case_elim, make_proj_elim};
+use crate::{
+    comp::{
+        fluid_dynamics::{Drag, WingShape, WingState, Glide},
+        Ori,
+    },
+    make_case_elim, make_proj_elim,
+};
 use rand::{seq::SliceRandom, thread_rng};
 use serde::{Deserialize, Serialize};
+use vek::*;
 
 make_proj_elim!(
     body,
@@ -23,6 +30,18 @@ impl Body {
         let body_type = *(&ALL_BODY_TYPES).choose(rng).unwrap();
         Self { species, body_type }
     }
+
+    /// Dimensions of the body (wings folded)
+    pub const fn dimensions(&self) -> Vec3<f32> { Vec3::new(1.0, 5.0, 2.4) }
+
+    /// Distance from wing tip to wing tip and leading edge to trailing edge
+    /// respectively
+    // TODO: Check
+    pub const fn wing_dimensions(&self) -> Vec2<f32> { Vec2::new(7.0, 1.0) }
+
+    pub fn flying<'a>(&'a self, ori: &'a Ori) -> FlyingBirdLarge<'a> {
+        FlyingBirdLarge::from((self, ori))
+    }
 }
 
 impl From<Body> for super::Body {
@@ -83,3 +102,55 @@ make_case_elim!(
 );
 
 pub const ALL_BODY_TYPES: [BodyType; 2] = [BodyType::Female, BodyType::Male];
+
+#[derive(Copy, Clone)]
+pub struct FlyingBirdLarge<'a> {
+    wing_shape: WingShape,
+    wing_state: WingState,
+    planform_area: f32,
+    body: &'a Body,
+    ori: &'a Ori,
+}
+
+impl<'a> From<(&'a Body, &'a Ori)> for FlyingBirdLarge<'a> {
+    fn from((body, ori): (&'a Body, &'a Ori)) -> Self {
+        let Vec2 {
+            x: span_length,
+            y: chord_length,
+        } = body.wing_dimensions();
+        let planform_area = WingShape::elliptical_planform_area(span_length, chord_length);
+        FlyingBirdLarge {
+            wing_shape: WingShape::elliptical(span_length, chord_length),
+            wing_state: WingState::Flapping,
+            planform_area,
+            body,
+            ori,
+        }
+    }
+}
+
+impl Drag for Body {
+    fn parasite_drag_coefficient(&self) -> f32 {
+        let radius = self.dimensions().map(|a| a * 0.5);
+        // "Field Estimates of body::Body Drag Coefficient on the Basis of
+        // Dives in Passerine Birds", Anders Hedenström and Felix Liechti, 2001
+        const CD: f32 = 0.2;
+        CD * std::f32::consts::PI * radius.x * radius.z
+    }
+}
+
+impl Drag for FlyingBirdLarge<'_> {
+    fn parasite_drag_coefficient(&self) -> f32 {
+        self.body.parasite_drag_coefficient() + self.planform_area * 0.004
+    }
+}
+
+impl Glide for FlyingBirdLarge<'_> {
+    fn wing_shape(&self) -> &WingShape { &self.wing_shape }
+
+    fn is_gliding(&self) -> bool { matches!(self.wing_state, WingState::Fixed) }
+
+    fn planform_area(&self) -> f32 { self.planform_area }
+
+    fn ori(&self) -> &Ori { self.ori }
+}

common/src/comp/body/bird_medium.rs

@@ -1,6 +1,13 @@
-use crate::{make_case_elim, make_proj_elim};
+use crate::{
+    comp::{
+        fluid_dynamics::{Drag, WingShape, WingState, Glide},
+        Ori,
+    },
+    make_case_elim, make_proj_elim,
+};
 use rand::{seq::SliceRandom, thread_rng};
 use serde::{Deserialize, Serialize};
+use vek::*;
 
 make_proj_elim!(
     body,
@@ -23,6 +30,18 @@ impl Body {
         let body_type = *(&ALL_BODY_TYPES).choose(rng).unwrap();
         Self { species, body_type }
     }
+
+    /// Dimensions of the body (wings folded)
+    pub const fn dimensions(&self) -> Vec3<f32> { Vec3::new(0.5, 1.0, 1.1) }
+
+    /// Distance from wing tip to wing tip and leading edge to trailing edge
+    /// respectively
+    // TODO: Check
+    pub const fn wing_dimensions(&self) -> Vec2<f32> { Vec2::new(2.0, 0.4) }
+
+    pub fn flying<'a>(&'a self, ori: &'a Ori) -> FlyingBirdMedium<'a> {
+        FlyingBirdMedium::from((self, ori))
+    }
 }
 
 impl From<Body> for super::Body {
@@ -102,3 +121,55 @@ make_case_elim!(
     }
 );
 pub const ALL_BODY_TYPES: [BodyType; 2] = [BodyType::Female, BodyType::Male];
+
+#[derive(Copy, Clone)]
+pub struct FlyingBirdMedium<'a> {
+    wing_shape: WingShape,
+    wing_state: WingState,
+    planform_area: f32,
+    body: &'a Body,
+    ori: &'a Ori,
+}
+
+impl<'a> From<(&'a Body, &'a Ori)> for FlyingBirdMedium<'a> {
+    fn from((body, ori): (&'a Body, &'a Ori)) -> Self {
+        let Vec2 {
+            x: span_length,
+            y: chord_length,
+        } = body.wing_dimensions();
+        let planform_area = WingShape::elliptical_planform_area(span_length, chord_length);
+        FlyingBirdMedium {
+            wing_shape: WingShape::elliptical(span_length, chord_length),
+            wing_state: WingState::Flapping,
+            planform_area,
+            body,
+            ori,
+        }
+    }
+}
+
+impl Drag for Body {
+    fn parasite_drag_coefficient(&self) -> f32 {
+        let radius = self.dimensions().map(|a| a * 0.5);
+        // "Field Estimates of body::Body Drag Coefficient on the Basis of
+        // Dives in Passerine Birds", Anders Hedenström and Felix Liechti, 2001
+        const CD: f32 = 0.2;
+        CD * std::f32::consts::PI * radius.x * radius.z
+    }
+}
+
+impl Drag for FlyingBirdMedium<'_> {
+    fn parasite_drag_coefficient(&self) -> f32 {
+        self.body.parasite_drag_coefficient() + self.planform_area * 0.004
+    }
+}
+
+impl Glide for FlyingBirdMedium<'_> {
+    fn wing_shape(&self) -> &WingShape { &self.wing_shape }
+
+    fn is_gliding(&self) -> bool { matches!(self.wing_state, WingState::Fixed) }
+
+    fn planform_area(&self) -> f32 { self.planform_area }
+
+    fn ori(&self) -> &Ori { self.ori }
+}

common/src/comp/body/dragon.rs

@@ -1,6 +1,13 @@
-use crate::{make_case_elim, make_proj_elim};
+use crate::{
+    comp::{
+        fluid_dynamics::{Drag, Glide, WingShape, WingState},
+        Ori,
+    },
+    make_case_elim, make_proj_elim,
+};
 use rand::{seq::SliceRandom, thread_rng};
 use serde::{Deserialize, Serialize};
+use vek::*;
 
 make_proj_elim!(
     body,
@@ -23,6 +30,18 @@ impl Body {
         let body_type = *(&ALL_BODY_TYPES).choose(rng).unwrap();
         Self { species, body_type }
     }
+
+    /// Dimensions of the body (wings folded)
+    pub const fn dimensions(&self) -> Vec3<f32> { Vec3::new(5.0, 10.0, 16.0) }
+
+    /// Distance from wing tip to wing tip and leading edge to trailing edge
+    /// respectively
+    // TODO: Check
+    pub const fn wing_dimensions(&self) -> Vec2<f32> { Vec2::new(16.0, 5.0) }
+
+    pub fn flying<'a>(&'a self, ori: &'a Ori) -> FlyingDragon<'a> {
+        FlyingDragon::from((self, ori))
+    }
 }
 
 impl From<Body> for super::Body {
@@ -77,3 +96,55 @@ make_case_elim!(
 );
 
 pub const ALL_BODY_TYPES: [BodyType; 2] = [BodyType::Female, BodyType::Male];
+
+#[derive(Copy, Clone)]
+pub struct FlyingDragon<'a> {
+    wing_shape: WingShape,
+    wing_state: WingState,
+    planform_area: f32,
+    body: &'a Body,
+    ori: &'a Ori,
+}
+
+impl<'a> From<(&'a Body, &'a Ori)> for FlyingDragon<'a> {
+    fn from((body, ori): (&'a Body, &'a Ori)) -> Self {
+        let Vec2 {
+            x: span_length,
+            y: chord_length,
+        } = body.wing_dimensions();
+        let planform_area = WingShape::elliptical_planform_area(span_length, chord_length);
+        FlyingDragon {
+            wing_shape: WingShape::elliptical(span_length, chord_length),
+            wing_state: WingState::Flapping,
+            planform_area,
+            body,
+            ori,
+        }
+    }
+}
+
+impl Drag for Body {
+    fn parasite_drag_coefficient(&self) -> f32 {
+        let radius = self.dimensions().map(|a| a * 0.5);
+        // "Field Estimates of body::Body Drag Coefficient on the Basis of
+        // Dives in Passerine Birds", Anders Hedenström and Felix Liechti, 2001
+        const CD: f32 = 0.2;
+        CD * std::f32::consts::PI * radius.x * radius.z
+    }
+}
+
+impl Drag for FlyingDragon<'_> {
+    fn parasite_drag_coefficient(&self) -> f32 {
+        self.body.parasite_drag_coefficient() + self.planform_area * 0.004
+    }
+}
+
+impl Glide for FlyingDragon<'_> {
+    fn wing_shape(&self) -> &WingShape { &self.wing_shape }
+
+    fn is_gliding(&self) -> bool { matches!(self.wing_state, WingState::Fixed) }
+
+    fn planform_area(&self) -> f32 { self.planform_area }
+
+    fn ori(&self) -> &Ori { self.ori }
+}

common/src/comp/fluid_dynamics.rs

@@ -1,5 +1,3 @@
-#[cfg(not(target_arch = "wasm32"))]
-use super::body::{object, Body};
 use super::{Density, Ori, Vel};
 use crate::{
     consts::{AIR_DENSITY, LAVA_DENSITY, WATER_DENSITY},
@@ -120,20 +118,12 @@ impl Default for Fluid {
     }
 }
 
-pub struct Wings {
-    pub aspect_ratio: f32,
-    pub planform_area: f32,
-    pub ori: Ori,
-}
+pub trait Drag: Clone {
+    /// Drag coefficient for skin friction and flow separation.
+    /// (Multiplied by reference area.)
+    fn parasite_drag_coefficient(&self) -> f32;
 
-#[cfg(not(target_arch = "wasm32"))]
-impl Body {
-    pub fn aerodynamic_forces(
-        &self,
-        rel_flow: &Vel,
-        fluid_density: f32,
-        wings: Option<&Wings>,
-    ) -> Vec3<f32> {
+    fn drag(&self, rel_flow: &Vel, fluid_density: f32) -> Vec3<f32> {
         let v_sq = rel_flow.0.magnitude_squared();
         if v_sq < 0.25 {
             // don't bother with miniscule forces
@@ -141,180 +131,99 @@ impl Body {
         } else {
             let rel_flow_dir = Dir::new(rel_flow.0 / v_sq.sqrt());
             // All the coefficients come pre-multiplied by their reference area
-            0.5 * fluid_density
-                * v_sq
-                * match wings {
-                    Some(&Wings {
-                        aspect_ratio,
-                        planform_area,
-                        ori,
-                    }) => {
-                        if aspect_ratio > 25.0 {
-                            tracing::warn!(
-                                "Calculating lift for wings with an aspect ratio of {}. The \
-                                 formulas are only valid for aspect ratios below 25.",
-                                aspect_ratio
-                            )
-                        };
-                        let ar = aspect_ratio.min(24.0);
-                        // We have an elliptical wing; proceed to calculate its lift and drag
+            0.5 * fluid_density * v_sq * self.parasite_drag_coefficient() * *rel_flow_dir
+        }
+    }
+}
 
+pub trait Glide: Drag {
+    const STALL_ANGLE: f32 = PI * 0.1;
+
+    fn wing_shape(&self) -> &WingShape;
+
+    fn planform_area(&self) -> f32;
+
+    fn ori(&self) -> &Ori;
+
+    fn is_gliding(&self) -> bool;
+
+    /// Total lift coefficient for a finite wing of symmetric aerofoil shape and
+    /// elliptical pressure distribution. (Multiplied by reference area.)
+    fn lift_coefficient(&self, angle_of_attack: f32) -> f32 {
+        let aoa_abs = angle_of_attack.abs();
+        self.planform_area()
+            * if self.is_gliding() { 1.0 } else { 0.2 }
+            * if aoa_abs <= Self::STALL_ANGLE {
+                self.wing_shape().lift_slope() * angle_of_attack
+            } else {
+                // 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 = angle_of_attack.signum();
+                let c_l_max = self.wing_shape().lift_slope() * Self::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
+                }
+            }
+    }
+
+    /// Total drag coefficient (multiplied by reference area)
+    fn drag_coefficient(&self, aspect_ratio: f32, lift_coefficient: f32) -> f32 {
+        self.parasite_drag_coefficient() + induced_drag_coefficient(aspect_ratio, lift_coefficient)
+    }
+
+    fn aerodynamic_forces(&self, rel_flow: &Vel, fluid_density: f32) -> Vec3<f32> {
+        let v_sq = rel_flow.0.magnitude_squared();
+        if v_sq < 0.1 {
+            // don't bother with miniscule forces
+            Vec3::zero()
+        } else {
+            let q = 0.5 * fluid_density * v_sq;
+            let rel_flow_dir = Dir::new(rel_flow.0 / v_sq.sqrt());
+            if self.is_gliding() {
+                let ori = self.ori();
+                let ar = self.wing_shape().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 = lift_coefficient(ar, planform_area, aoa);
+                let c_l = self.lift_coefficient(aoa);
 
-                        // lift dir will be orthogonal to the local relative flow vector.
-                        // Local relative flow is the resulting vector of (relative) freestream
-                        // 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
-                            let aoa_eff = aoa - aoa_i;
-                            // Angle between chord line and local relative wind is aoa_eff
-                            // radians. 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()
-                        };
+                let lift = q * c_l * *lift_dir(ori, ar, c_l, aoa);
+                let drag = q * self.drag_coefficient(ar, c_l) * *rel_flow_dir;
 
-                        // drag coefficient
-                        let c_d = {
-                            // Oswald's efficiency factor (empirically derived--very magical)
-                            // (this definition should not be used for aspect ratios > 25)
-                            let e = 1.78 * (1.0 - 0.045 * ar.powf(0.68)) - 0.64;
-                            // induced drag coefficient (drag due to lift)
-                            let cdi = c_l.powi(2) / (PI * e * ar);
-
-                            zero_lift_drag_coefficient(planform_area)
-                                + self.parasite_drag_coefficient()
-                                + cdi
-                        };
-                        debug_assert!(c_d.is_sign_positive());
-                        debug_assert!(c_l.is_sign_positive() || aoa.is_sign_negative());
-
-                        c_l * *lift_dir + c_d * *rel_flow_dir
-                    },
-
-                    _ => self.parasite_drag_coefficient() * *rel_flow_dir,
-                }
-        }
-    }
-
-    /// Parasite drag is the sum of pressure drag and skin friction.
-    /// Skin friction is the drag arising from the shear forces between a fluid
-    /// and a surface, while pressure drag is due to flow separation. Both are
-    /// viscous effects.
-    ///
-    /// This coefficient includes the reference area.
-    fn parasite_drag_coefficient(&self) -> f32 {
-        // Reference area and drag coefficient assumes best-case scenario of the
-        // orientation producing least amount of drag
-        match self {
-            // Cross-section, head/feet first
-            Body::BipedLarge(_) | Body::BipedSmall(_) | Body::Golem(_) | Body::Humanoid(_) => {
-                let dim = self.dimensions().xy().map(|a| 0.7 * a * 0.5);
-                const CD: f32 = 0.7;
-                CD * PI * dim.x * dim.y
-            },
-
-            // Cross-section, nose/tail first
-            Body::Theropod(_)
-            | Body::QuadrupedMedium(_)
-            | Body::QuadrupedSmall(_)
-            | Body::QuadrupedLow(_) => {
-                let dim = self.dimensions().map(|a| a * 0.5);
-                let cd: f32 = if matches!(self, Body::QuadrupedLow(_)) {
-                    0.7
+                lift + drag
             } else {
-                    1.0
-                };
-                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);
-                let cd: f32 = match self {
-                    // "Field Estimates of Body Drag Coefficient on the Basis of Dives in Passerine
-                    // Birds", Anders Hedenström and Felix Liechti, 2001
-                    Body::BirdLarge(_) | Body::BirdMedium(_) => 0.2,
-                    // arbitrary
-                    _ => 0.7,
-                };
-                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);
-                // "A Simple Method to Determine Drag Coefficients in Aquatic Animals",
-                // D. Bilo and W. Nachtigall, 1980
-                const CD: f32 = 0.031;
-                CD * PI * dim.x * 0.2 * dim.z
-            },
-
-            Body::Object(object) => match object {
-                // very streamlined objects
-                object::Body::Arrow
-                | object::Body::ArrowSnake
-                | object::Body::ArrowTurret
-                | object::Body::FireworkBlue
-                | object::Body::FireworkGreen
-                | object::Body::FireworkPurple
-                | object::Body::FireworkRed
-                | object::Body::FireworkWhite
-                | object::Body::FireworkYellow
-                | object::Body::MultiArrow => {
-                    let dim = self.dimensions().map(|a| a * 0.5);
-                    const CD: f32 = 0.02;
-                    CD * PI * dim.x * dim.z
-                },
-
-                // spherical-ish objects
-                object::Body::BoltFire
-                | object::Body::BoltFireBig
-                | object::Body::BoltNature
-                | object::Body::Bomb
-                | object::Body::PotionBlue
-                | object::Body::PotionGreen
-                | object::Body::PotionRed
-                | object::Body::Pouch
-                | object::Body::Pumpkin
-                | object::Body::Pumpkin2
-                | object::Body::Pumpkin3
-                | object::Body::Pumpkin4
-                | object::Body::Pumpkin5 => {
-                    let dim = self.dimensions().map(|a| a * 0.5);
-                    const CD: f32 = 0.5;
-                    CD * PI * dim.x * dim.z
-                },
-
-                _ => {
-                    let dim = self.dimensions();
-                    const CD: f32 = 2.0;
-                    CD * (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0)
-                },
-            },
-
-            Body::Ship(_) => {
-                // Airships tend to use the square of the cube root of its volume for
-                // reference area
-                let dim = self.dimensions();
-                (PI / 6.0 * dim.x * dim.y * dim.z).powf(2.0 / 3.0)
-            },
+                q * self.parasite_drag_coefficient() * *rel_flow_dir
             }
         }
     }
+}
+
+
+pub fn lift_dir(ori: &Ori, aspect_ratio: f32, lift_coefficient: f32, angle_of_attack: f32) -> Dir {
+    // lift dir will be orthogonal to the local relative flow vector. Local relative
+    // flow is the resulting vector of (relative) freestream flow + downwash
+    // (created by the vortices of the wing tips)
+
+    // induced angle of attack
+    let aoa_i = lift_coefficient / (PI * aspect_ratio);
+    // effective angle of attack; the aoa as seen by aerofoil after downwash
+    let aoa_eff = angle_of_attack - aoa_i;
+
+    // Angle between chord line and local relative wind is aoa_eff radians.
+    // 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()
+}
 
 /// Geometric angle of attack
 ///
@@ -329,35 +238,39 @@ pub fn angle_of_attack(ori: &Ori, rel_flow_dir: &Dir) -> f32 {
         .unwrap_or(0.0)
 }
 
-/// Total lift coefficient for a finite wing of symmetric aerofoil shape and
-/// elliptical pressure distribution.
-pub fn lift_coefficient(aspect_ratio: f32, planform_area: f32, aoa: f32) -> f32 {
-    let aoa_abs = aoa.abs();
-    let stall_angle = PI * 0.1;
-    planform_area
-        * if aoa_abs < stall_angle {
-            lift_slope(aspect_ratio, None) * aoa
-        } else {
-            // 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 = lift_slope(aspect_ratio, 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
+#[derive(Copy, Clone)]
+pub enum WingState {
+    Fixed,
+    Flapping,
+    Retracted,
 }
+
+#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub enum WingShape {
+    Elliptical { aspect_ratio: f32 },
+    // Tapered { aspect_ratio: f32, e: f32 },
+    Swept { aspect_ratio: f32, angle: f32 },
+    // Delta,
+}
+
+impl WingShape {
+    pub fn aspect_ratio(&self) -> f32 {
+        match self {
+            Self::Elliptical { aspect_ratio } => *aspect_ratio,
+            Self::Swept { aspect_ratio, .. } => *aspect_ratio,
         }
     }
 
-/// 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(planform_area: f32) -> f32 { planform_area * 0.004 }
+    pub fn elliptical_planform_area(span_length: f32, chord_length: f32) -> f32 {
+        std::f32::consts::PI * span_length * chord_length * 0.25
+    }
+
+    pub fn elliptical(span_length: f32, chord_length: f32) -> Self {
+        Self::Elliptical {
+            aspect_ratio: span_length.powi(2)
+                / Self::elliptical_planform_area(span_length, chord_length),
+        }
+    }
 
     /// 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).
@@ -374,17 +287,12 @@ pub fn zero_lift_drag_coefficient(planform_area: f32) -> f32 { planform_area * 0
     /// thereafter falls precipitously, causing a stall (this is the stall
     /// angle) 3. effective aoa, i.e. geometric aoa - induced aoa; assumes
     /// no sideslip
-// TODO: Look into handling tapered wings
-fn lift_slope(aspect_ratio: f32, sweep_angle: Option<f32>) -> f32 {
+    pub fn lift_slope(&self) -> 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 {
+        match self {
+            WingShape::Elliptical { aspect_ratio } => {
+                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)
@@ -393,4 +301,36 @@ fn lift_slope(aspect_ratio: f32, sweep_angle: Option<f32>) -> f32 {
                     // Prandtl's lifting-line theory
                     a0 / (1.0 + (a0 / (PI * aspect_ratio)))
                 }
+            },
+            // WingShape::Tapered { aspect_ratio, e } => todo!(),
+            WingShape::Swept {
+                aspect_ratio,
+                angle,
+            } => {
+                // for swept wings we use Kuchemann's modification to Helmbold's
+                // equation
+                let a0_cos_sweep = a0 * angle.cos();
+                let x = a0_cos_sweep / (PI * aspect_ratio);
+                a0_cos_sweep / ((1.0 + x.powi(2)).sqrt() + x)
+            },
+        }
+    }
+}
+
+/// Induced drag coefficient (drag due to lift)
+pub fn induced_drag_coefficient(aspect_ratio: f32, lift_coefficient: f32) -> f32 {
+    let ar = aspect_ratio;
+    if ar > 25.0 {
+        tracing::warn!(
+            "Calculating induced drag for wings with a given aspect ratio of {}. The formulas are \
+             only valid for aspect ratios below 25, so it will be substituted.",
+            ar
+        )
+    };
+    let ar = ar.min(24.0);
+    // Oswald's efficiency factor (empirically derived--very magical)
+    // (this definition should not be used for aspect ratios > 25)
+    let e = 1.78 * (1.0 - 0.045 * ar.powf(0.68)) - 0.64;
+    // induced drag coefficient (drag due to lift)
+    lift_coefficient.powi(2) / (PI * e * ar)
 }

common/src/consts.rs

@@ -2,7 +2,7 @@
 pub const MAX_PICKUP_RANGE: f32 = 8.0;
 pub const MAX_MOUNT_RANGE: f32 = 14.0;
 
-pub const GRAVITY: f32 = 9.8;
+pub const GRAVITY: f32 = 25.0;
 pub const FRIC_GROUND: f32 = 0.15;
 
 // Values for air taken from http://www-mdp.eng.cam.ac.uk/web/library/enginfo/aerothermal_dvd_only/aero/atmos/atmos.html

common/src/lib.rs

@@ -6,6 +6,7 @@
 #![feature(
     arbitrary_enum_discriminant,
     associated_type_defaults,
+    bindings_after_at,
     bool_to_option,
     const_generics,
     fundamental,

common/src/states/glide.rs

@@ -1,8 +1,9 @@
 use super::utils::handle_climb;
 use crate::{
     comp::{
-        fluid_dynamics::angle_of_attack, inventory::slot::EquipSlot, CharacterState, Ori,
-        StateUpdate, Vel,
+        fluid_dynamics::{angle_of_attack, Drag, Glide, WingShape},
+        inventory::slot::EquipSlot,
+        CharacterState, Ori, StateUpdate, Vel,
     },
     states::behavior::{CharacterBehavior, JoinData},
     util::{Dir, Plane, Projection},
@@ -17,7 +18,7 @@ const PITCH_SLOW_TIME: f32 = 0.5;
 pub struct Data {
     /// The aspect ratio is the ratio of the span squared to actual planform
     /// area
-    pub aspect_ratio: f32,
+    pub wing_shape: WingShape,
     pub planform_area: f32,
     pub ori: Ori,
     last_vel: Vel,
@@ -36,7 +37,7 @@ impl Data {
         let planform_area = PI * chord_length * span_length * 0.25;
         let aspect_ratio = span_length.powi(2) / planform_area;
         Self {
-            aspect_ratio,
+            wing_shape: WingShape::Elliptical { aspect_ratio },
             planform_area,
             ori,
             last_vel: Vel::zero(),
@@ -69,6 +70,20 @@ impl Data {
     }
 }
 
+impl Drag for Data {
+    fn parasite_drag_coefficient(&self) -> f32 { self.planform_area * 0.004 }
+}
+
+impl Glide for Data {
+    fn wing_shape(&self) -> &WingShape { &self.wing_shape }
+
+    fn is_gliding(&self) -> bool { true }
+
+    fn planform_area(&self) -> f32 { self.planform_area }
+
+    fn ori(&self) -> &Ori { &self.ori }
+}
+
 impl CharacterBehavior for Data {
     fn behavior(&self, data: &JoinData) -> StateUpdate {
         let mut update = StateUpdate::from(data);

common/systems/src/phys.rs

@@ -1,7 +1,7 @@
 use common::{
     comp::{
         body::ship::figuredata::{VoxelCollider, VOXEL_COLLIDER_MANIFEST},
-        fluid_dynamics::{Fluid, LiquidKind, Wings},
+        fluid_dynamics::{Fluid, LiquidKind, Glide},
         BeamSegment, Body, CharacterState, Collider, Density, Mass, Mounting, Ori, PhysicsState,
         Pos, PosVelDefer, PreviousPhysCache, Projectile, Scale, Shockwave, Stats, Sticky, Vel,
     },
@@ -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,9 @@ fn fluid_density(height: f32, fluid: &Fluid) -> Density {
 fn integrate_forces(
     dt: &DeltaTime,
     mut vel: Vel,
-    (body, wings): (&Body, Option<&Wings>),
+    ori: Option<&Ori>,
+    body: &Body,
+    character_state: Option<&CharacterState>,
     density: &Density,
     mass: &Mass,
     fluid: &Fluid,
@@ -59,10 +60,18 @@ 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 impulse = dt.0 * body.aerodynamic_forces(&rel_flow, fluid_density.0, wings);
-        debug_assert!(!impulse.map(|a| a.is_nan()).reduce_or());
-        if !impulse.is_approx_zero() {
-            let new_v = vel.0 + impulse / mass.0;
+        let aerodynamic_forces = body.aerodynamic_forces(ori, &rel_flow, fluid_density.0)
+            + match character_state {
+                Some(CharacterState::Glide(glider)) => {
+                    glider.aerodynamic_forces(&rel_flow, fluid_density.0)
+                },
+                _ => Vec3::zero(),
+            };
+        // let impulse = dt.0 * body.aerodynamic_forces(ori, &rel_flow,
+        // fluid_density.0);
+        debug_assert!(!aerodynamic_forces.map(|a| a.is_nan()).reduce_or());
+        if !aerodynamic_forces.is_approx_zero() {
+            let new_v = vel.0 + dt.0 * aerodynamic_forces / mass.0;
             // If the new velocity is in the opposite direction, it's because the forces
             // involved are too high for the current tick to handle. We deal with this by
             // removing the component of our velocity vector along the direction of force.
@@ -71,7 +80,7 @@ fn integrate_forces(
             if new_v.dot(vel.0) < 0.0 {
                 // Multiply by a factor to prevent full stop, as this can cause things to get
                 // stuck in high-density medium
-                vel.0 -= vel.0.projected(&impulse) * 0.9;
+                vel.0 -= vel.0.projected(&aerodynamic_forces) * 0.9;
             } else {
                 vel.0 = new_v;
             }
@@ -81,8 +90,7 @@ fn integrate_forces(
 
     // Hydrostatic/aerostatic forces
     // modify gravity to account for the effective density as a result of buoyancy
-    let down_force = dt.0 * gravity * (density.0 - fluid_density.0) / density.0;
-    vel.0.z -= down_force;
+    vel.0.z -= dt.0 * gravity * (density.0 - fluid_density.0) / density.0;
 
     vel
 }
@@ -576,6 +584,7 @@ impl<'a> PhysicsData<'a> {
             read.character_states.maybe(),
             &write.physics_states,
             &read.masses,
+            write.orientations.maybe(),
             &read.densities,
             !&read.mountings,
         )
@@ -594,6 +603,7 @@ impl<'a> PhysicsData<'a> {
                     character_state,
                     physics_state,
                     mass,
+                    ori,
                     density,
                     _,
                 )| {
@@ -617,24 +627,12 @@ impl<'a> PhysicsData<'a> {
                                 vel.0.z -= dt.0 * GRAVITY;
                             },
                             Some(fluid) => {
-                                let wings = match character_state {
-                                    Some(&CharacterState::Glide(states::glide::Data {
-                                        aspect_ratio,
-                                        planform_area,
-                                        ori,
-                                        ..
-                                    })) => Some(Wings {
-                                        aspect_ratio,
-                                        planform_area,
-                                        ori,
-                                    }),
-
-                                    _ => None,
-                                };
                                 vel.0 = integrate_forces(
                                     &dt,
                                     *vel,
-                                    (body, wings.as_ref()),
+                                    ori,
+                                    body,
+                                    character_state,
                                     density,
                                     mass,
                                     &fluid,