diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 8b9ba632..cc820141 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,6 +34,9 @@ jobs: arch: x64 version: 1 steps: + - name: Install matplotlib + run: if [ "$RUNNER_OS" = "Linux" ]; then sudo apt-get install -y python3-matplotlib; fi + shell: bash - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 with: diff --git a/CITATION.cff b/CITATION.cff index 3e23fd7e..828e5bf2 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -1,21 +1,20 @@ -cff-version: 0.1.0 +cff-version: 1.2.0 message: "If you are using this software, please cite it as shown below." authors: - - - family-names: "Cayon" + - family-names: "Cayon" given-names: "Oriol" orcid: "https://orcid.org/0000-0002-2065-8673" - family-names: "Poland" + - family-names: "Poland" given-names: "Jelle Agatho Wilhelm" orcid: "https://orcid.org/0000-0003-3164-5648" -title: "VortexStepMethod" +title: "VortexStepMethod.jl" keywords: - - Airborne Wind Energy -version: 0.1.0 + - Airborne Wind Energy, AWE, AWES, Vortex step methods +version: 1.0.0 # doi: "11.1111/11111" # date-released: YYYY-MM-DD license: MIT -url: "https://github.com/ocayon/Vortex-Step-Method" +url: "https://github.com/Albatross-Kite-Transport/VortexStepMethod.jl" preferred-citation: type: article authors: @@ -31,4 +30,3 @@ preferred-citation: issue: 7 volume: 16 year: 2023 - article-number: 3061 diff --git a/Project.toml b/Project.toml index 37774598..46007f39 100644 --- a/Project.toml +++ b/Project.toml @@ -6,20 +6,20 @@ version = "0.1.0" [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] +ControlPlots = "0.2.5" Measures = "0.3" StaticArrays = "1" -Plots = "1" Statistics = "1" DelimitedFiles = "1" BenchmarkTools = "1" diff --git a/README.md b/README.md index 6001ed34..9a4c2810 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ [![Build Status](https://github.com/Albatross-Kite-Transport/VortexStepMethod.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/Albatross-Kite-Transport/VortexStepMethod.jl/actions/workflows/CI.yml?query=branch%3Amain) [![Coverage](https://codecov.io/gh/Albatross-Kite-Transport/VortexStepMethod.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/Albatross-Kite-Transport/VortexStepMethod.jl) +======= + # Simulation of a 3D airfoil using the Vortex Step Method The Vortex Step Method (VSM) is an enhanced lifting line method that improves upon the classic approach by solving the @@ -116,7 +118,7 @@ Surfplan files can be converted to an input for `VortexStepMethod.jl` using the ## Output - CL, CD, CS (side force coefficient) - the spanwise distribution of forces - --> moment coefficients (not yet implemented) + --> moment coefficients (will be implemented in release 1.1) ## Citation If you use this project in your research, please consider citing it. diff --git a/Spanwise Distributions.pdf b/Spanwise Distributions.pdf deleted file mode 100644 index 79fea672..00000000 Binary files a/Spanwise Distributions.pdf and /dev/null differ diff --git a/docs/geometry.png b/docs/geometry.png new file mode 100644 index 00000000..4510f965 Binary files /dev/null and b/docs/geometry.png differ diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl index 2cf567d5..25dada2b 100644 --- a/examples/rectangular_wing.jl +++ b/examples/rectangular_wing.jl @@ -1,5 +1,5 @@ using LinearAlgebra -using Plots +using ControlPlots using VortexStepMethod # Step 1: Define wing parameters @@ -55,36 +55,33 @@ println("Projected area = $(round(results_vsm["projected_area"], digits=4)) m²" # Step 6: Plot geometry plot_geometry( wa, - "rectangular_wing_geometry"; + "Rectangular_wing_geometry"; data_type=".pdf", save_path=".", is_save=false, is_show=true, ) +nothing -# # Step 7: Plot spanwise distributions -# y_coordinates = [panel.aerodynamic_center[2] for panel in wa.panels] +# Step 7: Plot spanwise distributions +y_coordinates = [panel.aerodynamic_center[2] for panel in wa.panels] -# plot_distribution( -# [y_coordinates, y_coordinates], -# [results_vsm, results_llt], -# ["VSM", "LLT"], -# title="Spanwise Distributions" -# ) - -# # Step 8: Plot polar curves -# angle_range = range(0, 20, 20) -# plot_polars( -# [llt_solver, vsm_solver], -# [wa, wa], -# ["LLT", "VSM"], -# angle_range=angle_range, -# angle_type="angle_of_attack", -# Umag=Umag, -# title="Rectangular Wing Polars" -# ) +plot_distribution( + [y_coordinates, y_coordinates], + [results_vsm, results_llt], + ["VSM", "LLT"], + title="Spanwise Distributions" +) -# Save plots if needed -# savefig("geometry.pdf") -# savefig("distributions.pdf") -# savefig("polars.pdf") \ No newline at end of file +# Step 8: Plot polar curves +angle_range = range(0, 20, 20) +plot_polars( + [llt_solver, vsm_solver], + [wa, wa], + ["LLT", "VSM"], + angle_range=angle_range, + angle_type="angle_of_attack", + Umag=Umag, + title="Rectangular Wing Polars" +) +nothing diff --git a/examples/testing_stall_model.jl b/examples/testing_stall_model.jl index 25ab07ec..c11473e8 100644 --- a/examples/testing_stall_model.jl +++ b/examples/testing_stall_model.jl @@ -3,7 +3,6 @@ using VortexStepMethod using CSV using DataFrames using LinearAlgebra -using Plots # Find root directory root_dir = dirname(@__DIR__) diff --git a/results/TUDELFT_V3_LEI_KITE/polars/$C_L$ vs $C_D$.pdf b/results/TUDELFT_V3_LEI_KITE/polars/$C_L$ vs $C_D$.pdf new file mode 100644 index 00000000..a8b97f81 Binary files /dev/null and b/results/TUDELFT_V3_LEI_KITE/polars/$C_L$ vs $C_D$.pdf differ diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 42b4b51c..fdecf0eb 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -6,7 +6,7 @@ using Logging using Statistics using Colors using DelimitedFiles -using Plots +using ControlPlots using Measures using LaTeXStrings @@ -18,6 +18,7 @@ export calculate_results, solve_circulation_distribution export add_section!, set_va! export calculate_span, calculate_projected_area export plot_wing, plot_circulation_distribution, plot_geometry, plot_distribution, plot_polars +export show_plot """ const MVec3 = MVector{3, Float64} @@ -31,8 +32,8 @@ const MVec3 = MVector{3, Float64} Position vector, either a `MVec3` or a `Vector` for use in function signatures. """ -const PosVector=Union{MVec3, Vector} -const VelVector=Union{MVec3, Vector} +const PosVector=Union{MVec3, Vector, SizedVector{3, Float64, Vector{Float64}}} +const VelVector=Union{MVec3, Vector, SizedVector{3, Float64, Vector{Float64}}} # Include core functionality include("wing_geometry.jl") diff --git a/src/color_palette.jl b/src/color_palette.jl index 615b7b0a..5722f42f 100644 --- a/src/color_palette.jl +++ b/src/color_palette.jl @@ -69,41 +69,37 @@ end Set the default style for plots using LaTeX. """ -function set_plot_style() - default( - fontfamily="Computer Modern", - titlefontsize=14, - guidefontsize=12, - tickfontsize=10, - legendfontsize=10, - linewidth=1, - markersize=6, - size=(600, 400), - grid=true, - gridstyle=:dash, - gridcolor=:gray, - gridalpha=0.3 - ) - - # LaTeX settings if using PGFPlotsX backend - if backend() isa PGFPlotsX - pgfplotsx( - tex_output_standalone=true, - pdf_engine="pdflatex", - latex_engine="pdflatex" - ) - end +function set_plot_style(titel_size=16) + # plt.style.use('seaborn-whitegrid') + # plt.style.use("seaborn-v0_8-whitegrid") + rcParams = plt.PyDict(plt.matplotlib."rcParams") + rcParams["text.usetex"] = true + rcParams["font.family"] = "serif" + rcParams["font.serif"] = ["Computer Modern Roman"] + rcParams["axes.titlesize"] = titel_size + # rcParams["axes.ymargin"] = 0.1 + rcParams["axes.labelsize"] = 12 + rcParams["axes.linewidth"] = 1 + rcParams["lines.linewidth"] = 1 + rcParams["lines.markersize"] = 6 + rcParams["xtick.labelsize"] = 10 + rcParams["ytick.labelsize"] = 10 + rcParams["legend.fontsize"] = 10 + rcParams["figure.titlesize"] = 16 + rcParams["pgf.texsystem"] = "pdflatex" # Use pdflatex + rcParams["pgf.rcfonts"] = false + rcParams["figure.figsize"] = (10, 6) # Default figure size end -""" - apply_palette!(p::Plots.Plot, colors::Vector{String}) +# """ +# apply_palette!(p::Plots.Plot, colors::Vector{String}) -Apply the color palette to a plot. -""" -function apply_palette!(p::Plots.Plot, colors::Vector{String}) - for (i, series) in enumerate(p.series_list) - color_name = colors[mod1(i, length(colors))] - series.plotattributes[:linecolor] = get_color(color_name) - end - return p -end \ No newline at end of file +# Apply the color palette to a plot. +# """ +# function apply_palette!(p::Plots.Plot, colors::Vector{String}) +# for (i, series) in enumerate(p.series_list) +# color_name = colors[mod1(i, length(colors))] +# series.plotattributes[:linecolor] = get_color(color_name) +# end +# return p +# end \ No newline at end of file diff --git a/src/panel.jl b/src/panel.jl index f572cac0..49d803e5 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -25,10 +25,10 @@ Represents a panel in a vortex step method simulation. - `filaments::Vector{BoundFilament}`: Panel filaments """ mutable struct Panel - TE_point_1::Vector{Float64} - LE_point_1::Vector{Float64} - TE_point_2::Vector{Float64} - LE_point_2::Vector{Float64} + TE_point_1::MVec3 + LE_point_1::MVec3 + TE_point_2::MVec3 + LE_point_2::MVec3 chord::Float64 va::Union{Nothing,Vector{Float64}} corner_points::Matrix{Float64} @@ -37,13 +37,13 @@ mutable struct Panel cd_coefficients::Union{Nothing,Vector{Float64}} cm_coefficients::Union{Nothing,Vector{Float64}} panel_polar_data::Union{Nothing,Matrix{Float64}} - aerodynamic_center::Vector{Float64} - control_point::Vector{Float64} - bound_point_1::Vector{Float64} - bound_point_2::Vector{Float64} - x_airf::Vector{Float64} - y_airf::Vector{Float64} - z_airf::Vector{Float64} + aerodynamic_center::MVec3 + control_point::MVec3 + bound_point_1::MVec3 + bound_point_2::MVec3 + x_airf::MVec3 + y_airf::MVec3 + z_airf::MVec3 width::Float64 filaments::Vector{Union{BoundFilament, SemiInfiniteFilament}} @@ -54,9 +54,9 @@ mutable struct Panel control_point::PosVector, bound_point_1::PosVector, bound_point_2::PosVector, - x_airf::Vector{Float64}, - y_airf::Vector{Float64}, - z_airf::Vector{Float64} + x_airf::PosVector, + y_airf::PosVector, + z_airf::PosVector ) # Initialize basic geometry TE_point_1 = copy(section_1.TE_point) diff --git a/src/plotting.jl b/src/plotting.jl index 3f30f5c2..a00acedf 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -19,9 +19,9 @@ function save_plot(fig, save_path, title; data_type=".pdf") @debug "Current working directory: $(pwd())" try - savefig(fig, full_path) - @debug "Figure saved successfully" - + fig.savefig(full_path) + @debug "Figure saved as $data_type" + if isfile(full_path) @debug "File successfully saved to $full_path" @debug "File size: $(filesize(full_path)) bytes" @@ -33,6 +33,13 @@ function save_plot(fig, save_path, title; data_type=".pdf") @error "Error type: $(typeof(e))" rethrow(e) end + +# if os.path.exists(full_path): +# logging.debug(f"File successfully saved to {full_path}") +# logging.debug(f"File size: {os.path.getsize(full_path)} bytes") +# else: +# logging.info(f"File does not exist after save attempt: {full_path}") +# end end """ @@ -45,7 +52,7 @@ Display a plot at specified DPI. - `dpi`: Dots per inch for the figure (default: 130) """ function show_plot(fig; dpi=130) - display(fig) + plt.display(fig) end """ @@ -61,7 +68,7 @@ Plot a line segment in 3D with arrow. - `width`: Line width (default: 3) """ function plot_line_segment!(ax, segment, color, label; width=3) - plot!(ax, + ax.plot( [segment[1][1], segment[2][1]], [segment[1][2], segment[2][2]], [segment[1][3], segment[2][3]], @@ -69,9 +76,9 @@ function plot_line_segment!(ax, segment, color, label; width=3) ) dir = segment[2] - segment[1] - quiver!(ax, + ax.quiver( [segment[1][1]], [segment[1][2]], [segment[1][3]], - quiver=([dir[1]], [dir[2]], [dir[3]]), + [dir[1]], [dir[2]], [dir[3]], color=color ) end @@ -81,10 +88,10 @@ end Set 3D plot axes to equal scale. """ -function set_axes_equal!(ax) - x_lims = xlims(ax) - y_lims = ylims(ax) - z_lims = zlims(ax) +function set_axes_equal!(ax; zoom=1.8) + x_lims = ax.get_xlim3d() ./ zoom + y_lims = ax.get_ylim3d() ./ zoom + z_lims = ax.get_zlim3d() ./ zoom x_range = abs(x_lims[2] - x_lims[1]) y_range = abs(y_lims[2] - y_lims[1]) @@ -96,9 +103,9 @@ function set_axes_equal!(ax) y_mid = mean(y_lims) z_mid = mean(z_lims) - xlims!(ax, (x_mid - max_range/2, x_mid + max_range/2)) - ylims!(ax, (y_mid - max_range/2, y_mid + max_range/2)) - zlims!(ax, (z_mid - max_range/2, z_mid + max_range/2)) + ax.set_xlim3d((x_mid - max_range/2, x_mid + max_range/2)) + ax.set_ylim3d((y_mid - max_range/2, y_mid + max_range/2)) + ax.set_zlim3d((z_mid - max_range/2, z_mid + max_range/2)) end """ @@ -106,7 +113,9 @@ end Create a 3D plot of wing geometry including panels and filaments. """ -function create_geometry_plot(wing_aero, title, view_elevation, view_azimuth) +function create_geometry_plot(wing_aero, title, view_elevation, view_azimuth; zoom=1.8) + set_plot_style(28) + panels = wing_aero.panels va = isa(wing_aero.va, Tuple) ? wing_aero.va[1] : wing_aero.va @@ -116,11 +125,9 @@ function create_geometry_plot(wing_aero, title, view_elevation, view_azimuth) aerodynamic_centers = [panel.aerodynamic_center for panel in panels] # Create plot - plt = plot3d( - title=title, - size=(800, 800), - camera=(view_azimuth, view_elevation) - ) + fig = plt.figure(figsize=(14, 14)) + ax = fig.add_subplot(111, projection="3d") + ax.set_title(title) # Plot panels legend_used = Dict{String,Bool}() @@ -135,14 +142,17 @@ function create_geometry_plot(wing_aero, title, view_elevation, view_azimuth) push!(y_corners, y_corners[1]) push!(z_corners, z_corners[1]) - plot!(plt, x_corners, y_corners, z_corners, - color=:grey, linewidth=1, - label=i == 1 ? "Panel Edges" : "") + ax.plot(x_corners, + y_corners, + z_corners, + color=:grey, + linewidth=1, + label=i == 1 ? "Panel Edges" : "") # Plot control points and aerodynamic centers - scatter!(plt, [control_points[i][1]], [control_points[i][2]], [control_points[i][3]], + ax.scatter([control_points[i][1]], [control_points[i][2]], [control_points[i][3]], color=:green, label=i == 1 ? "Control Points" : "") - scatter!(plt, [aerodynamic_centers[i][1]], [aerodynamic_centers[i][2]], [aerodynamic_centers[i][3]], + ax.scatter([aerodynamic_centers[i][1]], [aerodynamic_centers[i][2]], [aerodynamic_centers[i][3]], color=:blue, label=i == 1 ? "Aerodynamic Centers" : "") # Plot filaments @@ -153,7 +163,7 @@ function create_geometry_plot(wing_aero, title, view_elevation, view_azimuth) x1, x2, color = filament @debug "Legend: $legend" show_legend = !get(legend_used, legend, false) - plot_line_segment!(plt, [x1, x2], color, show_legend ? legend : "") + plot_line_segment!(ax, [x1, x2], color, show_legend ? legend : "") legend_used[legend] = true end end @@ -163,15 +173,27 @@ function create_geometry_plot(wing_aero, title, view_elevation, view_azimuth) va_mag = norm(va) va_vector_begin = -2 * max_chord * va / va_mag va_vector_end = va_vector_begin + 1.5 * va / va_mag - plot_line_segment!(plt, [va_vector_begin, va_vector_end], :lightblue, "va") + plot_line_segment!(ax, [va_vector_begin, va_vector_end], :lightblue, "va") + + # Add legends for the first occurrence of each label + handles, labels = ax.get_legend_handles_labels() + by_label = Dict(zip(labels, handles)) + ax.legend(values(by_label), keys(by_label), bbox_to_anchor = (0,0,1.1,1)) # Set labels and make axes equal - xlabel!(plt, "x") - ylabel!(plt, "y") - zlabel!(plt, "z") - set_axes_equal!(plt) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + set_axes_equal!(ax; zoom) + + # Set the initial view + ax.view_init(elev=view_elevation, azim=view_azimuth) + + # Ensure the figure is fully rendered + # fig.canvas.draw() + plt.tight_layout(rect=(0,0,1,0.97)) - return plt + return fig end """ @@ -188,6 +210,7 @@ function plot_geometry(wing_aero, title; view_azimuth=-120) if is_save + plt.ioff() # Angled view fig = create_geometry_plot(wing_aero, "$(title)_angled_view", 15, -120) save_plot(fig, save_path, "$(title)_angled_view", data_type=data_type) @@ -206,8 +229,9 @@ function plot_geometry(wing_aero, title; end if is_show + plt.ion() fig = create_geometry_plot(wing_aero, title, view_elevation, view_azimuth) - display(fig) + plt.display(fig) end end @@ -237,97 +261,142 @@ function plot_distribution(y_coordinates_list, results_list, label_list; "Number of results ($(length(results_list))) must match number of labels ($(length(label_list)))" )) - # Create plot with layout - plt = plot( - layout=(3,3), - size=(1200, 800), - plot_title=title - ) + # Set the plot style + set_plot_style() - # CL Distribution - plot!(plt[1], - title=L"C_L Distribution", - xlabel="Spanwise Position y/b", - ylabel=L"Lift Coefficient C_L" - ) - for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - plot!(plt[1], y_coords, results["cl_distribution"], - label="$label CL: $(round(results["cl"], digits=2))") - end + # Initializing plot + fig, axs = plt.subplots(3, 3, figsize=(16, 10)) + fig.suptitle(title, fontsize=16) - # CD Distribution - plot!(plt[2], - title=L"C_D Distribution", - xlabel="Spanwise Position y/b", - ylabel=L"Drag Coefficient C_D" - ) - for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - plot!(plt[2], y_coords, results["cd_distribution"], - label="$label CD: $(round(results["cd"], digits=2))") + # CL plot + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + value = "$(round(result_i["cl"], digits=2))" + if label_i == "LLT" + label = label_i * L" $~C_L$: " * value + else + label = label_i * L" $C_L$: " * value + end + axs[1, 1].plot( + y_coordinates_i, + result_i["cl_distribution"], + label=label + ) end + axs[1, 1].set_title(L"$C_L$ Distribution", size=16) + axs[1, 1].set_xlabel(L"Spanwise Position $y/b$") + axs[1, 1].set_ylabel(L"Lift Coefficient $C_L$") + axs[1, 1].legend() + + # CD plot + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + value = "$(round(result_i["cl"], digits=2))" + if label_i == "LLT" + label = label_i * L" $~C_D$: " * value + else + label = label_i * L" $C_D$: " * value + end + axs[1, 2].plot( + y_coordinates_i, + result_i["cd_distribution"], + label=label + ) + end + axs[1, 2].set_title(L"$C_D$ Distribution", size=16) + axs[1, 2].set_xlabel(L"Spanwise Position $y/b$") + axs[1, 2].set_ylabel(L"Drag Coefficient $C_D$") + axs[1, 2].legend() # Gamma Distribution - plot!(plt[3], - title=L"\Gamma Distribution", - xlabel="Spanwise Position y/b", - ylabel=L"Circulation \Gamma" - ) - for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - plot!(plt[3], y_coords, results["gamma_distribution"], label=label) + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + axs[1, 3].plot( + y_coordinates_i, + result_i["gamma_distribution"], + label=label_i + ) end + axs[1, 3].set_title(L"\Gamma~Distribution", size=16) + axs[1, 3].set_xlabel(L"Spanwise Position $y/b$") + axs[1, 3].set_ylabel(L"Circulation~\Gamma") + axs[1, 3].legend() # Geometric Alpha - plot!(plt[4], - title=L"\alpha Geometric", - xlabel="Spanwise Position y/b", - ylabel=L"Angle of Attack \alpha (deg)" - ) - for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - plot!(plt[4], y_coords, rad2deg.(results["alpha_geometric"]), label=label) + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + axs[2, 1].plot( + y_coordinates_i, + result_i["alpha_geometric"], + label=label_i + ) end - - # Corrected Alpha - plot!(plt[5], - title=L"\alpha Corrected", - xlabel="Spanwise Position y/b", - ylabel=L"Angle of Attack \alpha (deg)" - ) - for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - plot!(plt[5], y_coords, rad2deg.(results["alpha_at_ac"]), label=label) + axs[2, 1].set_title(L"$\alpha$ Geometric", size=16) + axs[2, 1].set_xlabel(L"Spanwise Position $y/b$") + axs[2, 1].set_ylabel(L"Angle of Attack $\alpha$ (deg)") + axs[2, 1].legend() + + # Calculated/ Corrected Alpha + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + axs[2, 2].plot( + y_coordinates_i, + result_i["alpha_at_ac"], + label=label_i + ) + end + axs[2, 2].set_title(L"$\alpha$ result (corrected to aerodynamic center)", size=16) + axs[2, 2].set_xlabel(L"Spanwise Position $y/b$") + axs[2, 2].set_ylabel(L"Angle of Attack $\alpha$ (deg)") + axs[2, 2].legend() + + # Uncorrected Alpha plot + for (y_coordinates_i, result_i, label_i) in zip(y_coordinates_list, results_list, label_list) + axs[2, 3].plot( + y_coordinates_i, + result_i["alpha_uncorrected"], + label=label_i + ) end + axs[2, 3].set_title(L"$\alpha$ Uncorrected (if VSM, at the control point)", size=16) + axs[2, 3].set_xlabel(L"Spanwise Position $y/b$") + axs[2, 3].set_ylabel(L"Angle of Attack $\alpha$ (deg)") + axs[2, 3].legend() # Force Components for (idx, component) in enumerate(["x", "y", "z"]) - plot!(plt[6+idx], - title="Force in $component direction", - xlabel="Spanwise Position y/b", - ylabel=L"F_%$component" - ) + axs[3, idx].set_title("Force in $component direction", size=16) + axs[3, idx].set_xlabel(L"Spanwise Position $y/b$") + axs[3, idx].set_ylabel(raw"$F_\mathrm" * "{$component}" * raw"$") for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) # Extract force components for the current direction (idx) forces = results["F_distribution"][idx, :] - # Verify dimensions match if length(y_coords) != length(forces) @warn "Dimension mismatch in force plotting" length(y_coords) length(forces) component continue # Skip this component instead of throwing error end - - plot!(plt[6+idx], y_coords, forces, - label="$label ΣF_$component: $(round(results["F$component"], digits=2))N") + space="" + if label == "LLT" + space = "~" + end + axs[3, idx].plot( + y_coords, + forces, + label = "$label" * space * raw"$~\Sigma~F_\mathrm" * "{$component}:~" * + raw"$" * "$(round(results["F$component"], digits=2)) N" + ) + axs[3, idx].legend() end end - # Save and show plot - if is_save - save_plot(plt, save_path, title, data_type=data_type) - end + fig.tight_layout() + + # # Save and show plot + # if is_save + # save_plot(fig, save_path, title, data_type=data_type) + # end if is_show - show_plot(plt) + show_plot(fig) end - return plt + return fig end """ @@ -373,13 +442,15 @@ function generate_polar_data( # Previous gamma for initialization gamma = nothing + + set_plot_style() for (i, angle_i) in enumerate(angle_range) # Set angle based on type if angle_type == "angle_of_attack" α = deg2rad(angle_i) β = side_slip - elseif angle_type == "side_slip" + elseif angle_type == raw"side_slip" α = angle_of_attack β = deg2rad(angle_i) else @@ -486,92 +557,189 @@ function plot_polars( # Read all data first data = readdlm(path, ',') # Skip the header row by taking data from row 2 onwards + @show data[1, :] data = data[2:end, :] - push!(polar_data_list, [data[:,3], data[:,1], data[:,2]]) + push!(polar_data_list, [data[:,3], data[:,1], data[:,2], zeros(length(data[:,1]))]) end end - # Create plots with 2x2 layout - plt = plot( - layout=(2,2), - size=(1000, 1000), - plot_title=title - ) + # Initializing plot + fig, axs = plt.subplots(2, 2, figsize=(14, 14)) # Number of computational results (excluding literature) n_solvers = length(solver_list) - - # Plot CL vs angle - plot!(plt[1]) for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) - style = i ≤ n_solvers ? (:solid, :star, 7) : (:solid, :circle, 5) - plot!(plt[1], polar_data[1], polar_data[2], - label=label, linestyle=style[1], marker=style[2], markersize=style[3]) - + if i < n_solvers + linestyle = "-" + marker = "*" + markersize = 7 + else + linestyle = "-" + marker = "." + markersize = 5 + end + if contains(label, "LLT") + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => raw"~") + label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") + label = raw"$" * label * raw"$" + else + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => "~") + label = replace(label, "VSM" => raw"\mathrm{VSM}") + label = raw"$" * label * raw"$" + end + axs[1, 1].plot( + polar_data[1], + polar_data[2], + label=label, + linestyle=linestyle, + marker=marker, + markersize=markersize, + ) # Limit y-range if CL > 10 if maximum(polar_data[2]) > 10 - ylims!(plt[1], (-0.5, 2.0)) + axs[1, 1].set_ylim([-0.5, 2]) end + title = raw"$C_L" * raw"$" * " vs $angle_type [°]" + axs[1, 1].set_title(title) + axs[1, 1].set_xlabel("$angle_type [°]") + axs[1, 1].set_ylabel(L"$C_L$") + axs[1, 1].legend() end - title!(plt[1], L"C_L \textrm{ vs } %$angle_type") - xlabel!(plt[1], "$angle_type [deg]") - ylabel!(plt[1], L"C_L") - # Plot CD vs angle - plot!(plt[2]) for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) - style = i ≤ n_solvers ? (:solid, :star, 7) : (:solid, :circle, 5) - plot!(plt[2], polar_data[1], polar_data[3], - label=label, linestyle=style[1], marker=style[2], markersize=style[3]) - - # Limit y-range if CD > 10 - if maximum(polar_data[3]) > 10 - ylims!(plt[2], (-0.2, 0.5)) + if i < n_solvers + linestyle = "-" + marker = "*" + markersize = 7 + else + linestyle = "-" + marker = "." + markersize = 5 end + if contains(label, "LLT") + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => raw"~") + label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") + label = raw"$" * label * raw"$" + else + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => "~") + label = replace(label, "VSM" => raw"\mathrm{VSM}") + label = raw"$" * label * raw"$" + end + axs[1, 2].plot( + polar_data[1], + polar_data[3], + label=label, + linestyle=linestyle, + marker=marker, + markersize=markersize, + ) + # Limit y-range if CL > 10 + if maximum(polar_data[2]) > 10 + axs[1, 2].set_ylim([-0.5, 2]) + end + title = raw"$C_D" * raw"$" * " vs $angle_type [°]" + axs[1, 2].set_title(title) + axs[1, 2].set_xlabel("$angle_type [°]") + axs[1, 2].set_ylabel(L"$C_D$") + axs[1, 2].legend() end - title!(plt[2], L"C_D \textrm{ vs } %$angle_type") - xlabel!(plt[2], "$angle_type [deg]") - ylabel!(plt[2], L"C_D") - # Plot CS vs angle (if available) - plot!(plt[3]) + for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) - # Check if CS data is available (length > 3) - if length(polar_data) > 3 - style = i ≤ n_solvers ? (:solid, :star, 7) : (:solid, :circle, 5) - plot!(plt[3], polar_data[1], polar_data[4], - label=label, linestyle=style[1], marker=style[2], markersize=style[3]) + if i < n_solvers + linestyle = "-" + marker = "*" + markersize = 7 + else + linestyle = "-" + marker = "." + markersize = 5 + end + if contains(label, "LLT") + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => raw"~") + label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") + label = raw"$" * label * raw"$" + else + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => "~") + label = replace(label, "VSM" => raw"\mathrm{VSM}") + label = raw"$" * label * raw"$" end + axs[2, 1].plot( + polar_data[1], + polar_data[4], + label=label, + linestyle=linestyle, + marker=marker, + markersize=markersize, + ) + # Limit y-range if CL > 10 + if maximum(polar_data[2]) > 10 + axs[2, 1].set_ylim([-0.5, 2]) + end + title = raw"$C_S" * raw"$" * " vs $angle_type [°]" + axs[2, 1].set_title(title) + axs[2, 1].set_xlabel("$angle_type [°]") + axs[2, 1].set_ylabel(L"$C_S$") + axs[2, 1].legend() end - title!(plt[3], L"C_S \textrm{ vs } %$angle_type") - xlabel!(plt[3], "$angle_type [deg]") - ylabel!(plt[3], L"C_S") - # Plot CL vs CD - plot!(plt[4]) for (i, (polar_data, label)) in enumerate(zip(polar_data_list, label_list)) - style = i ≤ n_solvers ? (:solid, :star, 7) : (:solid, :circle, 5) - plot!(plt[4], polar_data[2], polar_data[3], # Note: CD on x-axis, CL on y-axis - label=label, linestyle=style[1], marker=style[2], markersize=style[3]) - - # Limit ranges if values > 10 + if i < n_solvers + linestyle = "-" + marker = "*" + markersize = 7 + else + linestyle = "-" + marker = "." + markersize = 5 + end + if contains(label, "LLT") + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => raw"~") + label = replace(label, "LLT" => raw"\mathrm{LLT}{~\,}") + label = raw"$" * label * raw"$" + else + label = replace(label, "e5" => raw"\cdot10^5") + label = replace(label, " " => "~") + label = replace(label, "VSM" => raw"\mathrm{VSM}") + label = raw"$" * label * raw"$" + end + axs[2, 2].plot( + polar_data[2], + polar_data[3], + label=label, + linestyle=linestyle, + marker=marker, + markersize=markersize, + ) + # Limit y-range if CL > 10 if maximum(polar_data[2]) > 10 || maximum(polar_data[3]) > 10 - ylims!(plt[4], (-0.5, 2.0)) - xlims!(plt[4], (-0.2, 0.5)) + axs[2, 2].set_ylim([-0.5, 2]) + axs[2, 2].set_xlim([-0.5, 2]) end + title = raw"$C_L" * raw"$" * " vs " * raw"$C_D" * raw"$" + axs[2, 2].set_title(title) + axs[2, 2].set_xlabel(L"$C_D$") + axs[2, 2].set_ylabel(L"$C_L$") + axs[2, 2].legend() end - title!(plt[4], L"C_L \textrm{ vs } C_D \textrm{ (over } %$angle_type \textrm{ range)}") - xlabel!(plt[4], L"C_D") - ylabel!(plt[4], L"C_L") + fig.tight_layout(h_pad=2.5, rect=(0.01,0.01,0.99,0.99)) + # Save and show plot - if is_save - save_plot(plt, save_path, title, data_type=data_type) + if is_save && !isnothing(save_path) + save_plot(fig, save_path, title, data_type=data_type) end - + if is_show - show_plot(plt) + show_plot(fig) end - return plt + return fig end \ No newline at end of file diff --git a/src/wing_aerodynamics.jl b/src/wing_aerodynamics.jl index 21bdec96..c6c74d61 100644 --- a/src/wing_aerodynamics.jl +++ b/src/wing_aerodynamics.jl @@ -502,11 +502,24 @@ function calculate_results(wa::WingAerodynamics, lift_wing_3D_sum += lift_prescribed_va * panel.width drag_wing_3D_sum += drag_prescribed_va * panel.width side_wing_3D_sum += side_prescribed_va * panel.width + + # TODO make this work + # fx_global_3D_sum += fx_global_3D + # fy_global_3D_sum += fy_global_3D + # fz_global_3D_sum += fz_global_3D # Store coefficients push!(cl_prescribed_va, lift_prescribed_va / (q_inf * panel.chord)) push!(cd_prescribed_va, drag_prescribed_va / (q_inf * panel.chord)) push!(cs_prescribed_va, side_prescribed_va / (q_inf * panel.chord)) + + # TODO translate this + # fx_global_3D_list.append(fx_global_3D) + # fy_global_3D_list.append(fy_global_3D) + # fz_global_3D_list.append(fz_global_3D) + # f_global_3D_list.append( + # np.array([fx_global_3D, fy_global_3D, fz_global_3D]) + # ) end if is_only_f_and_gamma_output diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index be662d19..1f305e8c 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -13,8 +13,8 @@ Represents a wing section with leading edge, trailing edge, and aerodynamic prop - `["lei_airfoil_breukels", [d_tube,camber]]`: LEI airfoil with Breukels parameters """ struct Section - LE_point::Vector{Float64} - TE_point::Vector{Float64} + LE_point::MVec3 + TE_point::MVec3 aero_input::Union{String, Tuple{String, Vector{Float64}}} function Section(LE_point::Vector{Float64}, TE_point::Vector{Float64}, aero_input::Union{String, Tuple{String, Vector{Float64}}}) @@ -48,12 +48,12 @@ Represents a wing composed of multiple sections with aerodynamic properties. mutable struct Wing n_panels::Int spanwise_panel_distribution::String - spanwise_direction::Vector{Float64} + spanwise_direction::PosVector sections::Vector{Section} function Wing(n_panels::Int; spanwise_panel_distribution::String="linear", - spanwise_direction::Vector{Float64}=[0.0, 1.0, 0.0]) + spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0])) new(n_panels, spanwise_panel_distribution, spanwise_direction, @@ -67,8 +67,8 @@ end Add a new section to the wing. """ -function add_section!(wing::Wing, LE_point::Vector{Float64}, - TE_point::Vector{Float64}, aero_input::Union{String, Tuple{String, Vector{Float64}}}) +function add_section!(wing::Wing, LE_point::PosVector, + TE_point::PosVector, aero_input::Union{String, Tuple{String, Vector{Float64}}}) push!(wing.sections, Section(LE_point, TE_point, aero_input)) end @@ -263,8 +263,8 @@ Returns: function refine_mesh_for_linear_cosine_distribution( spanwise_panel_distribution::String, n_sections::Int, - LE::Union{Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}, - TE::Union{Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}, + LE::Union{Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}, MMatrix{2, 3, Float64, 6}}, + TE::Union{Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}, MMatrix{2, 3, Float64, 6}}, aero_input::Union{Vector{String}, Vector{Tuple{String, Vector{Float64}}}}) # 1. Compute quarter chord line @@ -525,7 +525,7 @@ function calculate_projected_area(wing::Wing, z_plane_vector = z_plane_vector ./ norm(z_plane_vector) # Project point onto plane - function project_onto_plane(point::Vector{Float64}, normal::Vector{Float64}) + function project_onto_plane(point::PosVector, normal::Vector{Float64}) return point .- dot(point, normal) .* normal end diff --git a/test/runtests.jl b/test/runtests.jl index d893ec80..439caef8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,4 +8,5 @@ println("Running tests...") include("test_semi_infinite_filament.jl") include("test_wing_aerodynamics.jl") include("test_wing_geometry.jl") + include("test_plotting.jl") end \ No newline at end of file diff --git a/test/test_plotting.jl b/test/test_plotting.jl new file mode 100644 index 00000000..fd36c427 --- /dev/null +++ b/test/test_plotting.jl @@ -0,0 +1,11 @@ +using VortexStepMethod +using ControlPlots +using Test + +plt.ioff() +@testset "Plotting" begin + fig = plt.plot([1,2,3]) + @test fig isa Vector{plt.PyObject} + show_plot(fig) +end +plt.ion() \ No newline at end of file