Plot polars (#115)

* add xfoil and timers as normal dependency

* add plotting function for polar data

* fix aqua test

* test bin plot

* test without showing

* improve gitignore

* ignore right files

* add latex kwarg

Author: 1-Bart-1

.gitignore

@@ -2,7 +2,8 @@ Manifest.toml
 .vscode/settings.json
 venv
 results/TUDELFT_V3_LEI_KITE/polars/$C_L$ vs $C_D$.pdf
-*.bin
+data/*.bin
 results/TUDELFT_V3_LEI_KITE/polars/tutorial_testing_stall_model_n_panels_54_distribution_split_provided.pdf
 docs/build/
 results/TUDELFT_V3_LEI_KITE/polars/tutorial_testing_stall_model_n_panels_54_distribution_SPLIT_PROVIDED.pdf
+!test/data/*.bin

Project.toml

@@ -19,6 +19,8 @@ Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Timers = "21f18d07-b854-4dab-86f0-c15a3821819a"
+Xfoil = "19641d66-a62d-11e8-2441-8f57a969a9c4"
 
 [weakdeps]
 ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c"
@@ -29,9 +31,9 @@ VortexStepMethodControlPlotsExt = "ControlPlots"
 [compat]
 Aqua = "0.8"
 BenchmarkTools = "1"
+CSV = "0.10"
 Colors = "0.13"
 ControlPlots = "0.2.5"
-CSV = "0.10"
 DataFrames = "1.7"
 DefaultApplication = "1"
 DelimitedFiles = "1"
@@ -49,8 +51,8 @@ Serialization = "1"
 SharedArrays = "1"
 StaticArrays = "1"
 Statistics = "1"
-Timers = "0.1"
 Test = "1"
+Timers = "0.1"
 Xfoil = "0.5"
 julia = "1.10, 1.11"
 
@@ -63,8 +65,6 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Timers = "21f18d07-b854-4dab-86f0-c15a3821819a"
-Xfoil = "19641d66-a62d-11e8-2441-8f57a969a9c4"
 
 [targets]
-test = ["Test", "DataFrames", "CSV", "Distributed", "Documenter", "Xfoil", "Timers", "BenchmarkTools", "ControlPlots", "Aqua"]
+test = ["Test", "DataFrames", "CSV", "Distributed", "Documenter", "BenchmarkTools", "ControlPlots", "Aqua"]

examples/ram_air_kite.jl

@@ -1,14 +1,6 @@
 using ControlPlots
 using VortexStepMethod
 using LinearAlgebra
-using Pkg
-
-if !("CSV" ∈ keys(Pkg.project().dependencies))
-    using TestEnv
-    TestEnv.activate()
-end
-using CSV
-using DataFrames
 
 PLOT = true
 USE_TEX = false
@@ -44,6 +36,9 @@ vel_app = [
 ] * v_a
 set_va!(body_aero, vel_app)
 
+# Plotting polar data
+PLOT && plot_polar_data(body_aero)
+
 # Plotting geometry
 PLOT && plot_geometry(
     body_aero,

ext/VortexStepMethodControlPlotsExt.jl

@@ -3,7 +3,7 @@ using ControlPlots, LaTeXStrings, VortexStepMethod, LinearAlgebra, Statistics, D
 import ControlPlots: plt
 import VortexStepMethod: calculate_filaments_for_plotting
 
-export plot_wing, plot_circulation_distribution, plot_geometry, plot_distribution, plot_polars, save_plot, show_plot
+export plot_wing, plot_circulation_distribution, plot_geometry, plot_distribution, plot_polars, save_plot, show_plot, plot_polar_data
 
 include("../src/plotting.jl")
 

scripts/polars.jl

@@ -1,7 +1,6 @@
 using Distributed, Timers, Serialization, SharedArrays, StaticArrays
 using Interpolations
 using Xfoil
-using ControlPlots
 using Logging
 
 const SPEED_OF_SOUND = 343 # [m/s] at 20 °C, see: https://en.wikipedia.org/wiki/Speed_of_sound
@@ -133,42 +132,6 @@ try
         return lower_trailing_edge, upper_trailing_edge
     end
 
-    function replace_nan!(matrix)
-        rows, cols = size(matrix)
-        distance = max(rows, cols)
-        for i in 1:rows
-            for j in 1:cols
-                if isnan(matrix[i, j])
-                    neighbors = []
-                    for d in 1:distance
-                        found = false
-                        if i-d >= 1 && !isnan(matrix[i-d, j]);
-                            push!(neighbors, matrix[i-d, j])
-                            found = true
-                        end
-                        if i+d <= rows && !isnan(matrix[i+d, j])
-                            push!(neighbors, matrix[i+d, j])
-                            found = true
-                        end
-                        if j-d >= 1 && !isnan(matrix[i, j-d])
-                            push!(neighbors, matrix[i, j-d])
-                            found = true
-                        end
-                        if j+d <= cols && !isnan(matrix[i, j+d])
-                            push!(neighbors, matrix[i, j+d])
-                            found = true
-                        end
-                        if found; break; end
-                    end
-                    if !isempty(neighbors)
-                        matrix[i, j] = sum(neighbors) / length(neighbors)
-                    end
-                end
-            end
-        end
-        return nothing
-    end
-    
     kite_speed = v_wind
     chord_length = area / width
     local reynolds_number = kite_speed * chord_length / KINEMATIC_VISCOSITY # https://en.wikipedia.org/wiki/Reynolds_number
@@ -208,40 +171,12 @@ try
     cd_matrix = Matrix{Float64}(cd_matrix)
     cm_matrix = Matrix{Float64}(cm_matrix)
 
+    println("Generated lift matrix:")
     display(cl_matrix)
-
-    function plot_values(alphas, d_trailing_edge_angles, matrix, interp, name)
-        fig = plt.figure()
-        ax = fig.add_subplot(projection="3d")
-    
-        X_data = collect(d_trailing_edge_angles) .+ zeros(length(alphas))'
-        Y_data = collect(alphas)' .+ zeros(length(d_trailing_edge_angles))
-    
-        matrix = Matrix{Float64}(matrix)
-        interp_matrix = zeros(size(matrix)...)
-        int_alphas, int_d_trailing_edge_angles = alphas .+ deg2rad(0.5), d_trailing_edge_angles .+ deg2rad(0.5)
-        interp_matrix .= [interp(alpha, d_trailing_edge_angle) for alpha in int_alphas, d_trailing_edge_angle in int_d_trailing_edge_angles]
-        X_int = collect(int_d_trailing_edge_angles) .+ zeros(length(int_alphas))'
-        Y_int = collect(int_alphas)' .+ zeros(length(int_d_trailing_edge_angles))
-    
-        ax.plot_wireframe(X_data, Y_data, matrix, edgecolor="royalblue", lw=0.5, rstride=5, cstride=5, alpha=0.6)
-        ax.plot_wireframe(X_int, Y_int, interp_matrix, edgecolor="orange", lw=0.5, rstride=5, cstride=5, alpha=0.6)
-        plt.xlabel("Alpha")
-        plt.ylabel("Flap angle")
-        plt.zlabel("$name values")
-        plt.title("$name for different d_flap and angle")
-        plt.legend()
-        plt.grid(true)
-        plt.show()
-    end
-    
-    cl_interp = extrapolate(scale(interpolate(cl_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN)
-    cd_interp = extrapolate(scale(interpolate(cd_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN)
-    cm_interp = extrapolate(scale(interpolate(cm_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN)
-    
-    plot_values(alphas, d_trailing_edge_angles, cl_matrix, cl_interp, "Cl")
-    plot_values(alphas, d_trailing_edge_angles, cd_matrix, cd_interp, "Cd")
-    plot_values(alphas, d_trailing_edge_angles, cm_matrix, cm_interp, "Cm")
+    println("Generated drag matrix:")
+    display(cd_matrix)
+    println("Generated moment matrix:")
+    display(cm_matrix)
 
     @info "Relative trailing_edge height: $(upper - lower)"
     @info "Reynolds number for flying speed of $kite_speed is $reynolds_number"

src/VortexStepMethod.jl

@@ -31,7 +31,7 @@ export PanelDistribution, LINEAR, COSINE, COSINE_VAN_GARREL, SPLIT_PROVIDED, UNC
 export InitialGammaDistribution, ELLIPTIC, ZEROS
 export SolverStatus, FEASIBLE, INFEASIBLE, FAILURE
 
-export plot_geometry, plot_distribution, plot_circulation_distribution, plot_geometry, plot_polars, save_plot, show_plot
+export plot_geometry, plot_distribution, plot_circulation_distribution, plot_geometry, plot_polars, save_plot, show_plot, plot_polar_data
 
 # the following functions are defined in ext/VortexStepMethodExt.jl
 function plot_geometry end
@@ -41,6 +41,7 @@ function plot_geometry end
 function plot_polars end
 function save_plot end
 function show_plot end
+function plot_polar_data end
 
 """
    const MVec3    = MVector{3, Float64}

src/plotting.jl

@@ -836,4 +836,75 @@ function VortexStepMethod.plot_polars(
     end
 
     return fig
-end
+end
+
+"""
+    plot_polar_data(body_aero::BodyAerodynamics; alphas=collect(deg2rad.(-5:0.3:25)), beta_tes=collect(deg2rad.(-5:0.3:25)))
+
+Plot polar data (Cl, Cd, Cm) as 3D surfaces against alpha and beta_te angles. Beta_te is the trailing edge deflection angle
+relative to the 2d airfoil or panel chord line.
+
+# Arguments
+- `body_aero`: Wing aerodynamics struct
+
+# Keyword arguments
+- `alphas`: Range of angle of attack values in radians (default: -5° to 25° in 0.3° steps)
+- `beta_tes`: Range of trailing edge angles in radians (default: -5° to 25° in 0.3° steps)
+- `is_show`: Whether to display plots (default: true)
+- `use_tex`: if the external `pdflatex` command shall be used
+"""
+function VortexStepMethod.plot_polar_data(body_aero::BodyAerodynamics; 
+        alphas=collect(deg2rad.(-5:0.3:25)), 
+        beta_tes = collect(deg2rad.(-5:0.3:25)),
+        is_show = true,
+        use_tex = false
+        )
+    if body_aero.panels[1].aero_model === POLAR_MATRICES
+        set_plot_style()
+        
+        # Create figure with subplots
+        fig = plt.figure(figsize=(15, 6))
+        
+        # Get interpolation functions and labels
+        interp_data = [
+            (body_aero.panels[1].cl_interp, L"$C_l$"),
+            (body_aero.panels[1].cd_interp, L"$C_d$"),
+            (body_aero.panels[1].cm_interp, L"$C_m$")
+        ]
+        
+        # Create each subplot
+        for (idx, (interp, label)) in enumerate(interp_data)
+            ax = fig.add_subplot(1, 3, idx, projection="3d")
+            
+            # Create interpolation matrix
+            interp_matrix = zeros(length(alphas), length(beta_tes))
+            interp_matrix .= [interp(alpha, beta_te) for alpha in alphas, beta_te in beta_tes]
+            X = collect(beta_tes) .+ zeros(length(alphas))'
+            Y = collect(alphas)' .+ zeros(length(beta_tes))
+            
+            # Plot surface
+            ax.plot_wireframe(X, Y, interp_matrix, 
+                            edgecolor="blue", 
+                            lw=0.5, 
+                            rstride=5, 
+                            cstride=5, 
+                            alpha=0.6)
+            
+            # Set labels and title
+            ax.set_xlabel(L"$\beta$ [rad]")
+            ax.set_ylabel(L"$\alpha$ [rad]")
+            ax.set_zlabel(label)
+            ax.set_title(label * L" vs $\alpha$ and $\beta$")
+            ax.grid(true)
+        end
+        
+        # Adjust layout and display
+        plt.tight_layout(rect=(0.01, 0.01, 0.99, 0.99))
+        if is_show
+            show_plot(fig)
+        end
+        return fig
+    else
+        throw(ArgumentError("Plotting polar data for $(body_aero.panels[1].aero_model) is not implemented."))
+    end
+end

test/aqua.jl

@@ -1,2 +1,7 @@
 using Aqua
-Aqua.test_all(VortexStepMethod)
+@testset "Aqua.jl" begin
+    Aqua.test_all(
+      VortexStepMethod;
+      stale_deps=(ignore=[:Xfoil, :Timers],),
+    )
+  end