1-Bart-1 · Sep 23, 2026 · Sep 22, 2026 · Sep 22, 2026 · Sep 22, 2026 · Sep 22, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -82,6 +82,11 @@
   from its own wing's `spanwise_direction`, not the first wing's. `solve` computes
   `wing_span` and `aspect_ratio_projected` from the extent of all wings along the
   first wing's span, through the new `calculate_span(wings, spanwise_direction)`.
+- `plot_combined_analysis` saves the figure when `is_save` is true, spans its title over
+  the whole figure, frames the wing geometry, and draws literature polars dashed as
+  `plot_polars` does.
+- `plot_polars` and `plot_combined_analysis` draw a CS that is zero up to round-off as a
+  flat line instead of scaling the axis to the round-off noise.
 - `perpendicular_sections` and `obj_to_yaml` orient each cut by the leading edge through
   the neighbouring stations, not by the step that reached the station, which near a
   closing tip grazed the surface and returned an 80%-thick section. Every section moves

diff --git a/README.md b/README.md
@@ -16,7 +16,9 @@ This package was translated from the Python code version 1.0.0 available at http
 
 Example output is shown in the figure below, where the aerodynamic predictions of the TU Delft V3 LEI kite obtained with the Julia and Python implementations of the VSM are compared against CFD simulations and wind tunnel measurements.
 
-![Example output](docs/v3_example_output.png)
+![Example output](docs/TU_Delft_V3_Kite.png)
+
+The figure is rendered by [examples/readme_figure.jl](examples/readme_figure.jl).
 
 ## Installation
 Install [Julia 1.12](https://ufechner7.github.io/2024/08/09/installing-julia-with-juliaup.html) or later, 

diff --git a/docs/TU_Delft_V3_Kite.png b/docs/TU_Delft_V3_Kite.png
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
@@ -282,6 +282,8 @@ span_axis
 create_geometry_plot_makie
 plot_line_segment_makie!
 set_axes_equal_makie!
+polar_style
+widen_flat_ylims!
 map_airfoil_3d
 fitted_airfoil_3d
 generated_slices

diff --git a/docs/v3_example_output.png b/docs/v3_example_output.png
diff --git a/examples/V3_kite.jl b/examples/V3_kite.jl
@@ -16,34 +16,10 @@ NEURALFOIL = true
 # Rolling-ball radius of the shrink wrap; fillets the concave tube-canopy junction.
 MIN_CONCAVE_RADIUS = 0.4
 OUTPUT_DIR = joinpath(dirname(@__DIR__), "output")
-REFERENCE_POINT = [0.422646, 0.0, 9.3667]
 
-project_dir = dirname(@__DIR__)
-literature_paths = [
-    joinpath(project_dir, "data", "TUDELFT_V3_KITE", "literature_results",
-        "CFD_RANS_Rey_5e5_Poland2025_alpha_sweep_beta_0_NoStruts.csv"),
-    joinpath(project_dir, "data", "TUDELFT_V3_KITE", "literature_results",
-        "CFD_RANS_Rey_10e5_Poland2025_alpha_sweep_beta_0.csv"),
-    joinpath(project_dir, "data", "TUDELFT_V3_KITE", "literature_results",
-        "python_alpha_sweep.csv"),
-    joinpath(project_dir, "data", "TUDELFT_V3_KITE", "literature_results",
-        "windtunnel_alpha_sweep_beta_00_0_Poland_2025_Rey_5e5.csv"),
-]
-beta_literature_paths = [
-    joinpath(project_dir, "data", "TUDELFT_V3_KITE", "literature_results",
-        "windtunnel_beta_sweep_alpha_07_4_Poland_2025_Rey_5e5.csv"),
-]
-
-settings = VSMSettings(joinpath(project_dir, "data", "TUDELFT_V3_KITE",
-    "vsm_settings.yaml"); data_prefix=false)
-settings.wings[1].geometry_file = joinpath(project_dir,
-    settings.wings[1].geometry_file)
-wing = Wing(settings)
-refine!(wing)
-body_aero = BodyAerodynamics([wing])
-VortexStepMethod.reinit!(body_aero)
-solver = Solver(settings)
-solver.reference_point .= REFERENCE_POINT
+include("V3_kite_setup.jl")
+beta_literature_paths = [joinpath(literature_dir,
+    "windtunnel_beta_sweep_alpha_07_4_Poland_2025_Rey_5e5.csv")]
 
 if DEFORM
     VortexStepMethod.unrefined_deform!(

diff --git a/examples/V3_kite_setup.jl b/examples/V3_kite_setup.jl
@@ -0,0 +1,24 @@
+# TU Delft V3 kite from vsm_settings.yaml, shared by V3_kite.jl and readme_figure.jl.
+using VortexStepMethod
+
+REFERENCE_POINT = [0.422646, 0.0, 9.3667]
+
+project_dir = dirname(@__DIR__)
+literature_dir = joinpath(project_dir, "data", "TUDELFT_V3_KITE", "literature_results")
+literature_paths = joinpath.(literature_dir, [
+    "CFD_RANS_Rey_5e5_Poland2025_alpha_sweep_beta_0_NoStruts.csv",
+    "CFD_RANS_Rey_10e5_Poland2025_alpha_sweep_beta_0.csv",
+    "python_alpha_sweep.csv",
+    "windtunnel_alpha_sweep_beta_00_0_Poland_2025_Rey_5e5.csv",
+])
+
+settings = VSMSettings(joinpath(project_dir, "data", "TUDELFT_V3_KITE",
+    "vsm_settings.yaml"); data_prefix=false)
+settings.wings[1].geometry_file = joinpath(project_dir,
+    settings.wings[1].geometry_file)
+wing = Wing(settings)
+refine!(wing)
+body_aero = BodyAerodynamics([wing])
+VortexStepMethod.reinit!(body_aero)
+solver = Solver(settings)
+solver.reference_point .= REFERENCE_POINT
diff --git a/examples/readme_figure.jl b/examples/readme_figure.jl
@@ -0,0 +1,29 @@
+# Renders docs/TU_Delft_V3_Kite.png, the V3 kite figure shown in the README.
+using Pkg
+if Base.active_project() != joinpath(@__DIR__, "Project.toml")
+    Pkg.activate(@__DIR__)
+end
+using CairoMakie
+using MakieControlPlots
+CairoMakie.activate!()
+using VortexStepMethod
+
+include("V3_kite_setup.jl")
+labels = ["VSM Julia", "CFD Re=5e5", "CFD Re=10e5", "VSM Python Re=5e5",
+    "Wind tunnel Re=5e5"]
+
+set_va!(body_aero, settings)
+results = VortexStepMethod.solve(solver, body_aero)
+
+plot_combined_analysis(solver, body_aero, results;
+    labels,
+    literature_path_list=literature_paths,
+    angle_range=range(-5, 25, length=31),
+    va=settings.condition.va,
+    angle_of_attack_for_spanwise_distribution=10.0,
+    title="TU Delft V3 Kite",
+    is_show=false,
+    is_save=true,
+    save_path=joinpath(project_dir, "docs"),
+    data_type=".png",
+)
diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl
@@ -589,6 +589,31 @@ function set_axes_equal_makie!(ax, panels; zoom=1.8)
         z_mid - max_range / 2, z_mid + max_range / 2)
 end
 
+"""
+    polar_style(idx, n_solvers)
+
+Marker, marker size and line style of the `idx`-th polar: solver output for
+`idx <= n_solvers`, literature otherwise.
+"""
+polar_style(idx, n_solvers) = idx <= n_solvers ?
+    (marker=:star5, markersize=12, linestyle=:solid) :
+    (marker=:circle, markersize=8, linestyle=:dash)
+
+"""
+    widen_flat_ylims!(ax, series; min_span=0.1)
+
+Set the y-limits of `ax` to span at least `min_span` around the non-NaN values of
+all vectors in `series`.
+"""
+function widen_flat_ylims!(ax, series; min_span=0.1)
+    values = filter(!isnan, reduce(vcat, series))
+    isempty(values) && return
+    low, high = extrema(values)
+    high - low >= min_span && return
+    mid = (low + high) / 2
+    ylims!(ax, mid - min_span / 2, mid + min_span / 2)
+end
+
 """
     create_geometry_plot_makie(body_aero::BodyAerodynamics, title,
                                view_elevation, view_azimuth; zoom=1.8)
@@ -960,35 +985,22 @@ function VortexStepMethod.plot_polars(
 
         for (i, (polar_data, cm, label)) in enumerate(
                 zip(polar_data_list, cm_data_list, labels_with_re))
-            marker = i <= n_solvers ? :star5 : :circle
-            markersize = i <= n_solvers ? 12 : 8
-            linestyle = i <= n_solvers ? :solid : :dash
+            style = polar_style(i, n_solvers)
             angles = polar_data[1]
-            scatterlines!(ax_cl, angles, polar_data[2];
-                label=label, marker=marker, markersize=markersize, linestyle)
-            scatterlines!(ax_cd, angles, polar_data[3];
-                label=label, marker=marker, markersize=markersize, linestyle)
-            scatterlines!(ax_cs, angles, polar_data[4];
-                label=label, marker=marker, markersize=markersize, linestyle)
+            scatterlines!(ax_cl, angles, polar_data[2]; label, style...)
+            scatterlines!(ax_cd, angles, polar_data[3]; label, style...)
+            scatterlines!(ax_cs, angles, polar_data[4]; label, style...)
             if !all(isnan, cm.cmx)
-                scatterlines!(ax_cmx, angles,
-                    Float64.(cm.cmx);
-                    label=label, marker=marker,
-                    markersize=markersize, linestyle)
+                scatterlines!(ax_cmx, angles, Float64.(cm.cmx); label, style...)
             end
             if !all(isnan, cm.cmy)
-                scatterlines!(ax_cmy, angles,
-                    Float64.(cm.cmy);
-                    label=label, marker=marker,
-                    markersize=markersize, linestyle)
+                scatterlines!(ax_cmy, angles, Float64.(cm.cmy); label, style...)
             end
             if !all(isnan, cm.cmz)
-                scatterlines!(ax_cmz, angles,
-                    Float64.(cm.cmz);
-                    label=label, marker=marker,
-                    markersize=markersize, linestyle)
+                scatterlines!(ax_cmz, angles, Float64.(cm.cmz); label, style...)
             end
         end
+        widen_flat_ylims!(ax_cs, [pd[4] for pd in polar_data_list])
         Legend(fig[3, 1:3], ax_cl;
             orientation=:horizontal, tellwidth=false,
             tellheight=true)
@@ -1011,30 +1023,18 @@ function VortexStepMethod.plot_polars(
 
         for (i, (polar_data, label)) in enumerate(
                 zip(polar_data_list, labels_with_re))
-            marker = i <= n_solvers ? :star5 : :circle
-            markersize = i <= n_solvers ? 12 : 8
-            linestyle = i <= n_solvers ? :solid : :dash
-            scatterlines!(ax_cl, polar_data[1], polar_data[2];
-                label=label, marker=marker,
-                markersize=markersize, linestyle)
-            scatterlines!(ax_cd, polar_data[1], polar_data[3];
-                label=label, marker=marker,
-                markersize=markersize, linestyle)
-            scatterlines!(ax_cs, polar_data[1], polar_data[4];
-                label=label, marker=marker,
-                markersize=markersize, linestyle)
+            style = polar_style(i, n_solvers)
+            scatterlines!(ax_cl, polar_data[1], polar_data[2]; label, style...)
+            scatterlines!(ax_cd, polar_data[1], polar_data[3]; label, style...)
+            scatterlines!(ax_cs, polar_data[1], polar_data[4]; label, style...)
             if cl_over_cd
                 cl_cd = polar_data[2] ./ polar_data[3]
-                scatterlines!(ax_fourth, polar_data[1], cl_cd;
-                    label=label, marker=marker,
-                    markersize=markersize, linestyle)
+                scatterlines!(ax_fourth, polar_data[1], cl_cd; label, style...)
             else
-                scatterlines!(ax_fourth, polar_data[3],
-                    polar_data[2];
-                    label=label, marker=marker,
-                    markersize=markersize, linestyle)
+                scatterlines!(ax_fourth, polar_data[3], polar_data[2]; label, style...)
             end
         end
+        widen_flat_ylims!(ax_cs, [pd[4] for pd in polar_data_list])
         Legend(fig[3, :], ax_cl;
             orientation=:horizontal, tellwidth=false,
             tellheight=true)
@@ -1215,7 +1215,7 @@ function VortexStepMethod.plot_combined_analysis(
         # Fallback if screen detection fails
         Figure(size=(1800, 1200))
     end
-    Label(fig[0, :], title, fontsize=20, font=:bold)
+    Label(fig[0, 1:2], title, fontsize=20, font=:bold)
 
     # Use first body_aero for geometry and polar data display
     first_body = body_aeros[1]
@@ -1280,7 +1280,7 @@ function VortexStepMethod.plot_combined_analysis(
     plot_line_segment_makie!(ax_geo, [va_vector_begin, va_vector_end],
         :lightblue, "va")
 
-    set_axes_equal_makie!(ax_geo, panels; zoom=0.5)
+    set_axes_equal_makie!(ax_geo, panels; zoom=1.0)
     axislegend(ax_geo; position=:lt)
 
     # [1,2] Polar Data Surfaces or Curves
@@ -1363,34 +1363,18 @@ function VortexStepMethod.plot_combined_analysis(
         color = colors[mod1(si, length(colors))]
         y_si = [panel.aero_center[2] for panel in body_aeros[si].panels]
 
-        cl_val = round(rs["cl"], digits=2)
-        lines!(ax_cl, Vector(y_si), Vector(rs["cl_distribution"]);
-            label="$lbl CL: $cl_val", color)
-
-        cd_val = round(rs["cd"], digits=2)
-        lines!(ax_cd, Vector(y_si), Vector(rs["cd_distribution"]);
-            label="$lbl CD: $cd_val", color)
-
+        lines!(ax_cl, Vector(y_si), Vector(rs["cl_distribution"]); color)
+        lines!(ax_cd, Vector(y_si), Vector(rs["cd_distribution"]); color)
         lines!(ax_gamma, Vector(y_si),
             Vector(rs["gamma_distribution"]); label=lbl, color)
-
         lines!(ax_alpha_geo, Vector(y_si),
-            rad2deg.(Vector(rs["alpha_geometric"])); label=lbl, color)
-
+            rad2deg.(Vector(rs["alpha_geometric"])); color)
         lines!(ax_alpha_ac, Vector(y_si),
-            rad2deg.(Vector(rs["alpha_at_ac"])); label=lbl, color)
-
+            rad2deg.(Vector(rs["alpha_at_ac"])); color)
         lines!(ax_alpha_unc, Vector(y_si),
-            rad2deg.(Vector(rs["alpha_uncorrected"]));
-            label=lbl, color)
-
-        force_axes = [ax_fx, ax_fy, ax_fz]
-        components = ["x", "y", "z"]
-        for (idx, (ax, comp)) in enumerate(zip(force_axes, components))
-            forces = rs["F_distribution"][idx, :]
-            total_force = round(rs["F$comp"], digits=2)
-            lines!(ax, Vector(y_si), Vector(forces);
-                label="$lbl ΣF$comp: $total_force N", color)
+            rad2deg.(Vector(rs["alpha_uncorrected"])); color)
+        for (idx, ax) in enumerate((ax_fx, ax_fy, ax_fz))
+            lines!(ax, Vector(y_si), Vector(rs["F_distribution"][idx, :]); color)
         end
     end
 
@@ -1448,21 +1432,24 @@ function VortexStepMethod.plot_combined_analysis(
 
     for (idx, (pd, lbl)) in enumerate(polar_series)
         color = colors[mod1(idx, length(colors))]
-        marker = idx == 1 ? :star5 : :circle
-        markersize = idx == 1 ? 12 : 8
-
-        scatterlines!(ax_cl_polar, pd.angle, pd.cl; label=lbl, marker, markersize, color)
-        scatterlines!(ax_cd_polar, pd.angle, pd.cd; label=lbl, marker, markersize, color)
-        scatterlines!(ax_cs_polar, pd.angle, pd.cs; label=lbl, marker, markersize, color)
+        style = polar_style(idx, n_solvers)
+
+        scatterlines!(ax_cl_polar, pd.angle, pd.cl;
+            label=lbl, color, style...)
+        scatterlines!(ax_cd_polar, pd.angle, pd.cd;
+            label=lbl, color, style...)
+        scatterlines!(ax_cs_polar, pd.angle, pd.cs;
+            label=lbl, color, style...)
         if cl_over_cd
             cl_cd = pd.cl ./ pd.cd
             scatterlines!(ax_fourth_polar, pd.angle, cl_cd;
-                label=lbl, marker, markersize, color)
+                label=lbl, color, style...)
         else
             scatterlines!(ax_fourth_polar, pd.cd, pd.cl;
-                label=lbl, marker, markersize, color)
+                label=lbl, color, style...)
         end
     end
+    widen_flat_ylims!(ax_cs_polar, [pd.cs for (pd, _) in polar_series])
     Legend(fig[2, 2][3, :], ax_cl_polar;
         orientation=:horizontal, nbanks=2,
         tellwidth=false, tellheight=true)
@@ -1471,6 +1458,7 @@ function VortexStepMethod.plot_combined_analysis(
     colsize!(fig.layout, 1, Relative(0.6))
     colsize!(fig.layout, 2, Relative(0.4))
 
+    is_save && save_plot(fig, save_path, title; data_type)
     is_show && show_plot(fig; name=title)
 
     return fig

diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl
@@ -143,11 +143,14 @@ end
         angle_of_attack=30.0,
         va=va,
         title="Combined Analysis",
-        is_save=false,
+        data_type=".png",
+        save_path=save_dir,
+        is_save=true,
         is_show=false,
         cl_over_cd=true
     )
     @test fig isa Figure
+    @test isfile(joinpath(save_dir, "Combined_Analysis.png"))
 
     # Plot combined analysis with cl_over_cd=false
     fig = plot_combined_analysis(
@@ -163,6 +166,16 @@ end
     )
     @test fig isa Figure
 
+    @testset "round-off CS plots flat, a varying CS keeps its autoscale" begin
+        fig = Figure()
+        ax_noise = Axis(fig[1, 1])
+        makie_ext.widen_flat_ylims!(ax_noise, [[1e-16, -1e-16], [0.0, NaN]])
+        @test ax_noise.limits[][2] == (-0.05, 0.05)
+        ax_wide = Axis(fig[1, 2])
+        makie_ext.widen_flat_ylims!(ax_wide, [[0.0, 0.5]])
+        @test isnothing(ax_wide.limits[][2])
+    end
+
     # Test polar data plotting
     body_aero = BodyAerodynamics([ram_wing])
     fig = plot_polar_data(body_aero; is_show=false)