Improve coverage and fix direction

Author: 1-Bart-1

src/wing_geometry.jl

@@ -1213,44 +1213,22 @@ function refine_mesh_by_splitting_provided_sections!(
                 reuse_aero_data
             )
 
-            # Apply catenary billowing to the just-created sections
+            # Apply billowing by rotating chords around LE
             if billowing_percentage > 0 && idx > start_idx
-                LE_1 = LE[left_section_index]
-                TE_1 = TE[left_section_index]
-                LE_2 = LE[left_section_index + 1]
-                TE_2 = TE[left_section_index + 1]
-
-                chord_1 = TE_1 - LE_1
-                chord_2 = TE_2 - LE_2
-                x_hat = normalize((chord_1 + chord_2) / 2)
-                y_hat = normalize(LE_1 - LE_2)
-                z_hat = cross(x_hat, y_hat)
-                z_norm = norm(z_hat)
-                if z_norm > 1e-10
-                    z_hat = z_hat / z_norm
-
-                    d = abs(dot(TE_1 - TE_2, y_hat))
-
-                    if d > 1e-12
-                        a = catenary_parameter(
-                            billowing_percentage, d)
-                        u = d / (2a)
-                        span_len = norm(LE_1 - LE_2)
-                        TE_mid = (TE_1 + TE_2) / 2
-
-                        for si in start_idx:(idx - 1)
-                            sec = wing.refined_sections[si]
-                            t = dot(sec.LE_point - LE_2,
-                                    y_hat) / span_len
-                            x = d * (t - 0.5)
-                            sag = a * (cosh(x / a) -
-                                       cosh(u))
-                            arc_y = d * (t - 0.5)
-                            sec.TE_point .= TE_mid +
-                                arc_y * y_hat + sag * z_hat
-                        end
-                    end
-                end
+                y_hat = normalize(
+                    LE[left_section_index] -
+                    LE[left_section_index + 1])
+                span_len = norm(
+                    LE[left_section_index] -
+                    LE[left_section_index + 1])
+                apply_billowing_to_pair!(
+                    wing.refined_sections,
+                    start_idx, idx - 1,
+                    y_hat, span_len,
+                    LE[left_section_index + 1],
+                    TE[left_section_index],
+                    TE[left_section_index + 1],
+                    billowing_percentage)
             end
         end
     end
@@ -1279,52 +1257,112 @@ end
 
 
 """
-    catenary_parameter(percentage, d)
+    billowing_arc_length(sections, start_si, end_si, y_hat,
+                         span_len, le_ref, te_left, te_right,
+                         angle_max)
 
-Compute the catenary parameter `a` such that a catenary spanning distance `d`
-has an arc length that is `percentage`% longer than `d`.
+Compute the TE arc length that would result from rotating each
+section's chord around `y_hat` by `angle_max * sin(π t)` radians,
+where `t` is the normalised spanwise position within the rib pair.
 
-Solves `u / sinh(u) = 1 - percentage / 100` where `u = d / (2a)`,
-using Newton's method.
+Non-allocating: uses only stack-allocated MVec3 arithmetic.
+"""
+function billowing_arc_length(
+    sections, start_si, end_si,
+    y_hat, span_len, le_ref,
+    te_left, te_right, angle_max
+)
+    prev_te = MVec3(te_left)
+    arc = 0.0
+    for si in start_si:end_si
+        sec = sections[si]
+        t = dot(sec.LE_point - le_ref, y_hat) / span_len
+        θ = -angle_max * sin(π * t)
+        chord_vec = sec.TE_point - sec.LE_point
+        ct = cos(θ); st = sin(θ)
+        d_y = dot(chord_vec, y_hat)
+        rotated = ct * chord_vec +
+            st * cross(y_hat, chord_vec) +
+            (1 - ct) * d_y * y_hat
+        current_te = sec.LE_point + rotated
+        arc += norm(current_te - prev_te)
+        prev_te = current_te
+    end
+    arc += norm(te_right - prev_te)
+    return arc
+end
 
-# Arguments
-- `percentage::Real`: How much shorter the chord is relative to the arc
-    (0–100). 0 means no sag (flat), larger values mean deeper catenary.
-- `d::Real`: Span distance between the two endpoints.
+"""
+    apply_billowing_to_pair!(sections, start_si, end_si, y_hat,
+                             span_len, le_ref, te_left, te_right,
+                             percentage)
 
-# Returns
-- `Float64`: The catenary parameter `a`.
+Apply billowing to refined sections between two ribs by rotating
+each section's chord around `y_hat`.  Uses Newton iteration to
+find the rotation amplitude that matches the target TE arc-length
+`percentage`, then applies the rotations in-place.
+
+The angle profile is `angle_max * sin(π t)` (zero at ribs,
+maximum at centre).  All angles are in radians.
+
+Non-allocating: modifies `sections[start_si:end_si].TE_point`
+in-place using only stack-allocated MVec3 arithmetic.
 """
-function catenary_parameter(percentage::Real, d::Real)
-    percentage == 0 && return Inf
-    0 < percentage || throw(ArgumentError(
-        "percentage must be ≥ 0, got $percentage"))
-    percentage < 100 || throw(ArgumentError(
-        "percentage must be < 100, got $percentage"))
-    target = 1 - percentage / 100  # u / sinh(u) = target
-    # Initial guess from Taylor: u/sinh(u) ≈ 1 - u²/6
-    u = sqrt(6 * (1 - target))
-    for _ in 1:100
-        f  = u / sinh(u) - target
-        # d/du [u/sinh(u)] = (sinh(u) - u*cosh(u)) / sinh(u)^2
-        sh = sinh(u)
-        df = (sh - u * cosh(u)) / sh^2
-        δ  = f / df
-        u -= δ
-        abs(δ) < 1e-12 && break
+function apply_billowing_to_pair!(
+    sections, start_si, end_si,
+    y_hat, span_len, le_ref,
+    te_left, te_right, percentage
+)
+    percentage <= 0 && return nothing
+    straight = norm(te_right - te_left)
+    straight < 1e-12 && return nothing
+    target_arc = straight / (1 - percentage / 100)
+
+    # Newton iteration on angle_max (rad)
+    angle_max = sqrt(4 * percentage / (100 - percentage))
+    for _ in 1:50
+        arc = billowing_arc_length(
+            sections, start_si, end_si,
+            y_hat, span_len, le_ref,
+            te_left, te_right, angle_max)
+        f = arc - target_arc
+        abs(f) < 1e-12 * target_arc && break
+
+        δ = max(1e-8, abs(angle_max) * 1e-8)
+        arc_p = billowing_arc_length(
+            sections, start_si, end_si,
+            y_hat, span_len, le_ref,
+            te_left, te_right, angle_max + δ)
+        df = (arc_p - arc) / δ
+        abs(df) < 1e-30 && break
+        angle_max -= f / df
     end
-    return d / (2u)
+
+    # Apply converged rotations in-place
+    for si in start_si:end_si
+        sec = sections[si]
+        t = dot(sec.LE_point - le_ref, y_hat) / span_len
+        θ = -angle_max * sin(π * t)
+        chord_vec = sec.TE_point - sec.LE_point
+        ct = cos(θ); st = sin(θ)
+        d_y = dot(chord_vec, y_hat)
+        rotated = ct * chord_vec +
+            st * cross(y_hat, chord_vec) +
+            (1 - ct) * d_y * y_hat
+        sec.TE_point .= sec.LE_point + rotated
+    end
+    return nothing
 end
 
 """
     refine_mesh_with_billowing!(wing; reuse_aero_data)
 
-Refine wing mesh using SPLIT_PROVIDED spacing with catenary TE billowing.
+Refine wing mesh using SPLIT_PROVIDED spacing with TE billowing.
 
-Between each pair of unrefined (rib) sections, the trailing edge follows a
-catenary curve that sags perpendicular to the chord and span (simulating fabric
-billowing between ribs on a ram-air kite). The leading edge stays linearly
-interpolated.
+Between each pair of unrefined (rib) sections, chord vectors are rotated
+around the leading edge to simulate fabric billowing.  The rotation
+amplitude is found iteratively so that the TE arc length matches the
+wing's `billowing_percentage`.
 
 Delegates to [`refine_mesh_by_splitting_provided_sections!`](@ref).
 """

test/Project.toml

@@ -2,10 +2,10 @@
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
@@ -14,16 +14,17 @@ Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+VortexStepMethod = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 
 [compat]
 Aqua = "0.8"
 BenchmarkTools = "1"
 CSV = "0.10"
-DataFrames = "1.7"
-Documenter = "1.8"
 CairoMakie = "0"
 ControlPlots = "0.2.5"
+DataFrames = "1.7"
+Documenter = "1.8"
 Interpolations = "0.15, 0.16"
 LinearAlgebra = "1"
 Logging = "1"

test/wing_geometry/test_billowing.jl

@@ -1,8 +1,7 @@
 using Test
 using LinearAlgebra
 using VortexStepMethod
-using VortexStepMethod: Wing, Section, add_section!, refine!,
-    catenary_parameter
+using VortexStepMethod: Wing, Section, add_section!, refine!
 
 """
     build_flat_wing(n_panels, n_ribs; billowing_percentage=0.0)
@@ -75,7 +74,7 @@ function te_arc_length_between_ribs(wing, rib_left, rib_right)
     return (; arc, straight)
 end
 
-@testset "Billowing (catenary)" begin
+@testset "Billowing" begin
     n_ribs = 5   # 4 rib pairs
     n_panels = 4 * 8  # 8 panels per rib pair
 
@@ -110,6 +109,40 @@ end
         end
     end
 
+    @testset "Billowing sags outward" begin
+        wing = build_flat_wing(n_panels, n_ribs;
+            billowing_percentage=10.0)
+        # For this flat test wing the negative rotation around
+        # y_hat=[0,-1,0] produces -z sag in global coords.
+        # On a real kite this maps to +z in the local airfoil frame.
+        for pair in 1:(n_ribs - 1)
+            le_left = wing.unrefined_sections[pair].LE_point
+            le_right = wing.unrefined_sections[pair + 1].LE_point
+            mid_y = (le_left[2] + le_right[2]) / 2
+            best = nothing
+            best_dist = Inf
+            for sec in wing.refined_sections
+                dist = abs(sec.LE_point[2] - mid_y)
+                if dist < best_dist
+                    best_dist = dist
+                    best = sec
+                end
+            end
+            @test best.TE_point[3] < 0
+        end
+    end
+
+    @testset "Chord length preserved ($pct%)" for pct in
+            [5.0, 10.0, 15.0]
+        wing = build_flat_wing(n_panels, n_ribs;
+            billowing_percentage=pct)
+        for sec in wing.refined_sections
+            chord = norm(sec.TE_point - sec.LE_point)
+            # All ribs have chord=1.0, so interpolated chord=1.0
+            @test isapprox(chord, 1.0; atol=1e-10)
+        end
+    end
+
     @testset "TE arc length matches percentage ($pct%)" for pct in
             [1.0, 5.0, 10.0, 15.0]
         wing = build_flat_wing(n_panels, n_ribs;
@@ -121,28 +154,35 @@ end
             diff = measured - pct
             @info "Pair $pair: measured=$(round(measured; digits=3))%" *
                   " target=$(pct)% diff=$(round(diff; digits=3))%"
-            @test isapprox(measured, pct; rtol=0.05)
+            @test isapprox(measured, pct; rtol=0.01)
         end
     end
 
-    @testset "catenary_parameter edge cases" begin
-        # Zero percentage returns Inf (flat)
-        @test catenary_parameter(0.0, 1.0) == Inf
-        # Invalid inputs throw
-        @test_throws ArgumentError catenary_parameter(-1.0, 1.0)
-        @test_throws ArgumentError catenary_parameter(100.0, 1.0)
-    end
+    @testset "apply_billowing_to_pair! edge cases" begin
+        using VortexStepMethod: apply_billowing_to_pair!, MVec3
+        wing = build_flat_wing(n_panels, n_ribs;
+            billowing_percentage=0.0)
+        y_hat = MVec3(0.0, -1.0, 0.0)
+        le_ref = MVec3(0.0, 2.5, 0.0)
+        te_l = MVec3(-1.0, 0.0, 0.0)
+        te_r = MVec3(-1.0, 2.5, 0.0)
+        orig = [MVec3(s.TE_point)
+                for s in wing.refined_sections]
+
+        # Zero percentage is a no-op
+        apply_billowing_to_pair!(
+            wing.refined_sections, 2, 8,
+            y_hat, 2.5, le_ref, te_l, te_r, 0.0)
+        for (i, sec) in enumerate(wing.refined_sections)
+            @test sec.TE_point == orig[i]
+        end
 
-    @testset "catenary_parameter arc length identity" begin
-        # Verify that catenary_parameter produces the correct
-        # arc length for a given percentage and distance.
-        for pct in [1.0, 5.0, 10.0, 20.0, 40.0]
-            d = 3.0
-            a = catenary_parameter(pct, d)
-            u = d / (2a)
-            arc = 2a * sinh(u)
-            expected_pct = (arc - d) / arc * 100
-            @test isapprox(expected_pct, pct; atol=1e-8)
+        # Coincident ribs (straight ≈ 0) is a no-op
+        apply_billowing_to_pair!(
+            wing.refined_sections, 2, 8,
+            y_hat, 2.5, le_ref, te_l, te_l, 10.0)
+        for (i, sec) in enumerate(wing.refined_sections)
+            @test sec.TE_point == orig[i]
         end
     end
 end