Exact middle twist

Author: 1-Bart-1

src/wing_geometry.jl

@@ -552,22 +552,45 @@ function _apply_refined_section_thetas!(wing::Wing{P, T}, section_thetas) where
     local_y = zeros(MVector{3, T})
     chord = zeros(MVector{3, T})
     normal = zeros(MVector{3, T})
+    n_sec = wing.n_panels + 1
 
-    for i in 1:(wing.n_panels + 1)
+    for i in 1:n_sec
         theta = section_thetas[i]
         section = wing.non_deformed_sections[i]
-
-        if i < wing.n_panels + 1
-            section2 = wing.non_deformed_sections[i + 1]
-            local_y .= normalize(section.LE_point - section2.LE_point)
-        else
-            section_prev = wing.non_deformed_sections[i - 1]
-            local_y .= normalize(section_prev.LE_point - section.LE_point)
+        le = section.LE_point
+
+        # Spanwise axis: average of unit vectors to the left and right neighbour LE
+        # points so curvature does not bias the twist axis. Boundaries fall back to
+        # whichever side exists.
+        local_y .= zero(T)
+        if i > 1
+            le_prev = wing.non_deformed_sections[i - 1].LE_point
+            dx = le_prev[1] - le[1]
+            dy = le_prev[2] - le[2]
+            dz = le_prev[3] - le[3]
+            inv_n = 1 / sqrt(dx * dx + dy * dy + dz * dz)
+            local_y[1] += dx * inv_n
+            local_y[2] += dy * inv_n
+            local_y[3] += dz * inv_n
+        end
+        if i < n_sec
+            le_next = wing.non_deformed_sections[i + 1].LE_point
+            dx = le[1] - le_next[1]
+            dy = le[2] - le_next[2]
+            dz = le[3] - le_next[3]
+            inv_n = 1 / sqrt(dx * dx + dy * dy + dz * dz)
+            local_y[1] += dx * inv_n
+            local_y[2] += dy * inv_n
+            local_y[3] += dz * inv_n
         end
+        inv_n = 1 / sqrt(local_y[1]^2 + local_y[2]^2 + local_y[3]^2)
+        local_y[1] *= inv_n
+        local_y[2] *= inv_n
+        local_y[3] *= inv_n
 
-        chord .= section.TE_point .- section.LE_point
+        chord .= section.TE_point .- le
         normal .= chord × local_y
-        @. wing.refined_sections[i].TE_point = section.LE_point +
+        @. wing.refined_sections[i].TE_point = le +
             cos(theta) * chord - sin(theta) * normal
     end
     return nothing

test/ram_geometry/test_kite_geometry.jl

@@ -190,22 +190,23 @@ using Serialization
         wing = ObjWing(test_obj_path, test_dat_path; remove_nan=true)
         body_aero = BodyAerodynamics([wing])
 
-        # Store original TE point for comparison
-        i = length(body_aero.panels) ÷ 2
+        # Sample the apex panel: with even n_panels, refined section (n_panels/2 + 1)
+        # sits exactly at γ=0, so panels[n_panels/2 + 1].TE_point_1 is the apex TE
+        # where the LE tangent is purely along y and twist preserves y exactly.
+        @test iseven(wing.n_panels)
+        i = wing.n_panels ÷ 2 + 1
         original_te_point = copy(body_aero.panels[i].TE_point_1)
 
-        # Apply deformation with non-zero angles
-        theta_dist = fill(deg2rad(30.0), wing.n_panels)  # 30 degrees twist for all panels
-        delta_dist = fill(deg2rad(5.0), wing.n_panels)   # 5 degrees TE deflection for all panels
+        theta_dist = fill(deg2rad(30.0), wing.n_panels)
+        delta_dist = fill(deg2rad(5.0), wing.n_panels)
 
         VortexStepMethod.deform!(wing, theta_dist, delta_dist)
         VortexStepMethod.reinit!(body_aero)
 
-        # Check if TE point changed after deformation
         deformed_te_point = copy(body_aero.panels[i].TE_point_1)
         @test !isapprox(original_te_point, deformed_te_point, atol=1e-2)
         @test deformed_te_point[3] < original_te_point[3] # right hand rule
-        @test deformed_te_point[2] ≈ original_te_point[2] atol=1e-2 # right hand rule
+        @test deformed_te_point[2] ≈ original_te_point[2] atol=1e-5 # right hand rule
         @test deformed_te_point[1] < original_te_point[1] # right hand rule
         @test body_aero.panels[i].delta ≈ deg2rad(5.0)