feat(flight): enhance flight model with transonic drag and vortex lift features

src/main/java/dev/lumareth/improvedAerodynamics/flight/AeroState.java

@@ -205,10 +205,7 @@ public final class AeroState implements PhysicsBridge.Publishable {
     private void evaluate(Frame frame, Vector3dc centerOfMass, double[] v, double[] w, boolean improved, boolean commit, Wrench out) {
         out.clear();
         double weight = clock.weight();
-        double liftConstant = ImprovedAero.LIFT_CONSTANT;
-        double stallAngleScale = Config.STALL_ANGLE_SCALE.get();
-        double liftScale = Config.LIFT_SCALE.get();
-        double thinDrag = Config.STABILIZER_DRAG.get();
+        ImprovedAero.Tuning tuning = improved ? ImprovedAero.Tuning.read() : null;
         double[] r = new double[3];
         if (commit) {
             wingStalled = false;
@@ -227,7 +224,7 @@ public final class AeroState implements PhysicsBridge.Publishable {
             double input = s.type().controllable ? input(s, node, state, 0.0) : 0.0;
 
             if (improved) {
-                boolean nowStalled = ImprovedAero.evaluate(s, v, w, r, input, weight, liftConstant, liftScale, stallAngleScale, thinDrag, stalled[i], out);
+                boolean nowStalled = ImprovedAero.evaluate(s, v, w, r, input, weight, tuning, stalled[i], out);
                 if (commit) {
                     stalled[i] = nowStalled;
                     wingStalled |= nowStalled && s.type() == SurfaceType.WING;
@@ -242,7 +239,7 @@ public final class AeroState implements PhysicsBridge.Publishable {
                 probe.clear();
                 double base = 0.0;
                 if (improved) {
-                    ImprovedAero.evaluate(s, v, w, r, input, weight, liftConstant, liftScale, stallAngleScale, thinDrag, stalled[i], probe);
+                    ImprovedAero.evaluate(s, v, w, r, input, weight, tuning, stalled[i], probe);
                 } else {
                     WariumAero.evaluate(s, v, w, r, input, weight, probe);
                 }
@@ -250,7 +247,7 @@ public final class AeroState implements PhysicsBridge.Publishable {
                 probe.clear();
                 double nudged = input(s, node, state, delta);
                 if (improved) {
-                    ImprovedAero.evaluate(s, v, w, r, nudged, weight, liftConstant, liftScale, stallAngleScale, thinDrag, stalled[i], probe);
+                    ImprovedAero.evaluate(s, v, w, r, nudged, weight, tuning, stalled[i], probe);
                 } else {
                     WariumAero.evaluate(s, v, w, r, nudged, weight, probe);
                 }

src/main/java/dev/lumareth/improvedAerodynamics/flight/AeroSurface.java

@@ -13,16 +13,17 @@ package dev.lumareth.improvedAerodynamics.flight;
  * @param node            BlockPos.asLong of the Vehicle Control Node the surface listens to, {@link #NO_NODE} if none
  * @param brakeInput      air brakes only: Warium's input value, they ignore smoothing and trim
  * @param aspectRatio     span squared over area of the wing the block belongs to, 0 if not known yet
+ * @param wingGroup       whether the wing the block belongs to contains WING blocks, not only control surfaces or stabilizers
  */
 public record AeroSurface(long key, SurfaceType type, boolean slab, double x, double y, double z, int facingX, int facingZ,
                           double pressure, double fluidFactor, double speedMultiplier,
                           long node, ControlAxis axis1, double multiplier1, ControlAxis axis2, double multiplier2,
-                          double brakeInput, double aspectRatio) {
+                          double brakeInput, double aspectRatio, boolean wingGroup) {
     public static final long NO_NODE = Long.MIN_VALUE;
 
-    AeroSurface withAspectRatio(double aspectRatio) {
+    AeroSurface withWingShape(double aspectRatio, boolean wingGroup) {
         return new AeroSurface(key, type, slab, x, y, z, facingX, facingZ, pressure, fluidFactor, speedMultiplier,
-                node, axis1, multiplier1, axis2, multiplier2, brakeInput, aspectRatio);
+                node, axis1, multiplier1, axis2, multiplier2, brakeInput, aspectRatio, wingGroup);
     }
 
     /** Warium's "Direction.getClockWise(Axis.Y)" of the facing, the span direction of a horizontal surface. */

src/main/java/dev/lumareth/improvedAerodynamics/flight/ImprovedAero.java

@@ -1,5 +1,7 @@
 package dev.lumareth.improvedAerodynamics.flight;
 
+import dev.lumareth.improvedAerodynamics.Config;
+
 /**
  * Flight model based on each surface's angle of attack, used instead of {@link WariumAero} when enabled.
  * <p>
@@ -41,6 +43,18 @@ final class ImprovedAero {
      */
     static final double LIFT_CONSTANT = 160.0 / (REFERENCE_SPEED * LIFT_SLOPE * REFERENCE_ANGLE);
 
+    /** The config values the model reads, taken once per evaluation. */
+    record Tuning(double liftScale, double stallAngleScale, double thinDrag,
+                  boolean transonic, double speedOfSound, double transonicPeakDrag, double supersonicDrag,
+                  boolean vortex, double vortexStallBonus, double vortexLiftStrength) {
+        static Tuning read() {
+            return new Tuning(Config.LIFT_SCALE.get(), Config.STALL_ANGLE_SCALE.get(), Config.STABILIZER_DRAG.get(),
+                    Config.TRANSONIC_DRAG.get(), Config.SPEED_OF_SOUND.get(), Config.TRANSONIC_PEAK_DRAG.get(),
+                    Config.SUPERSONIC_DRAG.get(),
+                    Config.VORTEX_LIFT.get(), Config.VORTEX_STALL_BONUS.get(), Config.VORTEX_LIFT_STRENGTH.get());
+        }
+    }
+
     /**
      * Adds the surface's forces to out.
      *
@@ -49,12 +63,10 @@ final class ImprovedAero {
      * @param r       block center relative to the center of mass, ship axes
      * @param input   control input of the surface, -1 to 1
      * @param stalled whether the surface was stalled in the previous substep
-     * @param liftScale multiplier for the lift only, not the drag
-     * @param thinDrag parasitic drag coefficient of stabilizers and vertical wings
      * @return whether the surface is stalled now
      */
     static boolean evaluate(AeroSurface s, double[] v, double[] w, double[] r, double input, double weight,
-                            double liftConstant, double liftScale, double stallAngleScale, double thinDrag, boolean stalled, Wrench out) {
+                            Tuning tuning, boolean stalled, Wrench out) {
         // Airflow at the block: ship velocity plus w x r
         double vx = v[0] + w[1] * r[2] - w[2] * r[1];
         double vy = v[1] + w[2] * r[0] - w[0] * r[2];
@@ -71,8 +83,13 @@ final class ImprovedAero {
             return false;
         }
 
-        // Warium's own per-block drag, lower for thin surfaces since they stabilize with lift, not drag
-        double parasitic = -localSpeed * (s.type().thin ? thinDrag : s.type().parasiticDrag) * p / m;
+        // Warium's own per-block drag, lower for thin surfaces since they stabilize with lift, not drag. Shock waves
+        // multiply it near the speed of sound; this is the only zero-lift drag a craft has, induced and stall drag
+        // below are not affected
+        double shockDrag = tuning.transonic()
+                ? TransonicDrag.multiplier(localSpeed / tuning.speedOfSound(), tuning.transonicPeakDrag(), tuning.supersonicDrag())
+                : 1.0;
+        double parasitic = -localSpeed * (s.type().thin ? tuning.thinDrag() : s.type().parasiticDrag) * shockDrag * p / m;
         out.addForce(r, vx * parasitic, vy * parasitic, vz * parasitic);
 
         // Chord c along the facing, normal n up for horizontal surfaces and to the side for vertical ones
@@ -99,7 +116,13 @@ final class ImprovedAero {
         }
         alpha = Math.max(-Math.PI / 2.0, Math.min(Math.PI / 2.0, alpha));
 
-        double critical = Math.toRadians(s.slab() ? SurfaceType.SLAB_STALL_ANGLE : s.type().stallAngle) * stallAngleScale;
+        // Low aspect ratio wings keep their flow attached longer thanks to leading-edge vortices. Control surfaces
+        // in a wing's group stall with the wing instead of before it
+        double vortexFactor = tuning.vortex() && (s.type() == SurfaceType.WING || s.wingGroup() && s.type().controllable)
+                ? VortexLift.factor(s.aspectRatio()) : 0.0;
+        double baseStall = s.slab() ? SurfaceType.SLAB_STALL_ANGLE
+                : s.wingGroup() && s.type().controllable ? SurfaceType.WING.stallAngle : s.type().stallAngle;
+        double critical = Math.toRadians(baseStall) * (1.0 + tuning.vortexStallBonus() * vortexFactor) * tuning.stallAngleScale();
         double absAlpha = Math.abs(alpha);
         // Flips only where both blends agree, so the hysteresis never makes the force jump
         if (!stalled && absAlpha > critical + STALL_BLEND) {
@@ -115,11 +138,14 @@ final class ImprovedAero {
         double slope = aspectRatio > 0.0 ? 2.0 * Math.PI / (1.0 + 2.0 / aspectRatio) : LIFT_SLOPE;
         double inducedDrag = aspectRatio > 0.0 ? 1.0 / (Math.PI * EFFICIENCY * aspectRatio) : INDUCED_DRAG;
         double attachedLift = slope * alpha;
+        if (s.type() == SurfaceType.WING) {
+            attachedLift += VortexLift.liftCoefficient(alpha, vortexFactor, tuning.vortexLiftStrength());
+        }
         double lift = attached * attachedLift + (1.0 - attached) * SEPARATED_LIFT * Math.sin(2.0 * alpha);
         double drag = attached * inducedDrag * attachedLift * attachedLift
                 + (1.0 - attached) * SEPARATED_DRAG * Math.sin(alpha) * Math.sin(alpha);
 
-        double q = liftConstant * p * fd * planeSpeedSq;
+        double q = LIFT_CONSTANT * p * fd * planeSpeedSq;
         // Lift is perpendicular to the airflow in the chord/normal plane, on the side of n
         double side = u < 0.0 ? -1.0 : 1.0;
         double lx = (Math.abs(u) * nx - side * vn * cx) / planeSpeed;
@@ -129,13 +155,13 @@ final class ImprovedAero {
         double dx = -(u * cx + vn * nx) / planeSpeed;
         double dy = -(vn * ny) / planeSpeed;
         double dz = -(u * cz + vn * nz) / planeSpeed;
-        double liftForce = q * lift * liftScale;
+        double liftForce = q * lift * tuning.liftScale();
         double dragForce = q * drag / m;
         out.addForce(r, lx * liftForce + dx * dragForce, ly * liftForce + dy * dragForce, lz * liftForce + dz * dragForce);
         return stalled;
     }
 
-    private static double smoothstep(double edge0, double edge1, double x) {
+    static double smoothstep(double edge0, double edge1, double x) {
         double t = Math.max(0.0, Math.min(1.0, (x - edge0) / (edge1 - edge0)));
         return t * t * (3.0 - 2.0 * t);
     }

src/main/java/dev/lumareth/improvedAerodynamics/flight/TransonicDrag.java

@@ -0,0 +1,33 @@
+package dev.lumareth.improvedAerodynamics.flight;
+
+/**
+ * Rise of zero-lift drag around the speed of sound: shock waves build up from Mach 0.8, drag peaks at Mach 1.05 and
+ * settles to a lower, still high, level once the flow is fully supersonic from Mach 1.3.
+ */
+final class TransonicDrag {
+    private static final double DRAG_RISE_START = 0.80;
+    private static final double PEAK = 1.05;
+    private static final double SUPERSONIC_START = 1.30;
+
+    private TransonicDrag() {
+    }
+
+    /**
+     * @param mach           airspeed divided by the speed of sound
+     * @param peakDrag       multiplier at the peak
+     * @param supersonicDrag multiplier once fully supersonic
+     * @return factor for parasitic drag, 1 in subsonic flight
+     */
+    static double multiplier(double mach, double peakDrag, double supersonicDrag) {
+        if (mach <= DRAG_RISE_START) {
+            return 1.0;
+        }
+        if (mach < PEAK) {
+            return 1.0 + (peakDrag - 1.0) * ImprovedAero.smoothstep(DRAG_RISE_START, PEAK, mach);
+        }
+        if (mach < SUPERSONIC_START) {
+            return peakDrag + (supersonicDrag - peakDrag) * ImprovedAero.smoothstep(PEAK, SUPERSONIC_START, mach);
+        }
+        return supersonicDrag;
+    }
+}

src/main/java/dev/lumareth/improvedAerodynamics/flight/VortexLift.java

@@ -0,0 +1,29 @@
+package dev.lumareth.improvedAerodynamics.flight;
+
+/**
+ * Leading-edge vortices of low aspect ratio (delta) wings: the sharp edge rolls the flow up into a vortex that keeps
+ * the wing lifting far beyond the angle where a straight wing stalls.
+ */
+final class VortexLift {
+    // Polhamus' leading-edge suction analogy constant, about pi for a sharp edge
+    private static final double VORTEX_CONSTANT = Math.PI;
+    private static final double FULL_VORTEX_ASPECT_RATIO = 1.0;
+    private static final double NO_VORTEX_ASPECT_RATIO = 3.0;
+
+    private VortexLift() {
+    }
+
+    /** How much of the vortex effect a wing gets: 1 at aspect ratio 1 or less, 0 at 3 or more, 0 if its shape is unknown. */
+    static double factor(double aspectRatio) {
+        if (aspectRatio <= 0.0) {
+            return 0.0;
+        }
+        return 1.0 - ImprovedAero.smoothstep(FULL_VORTEX_ASPECT_RATIO, NO_VORTEX_ASPECT_RATIO, aspectRatio);
+    }
+
+    /** Lift coefficient added by the vortex, odd in the angle of attack like the normal lift. */
+    static double liftCoefficient(double alpha, double factor, double strength) {
+        double sin = Math.sin(alpha);
+        return strength * factor * VORTEX_CONSTANT * sin * Math.abs(sin) * Math.cos(alpha);
+    }
+}

src/main/java/dev/lumareth/improvedAerodynamics/flight/WariumCapture.java

@@ -97,7 +97,7 @@ public final class WariumCapture {
         return new AeroSurface(pos.asLong(), type, slab, pos.getX() + 0.5, pos.getY() + 0.5, pos.getZ() + 0.5,
                 facing.getStepX(), facing.getStepZ(), pressure, fluidFactor(level), speedMultiplier,
                 node, ControlAxis.of(key), ControlAxis.multiplier(key), ControlAxis.of(key2), ControlAxis.multiplier(key2),
-                brakeInput, 0.0);
+                brakeInput, 0.0, false);
     }
 
     /** Registers a node's inputs once per tick, and saves its trim to it now and then so it survives reloading. */

src/main/java/dev/lumareth/improvedAerodynamics/flight/WingShape.java

@@ -61,6 +61,7 @@ final class WingShape {
         int[] count = new int[n];
         int[] minSpan = new int[n];
         int[] maxSpan = new int[n];
+        boolean[] hasWing = new boolean[n];
         java.util.Arrays.fill(minSpan, Integer.MAX_VALUE);
         java.util.Arrays.fill(maxSpan, Integer.MIN_VALUE);
         for (int i = 0; i < n; i++) {
@@ -70,6 +71,7 @@ final class WingShape {
             int root = find(parent, i);
             int span = all[i].type().vertical ? cell[i * 3 + 1] : cell[i * 3] * Math.abs(all[i].rightX()) + cell[i * 3 + 2] * Math.abs(all[i].rightZ());
             count[root]++;
+            hasWing[root] |= all[i].type() == SurfaceType.WING;
             minSpan[root] = Math.min(minSpan[root], span);
             maxSpan[root] = Math.max(maxSpan[root], span);
         }
@@ -80,7 +82,7 @@ final class WingShape {
             int root = find(parent, i);
             double span = maxSpan[root] - minSpan[root] + 1;
             double aspectRatio = span * span / count[root];
-            all[i] = all[i].withAspectRatio(Math.max(MIN_ASPECT_RATIO, Math.min(MAX_ASPECT_RATIO, aspectRatio)));
+            all[i] = all[i].withWingShape(Math.max(MIN_ASPECT_RATIO, Math.min(MAX_ASPECT_RATIO, aspectRatio)), hasWing[root]);
         }
         return all;
     }