From 01ab50e167c84d4c0a94c737958a3f220b4812f0 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 25 Aug 2025 06:18:30 +0200
Subject: [PATCH 01/10] add Makie extension

---
 lib/ControlSystemsBase/Project.toml           |   7 +-
 .../ext/ControlSystemsBaseMakieExt.jl         | 828 ++++++++++++++++++
 lib/ControlSystemsBase/src/CSMakie.jl         | 164 ++++
 .../src/ControlSystemsBase.jl                 |   4 +
 .../test/test_makie_plots.jl                  | 160 ++++
 5 files changed, 1162 insertions(+), 1 deletion(-)
 create mode 100644 lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
 create mode 100644 lib/ControlSystemsBase/src/CSMakie.jl
 create mode 100644 lib/ControlSystemsBase/test/test_makie_plots.jl

diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml
index 3f4a4f32a..912034afa 100644
--- a/lib/ControlSystemsBase/Project.toml
+++ b/lib/ControlSystemsBase/Project.toml
@@ -23,10 +23,12 @@ UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 ImplicitDifferentiation = "57b37032-215b-411a-8a7c-41a003a55207"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 
 [extensions]
 ControlSystemsBaseDSPExt = ["DSP"]
 ControlSystemsBaseImplicitDifferentiationExt = ["ImplicitDifferentiation", "ComponentArrays"]
+ControlSystemsBaseMakieExt = ["Makie"]
 
 [compat]
 Aqua = "0.5"
@@ -37,6 +39,7 @@ Hungarian = "0.7.0"
 ImplicitDifferentiation = "0.7.2"
 LinearAlgebra = "<0.0.1, 1"
 MacroTools = "0.5"
+Makie = "0.24"
 MatrixEquations = "1, 2.1"
 MatrixPencils = "1.8.3"
 Polynomials = "3.0, 4.0"
@@ -54,13 +57,15 @@ ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
 ImplicitDifferentiation = "57b37032-215b-411a-8a7c-41a003a55207"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Aqua", "ComponentArrays", "Documenter", "DSP", "FiniteDifferences", "ImplicitDifferentiation", "GenericSchur", "GR", "Plots", "SparseArrays", "StaticArrays"]
+test = ["Test", "Aqua", "ComponentArrays", "Documenter", "DSP", "FiniteDifferences", "ImplicitDifferentiation", "Makie", "GLMakie", "GenericSchur", "GR", "Plots", "SparseArrays", "StaticArrays"]
diff --git a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
new file mode 100644
index 000000000..8742af43e
--- /dev/null
+++ b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
@@ -0,0 +1,828 @@
+module ControlSystemsBaseMakieExt
+
+using ControlSystemsBase
+using ControlSystemsBase.CSMakie
+using Makie
+using LinearAlgebra
+using Printf
+
+# Import necessary functions from ControlSystemsBase
+using ControlSystemsBase: downsample, _processfreqplot, _default_freq_vector,
+                          _same_io_dims, _get_plotlabel, _to1series,
+                          SimResult, StepInfo, RootLocusResult,
+                          poles, tzeros, bode, nyquist, sigma, margin, 
+                          sisomargin, relative_gain_array, rlocus,
+                          input_names, output_names, state_names, system_name,
+                          iscontinuous, isdiscrete, issiso, isrational,
+                          integrator_excess, balance_statespace, LTISystem
+
+# Global plot scale settings (matching Plots.jl version)
+const _PlotScale = Ref("log10")
+const _PlotScaleFunc = Ref(:log10)
+const _PlotScaleStr = Ref("")
+
+"""
+    CSMakie.setPlotScale(str)
+
+Set the default scale of magnitude in `bodeplot` and `sigmaplot` for Makie plots.
+`str` should be either `"dB"` or `"log10"`. The default scale if none is chosen is `"log10"`.
+"""
+function CSMakie.setPlotScale(str::AbstractString)
+    if str == "dB"
+        _PlotScale[] = str
+        _PlotScaleFunc[] = :identity
+        _PlotScaleStr[] = "(dB)"
+    elseif str == "log10"
+        _PlotScale[] = str
+        _PlotScaleFunc[] = :log10
+        _PlotScaleStr[] = ""
+    else
+        error("Scale must be set to either \"dB\" or \"log10\"")
+    end
+end
+
+# Helper function to get y-scale transform
+function get_yscale_transform()
+    _PlotScaleFunc[] == :log10 ? Makie.log10 : identity
+end
+
+# ====== PZMap (Pole-Zero Map) ======
+function CSMakie.pzmap(systems::Union{LTISystem, AbstractVector{<:LTISystem}}; 
+                       hz=false, kwargs...)
+    systems_vec = systems isa AbstractVector ? systems : [systems]
+    
+    fig = Figure()
+    ax = Axis(fig[1,1], aspect = DataAspect(),
+             title = "Pole-zero map",
+             xlabel = hz ? "Real [Hz]" : "Real",
+             ylabel = hz ? "Imag [Hz]" : "Imag")
+    
+    # Add grid at zero
+    vlines!(ax, 0, color=:gray, alpha=0.3, linewidth=0.5)
+    hlines!(ax, 0, color=:gray, alpha=0.3, linewidth=0.5)
+    
+    scale_factor = hz ? 1/(2π) : 1
+    
+    for (i, system) in enumerate(systems_vec)
+        p = poles(system) .* scale_factor
+        z = tzeros(system) .* scale_factor
+        
+        # Plot zeros as circles
+        if !isempty(z)
+            scatter!(ax, real(z), imag(z), 
+                    marker=:circle, markersize=15, 
+                    color=(:transparent, 0.5), strokewidth=2, 
+                    strokecolor=Cycled(i), label="Zeros")
+        end
+        
+        # Plot poles as x's
+        if !isempty(p)
+            scatter!(ax, real(p), imag(p), 
+                    marker=:xcross, markersize=15,
+                    color=Cycled(i), label="Poles")
+        end
+        
+        # Add unit circle for discrete systems
+        if isdiscrete(system) && i == 1  # Only draw once
+            θ = range(0, 2π, length=100)
+            circle_scale = scale_factor
+            lines!(ax, circle_scale .* cos.(θ), circle_scale .* sin.(θ), 
+                  color=:gray, linestyle=:dash, alpha=0.5)
+        end
+    end
+    
+    return fig
+end
+
+function CSMakie.pzmap!(ax::Axis, systems::Union{LTISystem, AbstractVector{<:LTISystem}}; 
+                        hz=false, kwargs...)
+    systems_vec = systems isa AbstractVector ? systems : [systems]
+    
+    # Add grid at zero
+    vlines!(ax, 0, color=:gray, alpha=0.3, linewidth=0.5)
+    hlines!(ax, 0, color=:gray, alpha=0.3, linewidth=0.5)
+    
+    scale_factor = hz ? 1/(2π) : 1
+    
+    for (i, system) in enumerate(systems_vec)
+        p = poles(system) .* scale_factor
+        z = tzeros(system) .* scale_factor
+        
+        # Plot zeros as circles
+        if !isempty(z)
+            scatter!(ax, real(z), imag(z), 
+                    marker=:circle, markersize=15, 
+                    color=(:transparent, 0.5), strokewidth=2, 
+                    strokecolor=Cycled(i))
+        end
+        
+        # Plot poles as x's
+        if !isempty(p)
+            scatter!(ax, real(p), imag(p), 
+                    marker=:xcross, markersize=15,
+                    color=Cycled(i))
+        end
+        
+        # Add unit circle for discrete systems
+        if isdiscrete(system) && i == 1  # Only draw once
+            θ = range(0, 2π, length=100)
+            circle_scale = scale_factor
+            lines!(ax, circle_scale .* cos.(θ), circle_scale .* sin.(θ), 
+                  color=:gray, linestyle=:dash, alpha=0.5)
+        end
+    end
+    
+    return ax
+end
+
+# ====== Bodeplot ======
+function CSMakie.bodeplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+                          w=nothing; plotphase=true, unwrap=true, hz=false, 
+                          balance=true, adjust_phase_start=true, adaptive=true, kwargs...)
+    systems_vec = systems isa AbstractVector ? systems : [systems]
+    systems, w = isnothing(w) ? _processfreqplot(Val{:bode}(), systems_vec; adaptive) : 
+                                _processfreqplot(Val{:bode}(), systems_vec, w; adaptive)
+    
+    ws = (hz ? 1/(2π) : 1) .* w
+    ny, nu = size(systems[1])
+    
+    fig = Figure()
+    gl = GridLayout(fig[1, 1])
+    
+    # Create axes grid
+    n_rows = plotphase ? 2 : 1
+    axes_mag = Matrix{Any}(undef, ny, nu)
+    axes_phase = plotphase ? Matrix{Any}(undef, ny, nu) : nothing
+    
+    for j in 1:nu
+        for i in 1:ny
+            row_mag = plotphase ? 2i - 1 : i
+            
+            # Magnitude axis
+            ax_mag = Axis(gl[row_mag, j],
+                         xscale = log10,
+                         yscale = get_yscale_transform(),
+                         xlabel = plotphase ? "" : (hz ? "Frequency [Hz]" : "Frequency [rad/s]"),
+                         ylabel = i == 1 ? "Magnitude $_PlotScaleStr[]" : "",
+                         title = j == 1 && i == 1 ? "Bode Plot" : "")
+            axes_mag[i, j] = ax_mag
+            
+            # Phase axis
+            if plotphase
+                ax_phase = Axis(gl[2i, j],
+                               xscale = log10,
+                               xlabel = hz ? "Frequency [Hz]" : "Frequency [rad/s]",
+                               ylabel = i == 1 ? "Phase (deg)" : "")
+                axes_phase[i, j] = ax_phase
+                linkxaxes!(ax_mag, ax_phase)
+            end
+        end
+    end
+    
+    # Plot data for each system
+    for (si, s) in enumerate(systems)
+        sbal = balance ? balance_statespace(s)[1] : s
+        
+        intexcess = 0
+        if plotphase && adjust_phase_start && isrational(sbal)
+            intexcess = integrator_excess(sbal)
+        end
+        
+        mag, phase = bode(sbal, w; unwrap=false)
+        
+        if _PlotScale[] == "dB"
+            mag = 20*log10.(mag)
+        elseif 0 ∈ mag
+            replace!(mag, 0 => NaN)
+        end
+        
+        for j in 1:nu
+            for i in 1:ny
+                magdata = vec(mag[i, j, :])
+                if all(isnan.(magdata) .| isinf.(magdata))
+                    continue
+                end
+                
+                # Plot magnitude
+                ax_mag = axes_mag[i, j]
+                lab = _get_plotlabel(s, i, j)
+                
+                if adaptive && eltype(magdata) <: AbstractFloat
+                    lmag = _PlotScale[] == "dB" ? magdata : log.(magdata)
+                    wsi, magi, _ = downsample(ws, lmag, maximum(lmag) - minimum(lmag)/500)
+                    lines!(ax_mag, wsi, magi, label=lab)
+                else
+                    lines!(ax_mag, ws, magdata, label=lab)
+                end
+                
+                # Plot phase
+                if plotphase
+                    phasedata = vec(phase[i, j, :])
+                    
+                    if adjust_phase_start && isrational(sbal) && intexcess != 0
+                        nineties = round(Int, phasedata[1] / 90)
+                        phasedata .+= ((90*(-intexcess-nineties)) ÷ 360) * 360
+                    end
+                    
+                    if unwrap
+                        phasedata = ControlSystemsBase.unwrap(phasedata .* (π/180)) .* (180/π)
+                    end
+                    
+                    ax_phase = axes_phase[i, j]
+                    if adaptive && eltype(phasedata) <: AbstractFloat
+                        wsp, phasep, _ = downsample(ws, phasedata, 
+                                                   maximum(phasedata) - minimum(phasedata)/500)
+                        lines!(ax_phase, wsp, phasep)
+                    else
+                        lines!(ax_phase, ws, phasedata)
+                    end
+                end
+            end
+        end
+    end
+    
+    return fig
+end
+
+# ====== Nyquistplot ======
+function CSMakie.nyquistplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+                             w=nothing; Ms_circles=Float64[], Mt_circles=Float64[], 
+                             unit_circle=false, hz=false, critical_point=-1, 
+                             balance=true, adaptive=true, kwargs...)
+    systems_vec = systems isa AbstractVector ? systems : [systems]
+    systems, w = isnothing(w) ? _processfreqplot(Val{:nyquist}(), systems_vec; adaptive) : 
+                                _processfreqplot(Val{:nyquist}(), systems_vec, w; adaptive)
+    
+    ny, nu = size(systems[1])
+    
+    fig = Figure()
+    gl = GridLayout(fig[1, 1])
+    
+    # Create axes grid
+    axes = Matrix{Any}(undef, ny, nu)
+    
+    for j in 1:nu
+        for i in 1:ny
+            ax = Axis(gl[i, j],
+                     aspect = DataAspect(),
+                     xlabel = j == ny ? "Real" : "",
+                     ylabel = i == 1 ? "Imaginary" : "",
+                     title = i == 1 && j == 1 ? "Nyquist Plot" : "")
+            axes[i, j] = ax
+            
+            # Add grid lines at zero
+            vlines!(ax, 0, color=:gray, alpha=0.3, linewidth=0.5)
+            hlines!(ax, 0, color=:gray, alpha=0.3, linewidth=0.5)
+        end
+    end
+    
+    θ = range(0, 2π, length=100)
+    
+    for (si, s) in enumerate(systems)
+        sbal = balance ? balance_statespace(s)[1] : s
+        re_resp, im_resp = nyquist(sbal, w)[1:2]
+        
+        for j in 1:nu
+            for i in 1:ny
+                redata = vec(re_resp[i, j, :])
+                imdata = vec(im_resp[i, j, :])
+                
+                ax = axes[i, j]
+                lab = _get_plotlabel(s, i, j)
+                
+                if adaptive
+                    indsre, _, _ = downsample(w, redata, 1/500)
+                    indsim, _, _ = downsample(w, imdata, 1/500)
+                    inds = sort!(union(indsre, indsim))
+                    lines!(ax, redata[inds], imdata[inds], label=lab)
+                else
+                    lines!(ax, redata, imdata, label=lab)
+                end
+                
+                # Add circles and critical point for last system
+                if si == length(systems)
+                    # Ms circles
+                    for Ms in Ms_circles
+                        lines!(ax, -1 .+ (1/Ms) .* cos.(θ), (1/Ms) .* sin.(θ),
+                              color=:gray, linestyle=:dash, alpha=0.5,
+                              label="Ms = $(round(Ms, digits=2))")
+                    end
+                    
+                    # Mt circles
+                    for Mt in Mt_circles
+                        ct = -Mt^2/(Mt^2-1)
+                        rt = Mt/(abs(Mt^2-1))
+                        lines!(ax, ct .+ rt .* cos.(θ), rt .* sin.(θ),
+                              color=:gray, linestyle=:dash, alpha=0.5,
+                              label="Mt = $(round(Mt, digits=2))")
+                    end
+                    
+                    # Unit circle
+                    if unit_circle
+                        lines!(ax, cos.(θ), sin.(θ),
+                              color=:gray, linestyle=:dash, alpha=0.5)
+                    end
+                    
+                    # Critical point
+                    scatter!(ax, [critical_point], [0],
+                            marker=:xcross, markersize=15, color=:red)
+                end
+            end
+        end
+    end
+    
+    return fig
+end
+
+# ====== Sigmaplot ======
+function CSMakie.sigmaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+                           w=nothing; hz=false, balance=true, extrema=false, kwargs...)
+    systems_vec = systems isa AbstractVector ? systems : [systems]
+    systems, w = isnothing(w) ? _processfreqplot(Val{:sigma}(), systems_vec) : 
+                                _processfreqplot(Val{:sigma}(), systems_vec, w)
+    
+    ws = (hz ? 1/(2π) : 1) .* w
+    
+    fig = Figure()
+    ax = Axis(fig[1,1],
+             xscale = log10,
+             yscale = get_yscale_transform(),
+             title = "Sigma Plot",
+             xlabel = hz ? "Frequency [Hz]" : "Frequency [rad/s]",
+             ylabel = "Singular Values $_PlotScaleStr[]")
+    
+    for (si, s) in enumerate(systems)
+        sbal = balance ? balance_statespace(s)[1] : s
+        sv = sigma(sbal, w)[1]'
+        
+        if extrema && size(sv, 2) > 2
+            sv = sv[:, [1, end]]
+        end
+        
+        if _PlotScale[] == "dB"
+            sv = 20*log10.(sv)
+        end
+        
+        # Use _to1series to efficiently plot multiple lines
+        ws_series, sv_series = _to1series(ws, sv)
+        lines!(ax, ws_series, sv_series, color=Cycled(si))
+    end
+    
+    return fig
+end
+
+# ====== Marginplot ======
+function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+                            w=nothing; plotphase=true, hz=false, balance=true, 
+                            adjust_phase_start=true, adaptive=true, kwargs...)
+    systems_vec = systems isa AbstractVector ? systems : [systems]
+    systems, w = isnothing(w) ? _processfreqplot(Val{:bode}(), systems_vec; adaptive) : 
+                                _processfreqplot(Val{:bode}(), systems_vec, w; adaptive)
+    
+    ws = (hz ? 1/(2π) : 1) .* w
+    ny, nu = size(systems[1])
+    
+    fig = Figure()
+    gl = GridLayout(fig[1, 1])
+    
+    # Create axes grid
+    n_rows = plotphase ? 2 : 1
+    axes_mag = Matrix{Any}(undef, ny, nu)
+    axes_phase = plotphase ? Matrix{Any}(undef, ny, nu) : nothing
+    
+    for j in 1:nu
+        for i in 1:ny
+            row_mag = plotphase ? 2i - 1 : i
+            
+            # Magnitude axis
+            ax_mag = Axis(gl[row_mag, j],
+                         xscale = log10,
+                         yscale = get_yscale_transform(),
+                         xlabel = plotphase ? "" : (hz ? "Frequency [Hz]" : "Frequency [rad/s]"),
+                         ylabel = i == 1 ? "Magnitude $_PlotScaleStr[]" : "")
+            axes_mag[i, j] = ax_mag
+            
+            # Phase axis
+            if plotphase
+                ax_phase = Axis(gl[2i, j],
+                               xscale = log10,
+                               xlabel = hz ? "Frequency [Hz]" : "Frequency [rad/s]",
+                               ylabel = i == 1 ? "Phase (deg)" : "")
+                axes_phase[i, j] = ax_phase
+                linkxaxes!(ax_mag, ax_phase)
+            end
+        end
+    end
+    
+    # Plot data for each system
+    for (si, s) in enumerate(systems)
+        sbal = balance ? balance_statespace(s)[1] : s
+        bmag, bphase = bode(sbal, w)
+        
+        for j in 1:nu
+            for i in 1:ny
+                wgm, gm, wpm, pm, fullPhase, phasedata = sisomargin(sbal[i,j], w; 
+                                                                     full=true, 
+                                                                     allMargins=true, 
+                                                                     adjust_phase_start)
+                
+                # Limit number of margins shown
+                if length(gm) > 5
+                    idx = sortperm(gm)
+                    wgm = wgm[idx[1:5]]
+                    gm = gm[idx[1:5]]
+                end
+                if length(pm) > 5
+                    idx = sortperm(pm)
+                    wpm = wpm[idx[1:5]]
+                    pm = pm[idx[1:5]]
+                end
+                
+                # Magnitude plot
+                ax_mag = axes_mag[i, j]
+                magdata = vec(bmag[i, j, :])
+                
+                if _PlotScale[] == "dB"
+                    magdata = 20*log10.(magdata)
+                    mag_margins = 20 .* log10.(1 ./ gm)
+                    oneLine = 0
+                else
+                    mag_margins = 1 ./ gm
+                    oneLine = 1
+                end
+                
+                # Plot magnitude response
+                if adaptive
+                    wsi, magi, _ = downsample(ws, magdata, 
+                                            maximum(magdata) - minimum(magdata)/500)
+                    lines!(ax_mag, wsi, magi, label=_get_plotlabel(s, i, j))
+                else
+                    lines!(ax_mag, ws, magdata, label=_get_plotlabel(s, i, j))
+                end
+                
+                # Unity gain line
+                hlines!(ax_mag, oneLine, color=:gray, linestyle=:dash, alpha=0.5)
+                
+                # Gain margins
+                wgm_display = hz ? wgm ./ (2π) : wgm
+                for k in 1:length(gm)
+                    lines!(ax_mag, [wgm_display[k], wgm_display[k]], 
+                          [oneLine, mag_margins[k]], color=:red)
+                    scatter!(ax_mag, [wgm_display[k]], [mag_margins[k]], 
+                            color=:red, markersize=8)
+                end
+                
+                # Add title with margins info
+                if si == length(systems)
+                    gm_str = "Gm: [" * join([Printf.@sprintf("%.2g", g) for g in gm], ", ") * "]"
+                    wgm_str = "ωgm: [" * join([Printf.@sprintf("%.2g", w) for w in wgm_display], ", ") * "]"
+                    ax_mag.title = gm_str * " " * wgm_str
+                end
+                
+                # Phase plot
+                if plotphase
+                    ax_phase = axes_phase[i, j]
+                    
+                    if adaptive
+                        wsp, phasep, _ = downsample(ws, phasedata, 
+                                                  maximum(phasedata) - minimum(phasedata)/500)
+                        lines!(ax_phase, wsp, phasep)
+                    else
+                        lines!(ax_phase, ws, phasedata)
+                    end
+                    
+                    # Phase margin lines
+                    wpm_display = hz ? wpm ./ (2π) : wpm
+                    for k in 1:length(pm)
+                        phase_line = fullPhase[k] - pm[k]
+                        hlines!(ax_phase, phase_line, color=:gray, linestyle=:dash, alpha=0.5)
+                        lines!(ax_phase, [wpm_display[k], wpm_display[k]], 
+                              [fullPhase[k], phase_line], color=:blue)
+                        scatter!(ax_phase, [wpm_display[k]], [fullPhase[k]], 
+                                color=:blue, markersize=8)
+                    end
+                    
+                    # Add title with phase margins info
+                    if si == length(systems)
+                        pm_str = "Pm: [" * join([Printf.@sprintf("%.2g°", p) for p in pm], ", ") * "]"
+                        wpm_str = "ωpm: [" * join([Printf.@sprintf("%.2g", w) for w in wpm_display], ", ") * "]"
+                        ax_phase.title = pm_str * " " * wpm_str
+                    end
+                end
+            end
+        end
+    end
+    
+    return fig
+end
+
+# ====== Root Locus Plot ======
+function CSMakie.rlocusplot(P::LTISystem, K=500; output=false, kwargs...)
+    # Compute root locus
+    result = rlocus(P, K; output=output)
+    roots, Z, K_vals = result.roots, result.Z, result.K
+    
+    array_K = eltype(K_vals) <: AbstractArray
+    
+    fig = Figure()
+    ax = Axis(fig[1,1],
+             aspect = DataAspect(),
+             title = "Root Locus",
+             xlabel = "Re(roots)",
+             ylabel = "Im(roots)")
+    
+    # Add grid at zero
+    vlines!(ax, 0, color=:gray, alpha=0.3, linewidth=0.5)
+    hlines!(ax, 0, color=:gray, alpha=0.3, linewidth=0.5)
+    
+    redata = real.(roots)
+    imdata = imag.(roots)
+    
+    # Plot root locus branches
+    for i in 1:size(redata, 2)
+        lines!(ax, redata[:, i], imdata[:, i], color=:blue)
+    end
+    
+    # Plot zeros
+    if !isempty(Z)
+        scatter!(ax, real.(Z), imag.(Z), 
+                marker=:circle, markersize=10, 
+                color=(:transparent, 0.5), strokewidth=2, 
+                strokecolor=:green, label="Zeros")
+    end
+    
+    # Plot open-loop poles
+    scatter!(ax, redata[1, :], imdata[1, :], 
+            marker=:xcross, markersize=10, 
+            color=:red, label="Open-loop poles")
+    
+    # Plot closed-loop poles for matrix K
+    if array_K
+        scatter!(ax, redata[end, :], imdata[end, :], 
+                marker=:diamond, markersize=10, 
+                color=:purple, label="Closed-loop poles")
+    end
+    
+    return fig
+end
+
+# ====== RGA Plot ======
+function CSMakie.rgaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+                         w=nothing; hz=false, balance=true, kwargs...)
+    systems_vec = systems isa AbstractVector ? systems : [systems]
+    systems, w = isnothing(w) ? _processfreqplot(Val{:sigma}(), systems_vec) : 
+                                _processfreqplot(Val{:sigma}(), systems_vec, w)
+    
+    ws = (hz ? 1/(2π) : 1) .* w
+    
+    fig = Figure()
+    ax = Axis(fig[1,1],
+             xscale = log10,
+             title = "RGA Plot",
+             xlabel = hz ? "Frequency [Hz]" : "Frequency [rad/s]",
+             ylabel = "Element magnitudes")
+    
+    for (si, s) in enumerate(systems)
+        sbal = balance ? balance_statespace(s)[1] : s
+        rga = abs.(relative_gain_array(sbal, w))
+        
+        for j in 1:size(rga, 1)
+            for i in 1:size(rga, 2)
+                label = "System $si, from $i to $j"
+                lines!(ax, ws, rga[j, i, :], label=label)
+            end
+        end
+    end
+    
+    Legend(fig[1, 2], ax)
+    
+    return fig
+end
+
+# ====== LeadLinkCurve ======
+function CSMakie.leadlinkcurve(start=1; kwargs...)
+    N = range(start, stop=10, length=50)
+    dph = map(Ni -> 180/π*atan(sqrt(Ni)) - atan(1/sqrt(Ni)), N)
+    
+    fig = Figure()
+    ax = Axis(fig[1,1],
+             xlabel = "N",
+             ylabel = "Phase advance [deg]",
+             title = "Lead Link Phase Advance")
+    
+    lines!(ax, N, dph, color=:blue, linewidth=2)
+    
+    return fig
+end
+
+# ====== Nicholsplot ======
+# Note: This is a simplified version. The full implementation would require
+# porting all the complex grid calculations from the Plots version
+function CSMakie.nicholsplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+                             w=nothing; text=true, 
+                             Gains=[12, 6, 3, 1, 0.5, -0.5, -1, -3, -6, -10, -20, -40, -60],
+                             pInc=30, kwargs...)
+    systems_vec = systems isa AbstractVector ? systems : [systems]
+    systems, w = isnothing(w) ? _processfreqplot(Val{:nyquist}(), systems_vec) : 
+                                _processfreqplot(Val{:nyquist}(), systems_vec, w)
+    
+    fig = Figure()
+    ax = Axis(fig[1,1],
+             title = "Nichols Chart",
+             xlabel = "Phase [deg]",
+             ylabel = "Magnitude [dB]")
+    
+    # Plot systems
+    for (sysi, s) in enumerate(systems)
+        ℜresp, ℑresp = nyquist(s, w)[1:2]
+        ℜdata = dropdims(ℜresp, dims=(1,2))
+        ℑdata = dropdims(ℑresp, dims=(1,2))
+        mag = 20*log10.(sqrt.(ℜdata.^2 + ℑdata.^2))
+        angles = 180/π*angle.(ℜdata .+ im*ℑdata)
+        
+        lines!(ax, angles, mag, linewidth=2, label="System $sysi")
+    end
+    
+    # Note: Full Nichols grid implementation would go here
+    # This is a placeholder for the complex grid calculations
+    
+    return fig
+end
+
+# ====== Recipes for Types (SimResult, StepInfo) ======
+# These use Makie's standard recipe system for plotting types
+
+@recipe(SimResultPlot, r) do scene
+    Attributes(
+        plotu = false,
+        plotx = false,
+        ploty = true,
+    )
+end
+
+function Makie.plot!(srp::SimResultPlot)
+    r = srp.r[]
+    
+    # Handle vector of SimResults
+    if r isa AbstractVector
+        for res in r
+            plot!(current_axis(), res; 
+                  plotu=srp.plotu[], plotx=srp.plotx[], ploty=srp.ploty[])
+        end
+        return srp
+    end
+    
+    plotu = srp.plotu[]
+    plotx = srp.plotx[]
+    ploty = srp.ploty[]
+    
+    ny, nu, nx = r.ny, r.nu, r.nx
+    t = r.t
+    n_series = size(r.y, 3)
+    
+    # Calculate number of subplots needed
+    nplots = ploty ? ny : 0
+    plotu && (nplots += nu)
+    plotx && (nplots += nx)
+    
+    # Get current figure or create new one
+    fig = current_figure()
+    
+    # Create grid layout for subplots
+    gl = GridLayout(fig[1, 1])
+    
+    plotind = 1
+    axes = []
+    
+    # Plot outputs
+    if ploty
+        for i in 1:ny
+            ax = Axis(gl[plotind, 1], 
+                     xlabel = plotind == nplots ? "Time (s)" : "",
+                     ylabel = output_names(r.sys)[i])
+            push!(axes, ax)
+            
+            for ms in 1:n_series
+                label = n_series > 1 ? "From $(input_names(r.sys)[ms])" : ""
+                if iscontinuous(r.sys)
+                    lines!(ax, t, r.y[i, :, ms], label=label)
+                else
+                    stairs!(ax, t, r.y[i, :, ms], label=label, step=:post)
+                end
+            end
+            plotind += 1
+        end
+    end
+    
+    # Plot inputs
+    if plotu
+        for i in 1:nu
+            ax = Axis(gl[plotind, 1],
+                     xlabel = plotind == nplots ? "Time (s)" : "",
+                     ylabel = input_names(r.sys)[i])
+            push!(axes, ax)
+            
+            if iscontinuous(r.sys)
+                lines!(ax, t, r.u[i, :])
+            else
+                stairs!(ax, t, r.u[i, :], step=:post)
+            end
+            plotind += 1
+        end
+    end
+    
+    # Plot states
+    if plotx
+        for i in 1:nx
+            ax = Axis(gl[plotind, 1],
+                     xlabel = plotind == nplots ? "Time (s)" : "",
+                     ylabel = state_names(r.sys)[i])
+            push!(axes, ax)
+            
+            if n_series > 1
+                for ms in 1:n_series
+                    if iscontinuous(r.sys)
+                        lines!(ax, t, r.x[i, :, ms])
+                    else
+                        stairs!(ax, t, r.x[i, :, ms], step=:post)
+                    end
+                end
+            else
+                if iscontinuous(r.sys)
+                    lines!(ax, t, r.x[i, :])
+                else
+                    stairs!(ax, t, r.x[i, :], step=:post)
+                end
+            end
+            plotind += 1
+        end
+    end
+    
+    # Link x-axes
+    if length(axes) > 1
+        linkxaxes!(axes...)
+    end
+    
+    srp
+end
+
+@recipe(StepInfoPlot, si) do scene
+    Attributes()
+end
+
+function Makie.plot!(sip::StepInfoPlot)
+    si = sip.si[]
+    
+    ax = current_axis()
+    
+    # Plot the step response
+    t = si.res.t
+    y = si.res.y[1, :]
+    
+    lines!(ax, t, y, color=:blue, linewidth=2)
+    
+    # Final value
+    hlines!(ax, si.yf, color=:black, linestyle=:dash, 
+           label=Printf.@sprintf("Final value: %.3f", si.yf))
+    
+    # Rise time
+    lines!(ax, t[si.i10:si.i90], y[si.i10:si.i90], 
+          color=:green, linewidth=3,
+          label=Printf.@sprintf("Rise time: %.3f", si.risetime))
+    
+    vlines!(ax, [t[si.i10], t[si.i90]], color=:green, linestyle=:dash, alpha=0.5)
+    
+    # Peak and overshoot
+    scatter!(ax, [si.peaktime], [si.peak], color=:red, markersize=10,
+            label=Printf.@sprintf("Peak: %.3f, Overshoot: %.1f%%", si.peak, si.overshoot))
+    lines!(ax, [si.peaktime, si.peaktime], [si.y0, si.peak], color=:red, alpha=0.5)
+    
+    # Settling time
+    scatter!(ax, [si.settlingtime], [y[si.settlingtimeind]], color=:purple, markersize=10,
+            label=Printf.@sprintf("Settling time: %.3f", si.settlingtime))
+    lines!(ax, [si.settlingtime, si.settlingtime], [si.y0, y[si.settlingtimeind]], 
+          color=:purple, alpha=0.5)
+    
+    # Settling threshold bands
+    hlines!(ax, [si.yf - si.stepsize*si.settling_th, si.yf + si.stepsize*si.settling_th],
+           color=:purple, linestyle=:dash, alpha=0.3,
+           label=Printf.@sprintf("Settling threshold: %.1f%%", 100*si.settling_th))
+    
+    # Undershoot if present
+    if si.undershoot != 0
+        t_under = t[si.lowerpeakind]
+        scatter!(ax, [t_under], [si.lowerpeak], color=:orange, markersize=10,
+                label=Printf.@sprintf("Undershoot: %.1f%%", si.undershoot))
+        lines!(ax, [t_under, t_under], [si.y0, si.lowerpeak], color=:orange, alpha=0.5)
+    end
+    
+    ax.xlabel = "Time (s)"
+    ax.ylabel = "Output"
+    ax.title = "Step Response Analysis"
+    
+    axislegend(ax, position=:rt)
+    
+    sip
+end
+
+end # module
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/src/CSMakie.jl b/lib/ControlSystemsBase/src/CSMakie.jl
new file mode 100644
index 000000000..c212c649a
--- /dev/null
+++ b/lib/ControlSystemsBase/src/CSMakie.jl
@@ -0,0 +1,164 @@
+"""
+    CSMakie
+
+Module providing Makie.jl plotting functions for ControlSystemsBase.
+These functions are loaded when Makie.jl is available through the package extension system.
+
+## Usage
+```julia
+using ControlSystemsBase, GLMakie  # or CairoMakie
+CSMakie.bodeplot(sys)
+CSMakie.nyquistplot(sys)
+# etc.
+```
+
+All functions will throw an informative error if called without Makie.jl loaded.
+"""
+module CSMakie
+
+using ..ControlSystemsBase: LTISystem
+
+# Error message for when Makie is not loaded
+const MAKIE_NOT_LOADED_ERROR = """
+Makie.jl must be loaded to use CSMakie plotting functions.
+Please run one of:
+  using GLMakie      # for interactive plots
+  using CairoMakie   # for static plots
+  using WGLMakie     # for web-based plots
+"""
+
+# Define stub functions that will be overloaded by the extension
+
+"""
+    bodeplot(sys::LTISystem, w=nothing; kwargs...)
+    bodeplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
+
+Create a Bode plot using Makie.jl. Requires Makie to be loaded.
+
+See the main `bodeplot` documentation for available options.
+"""
+function bodeplot(args...; kwargs...)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+"""
+    nyquistplot(sys::LTISystem, w=nothing; kwargs...)
+    nyquistplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
+
+Create a Nyquist plot using Makie.jl. Requires Makie to be loaded.
+
+See the main `nyquistplot` documentation for available options.
+"""
+function nyquistplot(args...; kwargs...)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+"""
+    sigmaplot(sys::LTISystem, w=nothing; kwargs...)
+    sigmaplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
+
+Create a singular value plot using Makie.jl. Requires Makie to be loaded.
+
+See the main `sigmaplot` documentation for available options.
+"""
+function sigmaplot(args...; kwargs...)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+"""
+    marginplot(sys::LTISystem, w=nothing; kwargs...)
+    marginplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
+
+Create a margin plot using Makie.jl. Requires Makie to be loaded.
+
+See the main `marginplot` documentation for available options.
+"""
+function marginplot(args...; kwargs...)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+"""
+    pzmap(sys::LTISystem; kwargs...)
+    pzmap(sys::Vector{<:LTISystem}; kwargs...)
+
+Create a pole-zero map using Makie.jl. Requires Makie to be loaded.
+
+See the main `pzmap` documentation for available options.
+"""
+function pzmap(args...; kwargs...)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+"""
+    pzmap!(ax, sys::LTISystem; kwargs...)
+    pzmap!(ax, sys::Vector{<:LTISystem}; kwargs...)
+
+Add a pole-zero map to an existing Makie axis. Requires Makie to be loaded.
+"""
+function pzmap!(args...; kwargs...)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+"""
+    nicholsplot(sys::LTISystem, w=nothing; kwargs...)
+    nicholsplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
+
+Create a Nichols chart using Makie.jl. Requires Makie to be loaded.
+
+See the main `nicholsplot` documentation for available options.
+"""
+function nicholsplot(args...; kwargs...)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+"""
+    rgaplot(sys::LTISystem, w=nothing; kwargs...)
+    rgaplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
+
+Create a relative gain array plot using Makie.jl. Requires Makie to be loaded.
+
+See the main `rgaplot` documentation for available options.
+"""
+function rgaplot(args...; kwargs...)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+"""
+    rlocusplot(P::LTISystem, K=500; kwargs...)
+
+Create a root locus plot using Makie.jl. Requires Makie to be loaded.
+
+See the main `rlocusplot` documentation for available options.
+"""
+function rlocusplot(args...; kwargs...)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+"""
+    leadlinkcurve(start=1; kwargs...)
+
+Plot the phase advance curve for a lead link using Makie.jl. Requires Makie to be loaded.
+
+See the main `leadlinkcurve` documentation for available options.
+"""
+function leadlinkcurve(args...; kwargs...)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+"""
+    setPlotScale(str)
+
+Set the default scale of magnitude in `bodeplot` and `sigmaplot` for Makie plots.
+`str` should be either `"dB"` or `"log10"`. The default scale is `"log10"`.
+
+Requires Makie to be loaded.
+"""
+function setPlotScale(str::AbstractString)
+    error(MAKIE_NOT_LOADED_ERROR)
+end
+
+# Also export the mutating versions that will be defined by the extension
+export bodeplot, nyquistplot, sigmaplot, marginplot, pzmap, pzmap!,
+       nicholsplot, rgaplot, rlocusplot, leadlinkcurve, setPlotScale
+
+end # module CSMakie
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/src/ControlSystemsBase.jl b/lib/ControlSystemsBase/src/ControlSystemsBase.jl
index 505e2f872..2524e4ac3 100644
--- a/lib/ControlSystemsBase/src/ControlSystemsBase.jl
+++ b/lib/ControlSystemsBase/src/ControlSystemsBase.jl
@@ -216,6 +216,10 @@ include("plotting.jl")
 
 include("root_locus.jl")
 
+# CSMakie module for Makie plotting support
+include("CSMakie.jl")
+export CSMakie
+
 @deprecate pole poles
 @deprecate tzero tzeros
 @deprecate num numvec
diff --git a/lib/ControlSystemsBase/test/test_makie_plots.jl b/lib/ControlSystemsBase/test/test_makie_plots.jl
new file mode 100644
index 000000000..5e320ae6a
--- /dev/null
+++ b/lib/ControlSystemsBase/test/test_makie_plots.jl
@@ -0,0 +1,160 @@
+using ControlSystemsBase
+using Test
+using LinearAlgebra
+
+# Check if Makie is available
+using Makie
+using CairoMakie # Use CairoMakie for non-interactive backend in tests
+
+@testset "Makie Plot Tests" begin
+    # Create test systems
+    P = tf([1], [1, 2, 1])
+    P2 = tf([1, 2], [1, 3, 2])
+    Pss = ss(P)
+    Pmimo = [P P2; P2 P]
+    Pdiscrete = c2d(P, 0.1)
+    
+    # Test frequency vector
+    w = exp10.(range(-2, stop=2, length=100))
+    
+    @testset "pzmap" begin
+        @test_nowarn begin
+            fig = CSMakie.pzmap(P)
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.pzmap([P, P2])
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.pzmap(Pdiscrete; hz=true)
+            @test fig isa Makie.Figure
+        end
+    end
+    
+    @testset "bodeplot" begin
+        @test_nowarn begin
+            fig = CSMakie.bodeplot(P)
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.bodeplot(P, w)
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.bodeplot([P, P2]; plotphase=false)
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.bodeplot(Pmimo; hz=true)
+            @test fig isa Makie.Figure
+        end
+        # Test dB scale
+        @test_nowarn begin
+            ControlSystemsBase.setPlotScale("dB")
+            fig = CSMakie.bodeplot(P)
+            @test fig isa Makie.Figure
+            ControlSystemsBase.setPlotScale("log10")
+        end
+    end
+    
+    @testset "nyquistplot" begin
+        @test_nowarn begin
+            fig = CSMakie.nyquistplot(P)
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.nyquistplot(P, w; Ms_circles=[1.5, 2.0])
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.nyquistplot([P, P2]; unit_circle=true)
+            @test fig isa Makie.Figure
+        end
+    end
+    
+    @testset "sigmaplot" begin
+        @test_nowarn begin
+            fig = CSMakie.sigmaplot(P)
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.sigmaplot(Pmimo; extrema=true)
+            @test fig isa Makie.Figure
+        end
+    end
+    
+    @testset "marginplot" begin
+        @test_nowarn begin
+            fig = CSMakie.marginplot(P)
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.marginplot(P; plotphase=false)
+            @test fig isa Makie.Figure
+        end
+    end
+    
+    @testset "rlocusplot" begin
+        @test_nowarn begin
+            fig = CSMakie.rlocusplot(P)
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.rlocusplot(P, 100)
+            @test fig isa Makie.Figure
+        end
+    end
+    
+    @testset "rgaplot" begin
+        @test_nowarn begin
+            fig = CSMakie.rgaplot(Pmimo)
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.rgaplot(Pmimo, w; hz=true)
+            @test fig isa Makie.Figure
+        end
+    end
+    
+    @testset "SimResult plot" begin
+        # Create a simple simulation result
+        t = 0:0.01:5
+        res = step(P, t)
+        @test_nowarn begin
+            fig = plot(res)
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = plot(res; plotu=true, plotx=true)
+            @test fig isa Makie.Figure
+        end
+    end
+    
+    @testset "StepInfo plot" begin
+        res = step(P, 10)
+        si = stepinfo(res)
+        @test_nowarn begin
+            fig = plot(si)
+            @test fig isa Makie.Figure
+        end
+    end
+    
+    @testset "leadlinkcurve" begin
+        @test_nowarn begin
+            fig = CSMakie.leadlinkcurve()
+            @test fig isa Makie.Figure
+        end
+        @test_nowarn begin
+            fig = CSMakie.leadlinkcurve(2)
+            @test fig isa Makie.Figure
+        end
+    end
+    
+    @testset "nicholsplot" begin
+        @test_nowarn begin
+            fig = nicholsCSMakie.plot(P)
+            @test fig isa Makie.Figure
+        end
+    end
+end
\ No newline at end of file

From 80019f341973fed46b8f34a5a056fc74212d382d Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 25 Aug 2025 08:00:02 +0200
Subject: [PATCH 02/10] fixes

---
 .../ext/ControlSystemsBaseMakieExt.jl         | 75 +++++++------------
 lib/ControlSystemsBase/src/CSMakie.jl         | 14 +---
 .../test/test_makie_plots.jl                  |  7 +-
 3 files changed, 29 insertions(+), 67 deletions(-)

diff --git a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
index 8742af43e..b0c1ae065 100644
--- a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
+++ b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
@@ -14,36 +14,12 @@ using ControlSystemsBase: downsample, _processfreqplot, _default_freq_vector,
                           sisomargin, relative_gain_array, rlocus,
                           input_names, output_names, state_names, system_name,
                           iscontinuous, isdiscrete, issiso, isrational,
-                          integrator_excess, balance_statespace, LTISystem
+                          integrator_excess, balance_statespace, LTISystem,
+                          _PlotScale, _PlotScaleFunc, _PlotScaleStr, _span  # Use existing plot scale settings
 
-# Global plot scale settings (matching Plots.jl version)
-const _PlotScale = Ref("log10")
-const _PlotScaleFunc = Ref(:log10)
-const _PlotScaleStr = Ref("")
-
-"""
-    CSMakie.setPlotScale(str)
-
-Set the default scale of magnitude in `bodeplot` and `sigmaplot` for Makie plots.
-`str` should be either `"dB"` or `"log10"`. The default scale if none is chosen is `"log10"`.
-"""
-function CSMakie.setPlotScale(str::AbstractString)
-    if str == "dB"
-        _PlotScale[] = str
-        _PlotScaleFunc[] = :identity
-        _PlotScaleStr[] = "(dB)"
-    elseif str == "log10"
-        _PlotScale[] = str
-        _PlotScaleFunc[] = :log10
-        _PlotScaleStr[] = ""
-    else
-        error("Scale must be set to either \"dB\" or \"log10\"")
-    end
-end
-
-# Helper function to get y-scale transform
+# Helper function to get y-scale transform for Makie
 function get_yscale_transform()
-    _PlotScaleFunc[] == :log10 ? Makie.log10 : identity
+    ControlSystemsBase._PlotScaleFunc == :log10 ? Makie.log10 : identity
 end
 
 # ====== PZMap (Pole-Zero Map) ======
@@ -163,7 +139,7 @@ function CSMakie.bodeplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}
                          xscale = log10,
                          yscale = get_yscale_transform(),
                          xlabel = plotphase ? "" : (hz ? "Frequency [Hz]" : "Frequency [rad/s]"),
-                         ylabel = i == 1 ? "Magnitude $_PlotScaleStr[]" : "",
+                         ylabel = i == 1 ? "Magnitude $(ControlSystemsBase._PlotScaleStr)" : "",
                          title = j == 1 && i == 1 ? "Bode Plot" : "")
             axes_mag[i, j] = ax_mag
             
@@ -190,7 +166,7 @@ function CSMakie.bodeplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}
         
         mag, phase = bode(sbal, w; unwrap=false)
         
-        if _PlotScale[] == "dB"
+        if ControlSystemsBase._PlotScale == "dB"
             mag = 20*log10.(mag)
         elseif 0 ∈ mag
             replace!(mag, 0 => NaN)
@@ -208,8 +184,7 @@ function CSMakie.bodeplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}
                 lab = _get_plotlabel(s, i, j)
                 
                 if adaptive && eltype(magdata) <: AbstractFloat
-                    lmag = _PlotScale[] == "dB" ? magdata : log.(magdata)
-                    wsi, magi, _ = downsample(ws, lmag, maximum(lmag) - minimum(lmag)/500)
+                    wsi, magi, _ = downsample(ws, magdata, _span(magdata)/500)
                     lines!(ax_mag, wsi, magi, label=lab)
                 else
                     lines!(ax_mag, ws, magdata, label=lab)
@@ -231,7 +206,7 @@ function CSMakie.bodeplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}
                     ax_phase = axes_phase[i, j]
                     if adaptive && eltype(phasedata) <: AbstractFloat
                         wsp, phasep, _ = downsample(ws, phasedata, 
-                                                   maximum(phasedata) - minimum(phasedata)/500)
+                                                   _span(phasedata)/500)
                         lines!(ax_phase, wsp, phasep)
                     else
                         lines!(ax_phase, ws, phasedata)
@@ -291,8 +266,8 @@ function CSMakie.nyquistplot(systems::Union{LTISystem, AbstractVector{<:LTISyste
                 lab = _get_plotlabel(s, i, j)
                 
                 if adaptive
-                    indsre, _, _ = downsample(w, redata, 1/500)
-                    indsim, _, _ = downsample(w, imdata, 1/500)
+                    indsre = downsample(w, redata, 1/500)[3]
+                    indsim = downsample(w, imdata, 1/500)[3]
                     inds = sort!(union(indsre, indsim))
                     lines!(ax, redata[inds], imdata[inds], label=lab)
                 else
@@ -349,7 +324,7 @@ function CSMakie.sigmaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}
              yscale = get_yscale_transform(),
              title = "Sigma Plot",
              xlabel = hz ? "Frequency [Hz]" : "Frequency [rad/s]",
-             ylabel = "Singular Values $_PlotScaleStr[]")
+             ylabel = "Singular Values $(ControlSystemsBase._PlotScaleStr)")
     
     for (si, s) in enumerate(systems)
         sbal = balance ? balance_statespace(s)[1] : s
@@ -359,10 +334,9 @@ function CSMakie.sigmaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}
             sv = sv[:, [1, end]]
         end
         
-        if _PlotScale[] == "dB"
+        if ControlSystemsBase._PlotScale == "dB"
             sv = 20*log10.(sv)
         end
-        
         # Use _to1series to efficiently plot multiple lines
         ws_series, sv_series = _to1series(ws, sv)
         lines!(ax, ws_series, sv_series, color=Cycled(si))
@@ -442,7 +416,7 @@ function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem
                 ax_mag = axes_mag[i, j]
                 magdata = vec(bmag[i, j, :])
                 
-                if _PlotScale[] == "dB"
+                if ControlSystemsBase._PlotScale == "dB"
                     magdata = 20*log10.(magdata)
                     mag_margins = 20 .* log10.(1 ./ gm)
                     oneLine = 0
@@ -454,7 +428,7 @@ function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem
                 # Plot magnitude response
                 if adaptive
                     wsi, magi, _ = downsample(ws, magdata, 
-                                            maximum(magdata) - minimum(magdata)/500)
+                                            _span(magdata)/500)
                     lines!(ax_mag, wsi, magi, label=_get_plotlabel(s, i, j))
                 else
                     lines!(ax_mag, ws, magdata, label=_get_plotlabel(s, i, j))
@@ -463,13 +437,12 @@ function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem
                 # Unity gain line
                 hlines!(ax_mag, oneLine, color=:gray, linestyle=:dash, alpha=0.5)
                 
-                # Gain margins
+                # Gain margins - draw vertical lines from unity to margin value
                 wgm_display = hz ? wgm ./ (2π) : wgm
-                for k in 1:length(gm)
+                for k in eachindex(wgm_display)
+                    # Draw vertical line from unity gain to the margin value
                     lines!(ax_mag, [wgm_display[k], wgm_display[k]], 
-                          [oneLine, mag_margins[k]], color=:red)
-                    scatter!(ax_mag, [wgm_display[k]], [mag_margins[k]], 
-                            color=:red, markersize=8)
+                          [oneLine, mag_margins[k]], color=:red, linewidth=2)
                 end
                 
                 # Add title with margins info
@@ -485,7 +458,7 @@ function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem
                     
                     if adaptive
                         wsp, phasep, _ = downsample(ws, phasedata, 
-                                                  maximum(phasedata) - minimum(phasedata)/500)
+                                                  _span(phasedata)/500)
                         lines!(ax_phase, wsp, phasep)
                     else
                         lines!(ax_phase, ws, phasedata)
@@ -493,13 +466,17 @@ function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem
                     
                     # Phase margin lines
                     wpm_display = hz ? wpm ./ (2π) : wpm
+                    
+                    # Draw horizontal lines at phase margin crossings
                     for k in 1:length(pm)
                         phase_line = fullPhase[k] - pm[k]
                         hlines!(ax_phase, phase_line, color=:gray, linestyle=:dash, alpha=0.5)
+                    end
+                    
+                    # Draw vertical lines showing the phase margins
+                    for k in 1:length(pm)
                         lines!(ax_phase, [wpm_display[k], wpm_display[k]], 
-                              [fullPhase[k], phase_line], color=:blue)
-                        scatter!(ax_phase, [wpm_display[k]], [fullPhase[k]], 
-                                color=:blue, markersize=8)
+                              [fullPhase[k], fullPhase[k] - pm[k]], color=:blue, linewidth=2)
                     end
                     
                     # Add title with phase margins info
diff --git a/lib/ControlSystemsBase/src/CSMakie.jl b/lib/ControlSystemsBase/src/CSMakie.jl
index c212c649a..9403faef0 100644
--- a/lib/ControlSystemsBase/src/CSMakie.jl
+++ b/lib/ControlSystemsBase/src/CSMakie.jl
@@ -145,20 +145,8 @@ function leadlinkcurve(args...; kwargs...)
     error(MAKIE_NOT_LOADED_ERROR)
 end
 
-"""
-    setPlotScale(str)
-
-Set the default scale of magnitude in `bodeplot` and `sigmaplot` for Makie plots.
-`str` should be either `"dB"` or `"log10"`. The default scale is `"log10"`.
-
-Requires Makie to be loaded.
-"""
-function setPlotScale(str::AbstractString)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
-
 # Also export the mutating versions that will be defined by the extension
 export bodeplot, nyquistplot, sigmaplot, marginplot, pzmap, pzmap!,
-       nicholsplot, rgaplot, rlocusplot, leadlinkcurve, setPlotScale
+       nicholsplot, rgaplot, rlocusplot, leadlinkcurve
 
 end # module CSMakie
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/test/test_makie_plots.jl b/lib/ControlSystemsBase/test/test_makie_plots.jl
index 5e320ae6a..2ce9bffa1 100644
--- a/lib/ControlSystemsBase/test/test_makie_plots.jl
+++ b/lib/ControlSystemsBase/test/test_makie_plots.jl
@@ -1,10 +1,7 @@
 using ControlSystemsBase
 using Test
 using LinearAlgebra
-
-# Check if Makie is available
-using Makie
-using CairoMakie # Use CairoMakie for non-interactive backend in tests
+using GLMakie
 
 @testset "Makie Plot Tests" begin
     # Create test systems
@@ -153,7 +150,7 @@ using CairoMakie # Use CairoMakie for non-interactive backend in tests
     
     @testset "nicholsplot" begin
         @test_nowarn begin
-            fig = nicholsCSMakie.plot(P)
+            fig = CSMakie.nicholsplot(P)
             @test fig isa Makie.Figure
         end
     end

From b860cc8289f2477e029f12a7ce1a25b3a120d57e Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 25 Aug 2025 08:23:58 +0200
Subject: [PATCH 03/10] fix simresultplot

---
 .../ext/ControlSystemsBaseMakieExt.jl         | 59 +++++--------------
 1 file changed, 15 insertions(+), 44 deletions(-)

diff --git a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
index b0c1ae065..32ab76ca3 100644
--- a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
+++ b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
@@ -626,33 +626,11 @@ function CSMakie.nicholsplot(systems::Union{LTISystem, AbstractVector{<:LTISyste
     return fig
 end
 
-# ====== Recipes for Types (SimResult, StepInfo) ======
-# These use Makie's standard recipe system for plotting types
+# ====== Direct plot methods for types ======
+# These allow direct plotting of types with plot()
 
-@recipe(SimResultPlot, r) do scene
-    Attributes(
-        plotu = false,
-        plotx = false,
-        ploty = true,
-    )
-end
-
-function Makie.plot!(srp::SimResultPlot)
-    r = srp.r[]
-    
-    # Handle vector of SimResults
-    if r isa AbstractVector
-        for res in r
-            plot!(current_axis(), res; 
-                  plotu=srp.plotu[], plotx=srp.plotx[], ploty=srp.ploty[])
-        end
-        return srp
-    end
-    
-    plotu = srp.plotu[]
-    plotx = srp.plotx[]
-    ploty = srp.ploty[]
-    
+# Direct plot method for SimResult
+function Makie.plot(r::SimResult; plotu=false, plotx=false, ploty=true)
     ny, nu, nx = r.ny, r.nu, r.nx
     t = r.t
     n_series = size(r.y, 3)
@@ -662,10 +640,8 @@ function Makie.plot!(srp::SimResultPlot)
     plotu && (nplots += nu)
     plotx && (nplots += nx)
     
-    # Get current figure or create new one
-    fig = current_figure()
-    
-    # Create grid layout for subplots
+    # Create figure with grid layout
+    fig = Figure()
     gl = GridLayout(fig[1, 1])
     
     plotind = 1
@@ -740,17 +716,16 @@ function Makie.plot!(srp::SimResultPlot)
         linkxaxes!(axes...)
     end
     
-    srp
-end
-
-@recipe(StepInfoPlot, si) do scene
-    Attributes()
+    return fig
 end
 
-function Makie.plot!(sip::StepInfoPlot)
-    si = sip.si[]
-    
-    ax = current_axis()
+# Direct plot method for StepInfo
+function Makie.plot(si::StepInfo)
+    fig = Figure()
+    ax = Axis(fig[1,1],
+             xlabel = "Time (s)",
+             ylabel = "Output",
+             title = "Step Response Analysis")
     
     # Plot the step response
     t = si.res.t
@@ -793,13 +768,9 @@ function Makie.plot!(sip::StepInfoPlot)
         lines!(ax, [t_under, t_under], [si.y0, si.lowerpeak], color=:orange, alpha=0.5)
     end
     
-    ax.xlabel = "Time (s)"
-    ax.ylabel = "Output"
-    ax.title = "Step Response Analysis"
-    
     axislegend(ax, position=:rt)
     
-    sip
+    return fig
 end
 
 end # module
\ No newline at end of file

From 037e4064e137671ce561f5be92dd073f10f96d5a Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 25 Aug 2025 08:32:29 +0200
Subject: [PATCH 04/10] add docs

---
 docs/src/lib/plotting.md                | 45 +++++++++++++++++++++++++
 lib/ControlSystemsBase/test/runtests.jl |  1 +
 2 files changed, 46 insertions(+)

diff --git a/docs/src/lib/plotting.md b/docs/src/lib/plotting.md
index 74c1ed311..faab520c3 100644
--- a/docs/src/lib/plotting.md
+++ b/docs/src/lib/plotting.md
@@ -8,6 +8,51 @@ Pages = ["plotting.md"]
 !!! note "Time-domain responses"
     There are no special functions to plot time-domain results, such as step and impulse responses, instead, simply call `plot` on the result structure (`ControlSystemsBase.SimResult`) returned by [`lsim`](@ref), [`step`](@ref), [`impulse`](@ref) etc.
 
+## Makie support
+!!! danger "Experimental"
+    
+    The support for plotting with Makie is currently experimental and at any time subject to breaking changes or removal **not** respecting semantic versioning.
+
+ControlSystemsBase provides experimental support for plotting with [Makie.jl](https://docs.makie.org/) through the `CSMakie` module. This support is loaded automatically when you load a Makie backend (GLMakie, CairoMakie, or WGLMakie).
+
+### Usage
+
+```julia
+using ControlSystemsBase, GLMakie  # or CairoMakie, WGLMakie
+
+# Create a system
+P = tf([1], [1, 2, 1])
+
+# Use CSMakie plotting functions
+CSMakie.bodeplot(P)
+CSMakie.nyquistplot(P)
+CSMakie.pzmap(P)
+# ... and more
+
+# Direct plotting of simulation results
+res = step(P, 10)
+plot(res)  # Creates a figure with time-domain response
+
+si = stepinfo(res)
+plot(si)  # Visualizes step response characteristics
+```
+
+### Available functions
+
+The `CSMakie` module provides Makie implementations of the following plotting functions:
+
+- `CSMakie.bodeplot` - Bode magnitude and phase plots
+- `CSMakie.nyquistplot` - Nyquist plots with optional M and Mt circles
+- `CSMakie.sigmaplot` - Singular value plots
+- `CSMakie.marginplot` - Gain and phase margin plots
+- `CSMakie.pzmap` - Pole-zero maps
+- `CSMakie.nicholsplot` - Nichols charts
+- `CSMakie.rgaplot` - Relative gain array plots
+- `CSMakie.rlocusplot` - Root locus plots
+- `CSMakie.leadlinkcurve` - Lead-link design curves
+
+Additionally, `SimResult` and `StepInfo` types can be plotted directly using Makie's `plot` function.
+
 # Plotting functions
 
 ```@autodocs
diff --git a/lib/ControlSystemsBase/test/runtests.jl b/lib/ControlSystemsBase/test/runtests.jl
index a39dc8566..34342de32 100644
--- a/lib/ControlSystemsBase/test/runtests.jl
+++ b/lib/ControlSystemsBase/test/runtests.jl
@@ -46,6 +46,7 @@ my_tests = [
             "test_demo_systems",
             "test_autovec",
             "test_plots",
+            "test_makie_plots",
             "test_dsp",
             "test_implicit_diff",
             "test_rootlocus",

From b49025aa718ec365c402c3a5a730cabdbb968ded Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 25 Aug 2025 08:44:14 +0200
Subject: [PATCH 05/10] add angle unwrapping in nicholsplot

---
 .../ext/ControlSystemsBaseMakieExt.jl            | 16 +++-------------
 lib/ControlSystemsBase/test/test_makie_plots.jl  |  4 ++--
 2 files changed, 5 insertions(+), 15 deletions(-)

diff --git a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
index 32ab76ca3..af967da4b 100644
--- a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
+++ b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
@@ -111,7 +111,6 @@ function CSMakie.pzmap!(ax::Axis, systems::Union{LTISystem, AbstractVector{<:LTI
     return ax
 end
 
-# ====== Bodeplot ======
 function CSMakie.bodeplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                           w=nothing; plotphase=true, unwrap=true, hz=false, 
                           balance=true, adjust_phase_start=true, adaptive=true, kwargs...)
@@ -219,7 +218,6 @@ function CSMakie.bodeplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}
     return fig
 end
 
-# ====== Nyquistplot ======
 function CSMakie.nyquistplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                              w=nothing; Ms_circles=Float64[], Mt_circles=Float64[], 
                              unit_circle=false, hz=false, critical_point=-1, 
@@ -309,7 +307,6 @@ function CSMakie.nyquistplot(systems::Union{LTISystem, AbstractVector{<:LTISyste
     return fig
 end
 
-# ====== Sigmaplot ======
 function CSMakie.sigmaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                            w=nothing; hz=false, balance=true, extrema=false, kwargs...)
     systems_vec = systems isa AbstractVector ? systems : [systems]
@@ -345,7 +342,6 @@ function CSMakie.sigmaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}
     return fig
 end
 
-# ====== Marginplot ======
 function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                             w=nothing; plotphase=true, hz=false, balance=true, 
                             adjust_phase_start=true, adaptive=true, kwargs...)
@@ -493,7 +489,6 @@ function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem
     return fig
 end
 
-# ====== Root Locus Plot ======
 function CSMakie.rlocusplot(P::LTISystem, K=500; output=false, kwargs...)
     # Compute root locus
     result = rlocus(P, K; output=output)
@@ -543,7 +538,6 @@ function CSMakie.rlocusplot(P::LTISystem, K=500; output=false, kwargs...)
     return fig
 end
 
-# ====== RGA Plot ======
 function CSMakie.rgaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                          w=nothing; hz=false, balance=true, kwargs...)
     systems_vec = systems isa AbstractVector ? systems : [systems]
@@ -616,7 +610,8 @@ function CSMakie.nicholsplot(systems::Union{LTISystem, AbstractVector{<:LTISyste
         ℑdata = dropdims(ℑresp, dims=(1,2))
         mag = 20*log10.(sqrt.(ℜdata.^2 + ℑdata.^2))
         angles = 180/π*angle.(ℜdata .+ im*ℑdata)
-        
+        # unwrap angles
+        ControlSystemsBase.unwrap!(angles)
         lines!(ax, angles, mag, linewidth=2, label="System $sysi")
     end
     
@@ -626,10 +621,6 @@ function CSMakie.nicholsplot(systems::Union{LTISystem, AbstractVector{<:LTISyste
     return fig
 end
 
-# ====== Direct plot methods for types ======
-# These allow direct plotting of types with plot()
-
-# Direct plot method for SimResult
 function Makie.plot(r::SimResult; plotu=false, plotx=false, ploty=true)
     ny, nu, nx = r.ny, r.nu, r.nx
     t = r.t
@@ -684,7 +675,7 @@ function Makie.plot(r::SimResult; plotu=false, plotx=false, ploty=true)
         end
     end
     
-    # Plot states
+    # Plot state
     if plotx
         for i in 1:nx
             ax = Axis(gl[plotind, 1],
@@ -719,7 +710,6 @@ function Makie.plot(r::SimResult; plotu=false, plotx=false, ploty=true)
     return fig
 end
 
-# Direct plot method for StepInfo
 function Makie.plot(si::StepInfo)
     fig = Figure()
     ax = Axis(fig[1,1],
diff --git a/lib/ControlSystemsBase/test/test_makie_plots.jl b/lib/ControlSystemsBase/test/test_makie_plots.jl
index 2ce9bffa1..9651d66f0 100644
--- a/lib/ControlSystemsBase/test/test_makie_plots.jl
+++ b/lib/ControlSystemsBase/test/test_makie_plots.jl
@@ -5,7 +5,7 @@ using GLMakie
 
 @testset "Makie Plot Tests" begin
     # Create test systems
-    P = tf([1], [1, 2, 1])
+    P = tf([1.2], [1, 2, 1])*tf(1, [10, 1])
     P2 = tf([1, 2], [1, 3, 2])
     Pss = ss(P)
     Pmimo = [P P2; P2 P]
@@ -129,7 +129,7 @@ using GLMakie
     end
     
     @testset "StepInfo plot" begin
-        res = step(P, 10)
+        res = step(P, 100)
         si = stepinfo(res)
         @test_nowarn begin
             fig = plot(si)

From dd15a6b7a42b9e0e968128ffac2ed8c2339a4b4f Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 25 Aug 2025 11:26:54 +0200
Subject: [PATCH 06/10] add mutating versions

---
 .../ext/ControlSystemsBaseMakieExt.jl         | 116 +++++++-------
 lib/ControlSystemsBase/src/CSMakie.jl         | 142 +++---------------
 .../test/test_makie_plots.jl                  |  11 ++
 3 files changed, 86 insertions(+), 183 deletions(-)

diff --git a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
index af967da4b..e98076ae9 100644
--- a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
+++ b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
@@ -23,11 +23,14 @@ function get_yscale_transform()
 end
 
 # ====== PZMap (Pole-Zero Map) ======
-function CSMakie.pzmap(systems::Union{LTISystem, AbstractVector{<:LTISystem}}; 
+function CSMakie.pzmap(args...; kwargs...)
+    fig = Figure()
+    CSMakie.pzmap!(fig, args...; kwargs...)
+end
+function CSMakie.pzmap!(fig, systems::Union{LTISystem, AbstractVector{<:LTISystem}}; 
                        hz=false, kwargs...)
     systems_vec = systems isa AbstractVector ? systems : [systems]
     
-    fig = Figure()
     ax = Axis(fig[1,1], aspect = DataAspect(),
              title = "Pole-zero map",
              xlabel = hz ? "Real [Hz]" : "Real",
@@ -70,48 +73,13 @@ function CSMakie.pzmap(systems::Union{LTISystem, AbstractVector{<:LTISystem}};
     return fig
 end
 
-function CSMakie.pzmap!(ax::Axis, systems::Union{LTISystem, AbstractVector{<:LTISystem}}; 
-                        hz=false, kwargs...)
-    systems_vec = systems isa AbstractVector ? systems : [systems]
-    
-    # Add grid at zero
-    vlines!(ax, 0, color=:gray, alpha=0.3, linewidth=0.5)
-    hlines!(ax, 0, color=:gray, alpha=0.3, linewidth=0.5)
-    
-    scale_factor = hz ? 1/(2π) : 1
-    
-    for (i, system) in enumerate(systems_vec)
-        p = poles(system) .* scale_factor
-        z = tzeros(system) .* scale_factor
-        
-        # Plot zeros as circles
-        if !isempty(z)
-            scatter!(ax, real(z), imag(z), 
-                    marker=:circle, markersize=15, 
-                    color=(:transparent, 0.5), strokewidth=2, 
-                    strokecolor=Cycled(i))
-        end
-        
-        # Plot poles as x's
-        if !isempty(p)
-            scatter!(ax, real(p), imag(p), 
-                    marker=:xcross, markersize=15,
-                    color=Cycled(i))
-        end
-        
-        # Add unit circle for discrete systems
-        if isdiscrete(system) && i == 1  # Only draw once
-            θ = range(0, 2π, length=100)
-            circle_scale = scale_factor
-            lines!(ax, circle_scale .* cos.(θ), circle_scale .* sin.(θ), 
-                  color=:gray, linestyle=:dash, alpha=0.5)
-        end
-    end
-    
-    return ax
-end
 
-function CSMakie.bodeplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+
+function CSMakie.bodeplot(args...; kwargs...)
+    fig = Figure()
+    CSMakie.bodeplot!(fig, args...; kwargs...)
+end
+function CSMakie.bodeplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                           w=nothing; plotphase=true, unwrap=true, hz=false, 
                           balance=true, adjust_phase_start=true, adaptive=true, kwargs...)
     systems_vec = systems isa AbstractVector ? systems : [systems]
@@ -121,7 +89,6 @@ function CSMakie.bodeplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}
     ws = (hz ? 1/(2π) : 1) .* w
     ny, nu = size(systems[1])
     
-    fig = Figure()
     gl = GridLayout(fig[1, 1])
     
     # Create axes grid
@@ -218,7 +185,11 @@ function CSMakie.bodeplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}
     return fig
 end
 
-function CSMakie.nyquistplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+function CSMakie.nyquistplot(args...; kwargs...)
+    fig = Figure()
+    CSMakie.nyquistplot!(fig, args...; kwargs...)
+end
+function CSMakie.nyquistplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                              w=nothing; Ms_circles=Float64[], Mt_circles=Float64[], 
                              unit_circle=false, hz=false, critical_point=-1, 
                              balance=true, adaptive=true, kwargs...)
@@ -227,8 +198,7 @@ function CSMakie.nyquistplot(systems::Union{LTISystem, AbstractVector{<:LTISyste
                                 _processfreqplot(Val{:nyquist}(), systems_vec, w; adaptive)
     
     ny, nu = size(systems[1])
-    
-    fig = Figure()
+
     gl = GridLayout(fig[1, 1])
     
     # Create axes grid
@@ -307,7 +277,11 @@ function CSMakie.nyquistplot(systems::Union{LTISystem, AbstractVector{<:LTISyste
     return fig
 end
 
-function CSMakie.sigmaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+function CSMakie.sigmaplot(args...; kwargs...)
+    fig = Figure()
+    CSMakie.sigmaplot!(fig, args...; kwargs...)
+end
+function CSMakie.sigmaplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                            w=nothing; hz=false, balance=true, extrema=false, kwargs...)
     systems_vec = systems isa AbstractVector ? systems : [systems]
     systems, w = isnothing(w) ? _processfreqplot(Val{:sigma}(), systems_vec) : 
@@ -315,7 +289,6 @@ function CSMakie.sigmaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}
     
     ws = (hz ? 1/(2π) : 1) .* w
     
-    fig = Figure()
     ax = Axis(fig[1,1],
              xscale = log10,
              yscale = get_yscale_transform(),
@@ -342,7 +315,11 @@ function CSMakie.sigmaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}
     return fig
 end
 
-function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+function CSMakie.marginplot(args...; kwargs...)
+    fig = Figure()
+    CSMakie.marginplot!(fig, args...; kwargs...)
+end
+function CSMakie.marginplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                             w=nothing; plotphase=true, hz=false, balance=true, 
                             adjust_phase_start=true, adaptive=true, kwargs...)
     systems_vec = systems isa AbstractVector ? systems : [systems]
@@ -352,7 +329,6 @@ function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem
     ws = (hz ? 1/(2π) : 1) .* w
     ny, nu = size(systems[1])
     
-    fig = Figure()
     gl = GridLayout(fig[1, 1])
     
     # Create axes grid
@@ -489,14 +465,17 @@ function CSMakie.marginplot(systems::Union{LTISystem, AbstractVector{<:LTISystem
     return fig
 end
 
-function CSMakie.rlocusplot(P::LTISystem, K=500; output=false, kwargs...)
+function CSMakie.rlocusplot(args...; kwargs...)
+    fig = Figure()
+    CSMakie.rlocusplot!(fig, args...; kwargs...)
+end
+function CSMakie.rlocusplot!(fig, P::LTISystem, K=500; output=false, kwargs...)
     # Compute root locus
     result = rlocus(P, K; output=output)
     roots, Z, K_vals = result.roots, result.Z, result.K
     
     array_K = eltype(K_vals) <: AbstractArray
     
-    fig = Figure()
     ax = Axis(fig[1,1],
              aspect = DataAspect(),
              title = "Root Locus",
@@ -538,7 +517,11 @@ function CSMakie.rlocusplot(P::LTISystem, K=500; output=false, kwargs...)
     return fig
 end
 
-function CSMakie.rgaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+function CSMakie.rgaplot(args...; kwargs...)
+    fig = Figure()
+    CSMakie.rgaplot!(fig, args...; kwargs...)
+end
+function CSMakie.rgaplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                          w=nothing; hz=false, balance=true, kwargs...)
     systems_vec = systems isa AbstractVector ? systems : [systems]
     systems, w = isnothing(w) ? _processfreqplot(Val{:sigma}(), systems_vec) : 
@@ -546,7 +529,6 @@ function CSMakie.rgaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}},
     
     ws = (hz ? 1/(2π) : 1) .* w
     
-    fig = Figure()
     ax = Axis(fig[1,1],
              xscale = log10,
              title = "RGA Plot",
@@ -571,11 +553,14 @@ function CSMakie.rgaplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}},
 end
 
 # ====== LeadLinkCurve ======
-function CSMakie.leadlinkcurve(start=1; kwargs...)
+function CSMakie.leadlinkcurve(args...; kwargs...)
+    fig = Figure()
+    CSMakie.leadlinkcurve!(fig, args...; kwargs...)
+end
+function CSMakie.leadlinkcurve!(fig, start=1; kwargs...)
     N = range(start, stop=10, length=50)
     dph = map(Ni -> 180/π*atan(sqrt(Ni)) - atan(1/sqrt(Ni)), N)
     
-    fig = Figure()
     ax = Axis(fig[1,1],
              xlabel = "N",
              ylabel = "Phase advance [deg]",
@@ -589,7 +574,11 @@ end
 # ====== Nicholsplot ======
 # Note: This is a simplified version. The full implementation would require
 # porting all the complex grid calculations from the Plots version
-function CSMakie.nicholsplot(systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
+function CSMakie.nicholsplot(args...; kwargs...)
+    fig = Figure()
+    CSMakie.nicholsplot!(fig, args...; kwargs...)
+end
+function CSMakie.nicholsplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISystem}}, 
                              w=nothing; text=true, 
                              Gains=[12, 6, 3, 1, 0.5, -0.5, -1, -3, -6, -10, -20, -40, -60],
                              pInc=30, kwargs...)
@@ -597,7 +586,6 @@ function CSMakie.nicholsplot(systems::Union{LTISystem, AbstractVector{<:LTISyste
     systems, w = isnothing(w) ? _processfreqplot(Val{:nyquist}(), systems_vec) : 
                                 _processfreqplot(Val{:nyquist}(), systems_vec, w)
     
-    fig = Figure()
     ax = Axis(fig[1,1],
              title = "Nichols Chart",
              xlabel = "Phase [deg]",
@@ -621,7 +609,9 @@ function CSMakie.nicholsplot(systems::Union{LTISystem, AbstractVector{<:LTISyste
     return fig
 end
 
-function Makie.plot(r::SimResult; plotu=false, plotx=false, ploty=true)
+Makie.plot(r::SimResult; kwargs...) = Makie.plot!(Figure(), r; kwargs...)
+
+function Makie.plot!(fig, r::SimResult; plotu=false, plotx=false, ploty=true)
     ny, nu, nx = r.ny, r.nu, r.nx
     t = r.t
     n_series = size(r.y, 3)
@@ -632,7 +622,6 @@ function Makie.plot(r::SimResult; plotu=false, plotx=false, ploty=true)
     plotx && (nplots += nx)
     
     # Create figure with grid layout
-    fig = Figure()
     gl = GridLayout(fig[1, 1])
     
     plotind = 1
@@ -710,8 +699,9 @@ function Makie.plot(r::SimResult; plotu=false, plotx=false, ploty=true)
     return fig
 end
 
-function Makie.plot(si::StepInfo)
-    fig = Figure()
+Makie.plot(si::StepInfo; kwargs...) = Makie.plot!(Figure(), si; kwargs...)
+
+function Makie.plot!(fig, si::StepInfo)
     ax = Axis(fig[1,1],
              xlabel = "Time (s)",
              ylabel = "Output",
diff --git a/lib/ControlSystemsBase/src/CSMakie.jl b/lib/ControlSystemsBase/src/CSMakie.jl
index 9403faef0..f22157c65 100644
--- a/lib/ControlSystemsBase/src/CSMakie.jl
+++ b/lib/ControlSystemsBase/src/CSMakie.jl
@@ -18,135 +18,37 @@ module CSMakie
 
 using ..ControlSystemsBase: LTISystem
 
-# Error message for when Makie is not loaded
-const MAKIE_NOT_LOADED_ERROR = """
-Makie.jl must be loaded to use CSMakie plotting functions.
-Please run one of:
-  using GLMakie      # for interactive plots
-  using CairoMakie   # for static plots
-  using WGLMakie     # for web-based plots
-"""
-
-# Define stub functions that will be overloaded by the extension
-
-"""
-    bodeplot(sys::LTISystem, w=nothing; kwargs...)
-    bodeplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
-
-Create a Bode plot using Makie.jl. Requires Makie to be loaded.
-
-See the main `bodeplot` documentation for available options.
-"""
-function bodeplot(args...; kwargs...)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
-
-"""
-    nyquistplot(sys::LTISystem, w=nothing; kwargs...)
-    nyquistplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
-
-Create a Nyquist plot using Makie.jl. Requires Makie to be loaded.
-
-See the main `nyquistplot` documentation for available options.
-"""
-function nyquistplot(args...; kwargs...)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
-
-"""
-    sigmaplot(sys::LTISystem, w=nothing; kwargs...)
-    sigmaplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
-
-Create a singular value plot using Makie.jl. Requires Makie to be loaded.
-
-See the main `sigmaplot` documentation for available options.
-"""
-function sigmaplot(args...; kwargs...)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
-
-"""
-    marginplot(sys::LTISystem, w=nothing; kwargs...)
-    marginplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
-
-Create a margin plot using Makie.jl. Requires Makie to be loaded.
-
-See the main `marginplot` documentation for available options.
-"""
-function marginplot(args...; kwargs...)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
+function bodeplot end
+function bodeplot! end
 
-"""
-    pzmap(sys::LTISystem; kwargs...)
-    pzmap(sys::Vector{<:LTISystem}; kwargs...)
+function nyquistplot end
+function nyquistplot! end
 
-Create a pole-zero map using Makie.jl. Requires Makie to be loaded.
 
-See the main `pzmap` documentation for available options.
-"""
-function pzmap(args...; kwargs...)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
+function sigmaplot end
+function sigmaplot! end
 
-"""
-    pzmap!(ax, sys::LTISystem; kwargs...)
-    pzmap!(ax, sys::Vector{<:LTISystem}; kwargs...)
+function marginplot end
+function marginplot! end
 
-Add a pole-zero map to an existing Makie axis. Requires Makie to be loaded.
-"""
-function pzmap!(args...; kwargs...)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
+function pzmap end
+function pzmap! end
 
-"""
-    nicholsplot(sys::LTISystem, w=nothing; kwargs...)
-    nicholsplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
+function nicholsplot end
+function nicholsplot! end
 
-Create a Nichols chart using Makie.jl. Requires Makie to be loaded.
+function rgaplot end
+function rgaplot! end
 
-See the main `nicholsplot` documentation for available options.
-"""
-function nicholsplot(args...; kwargs...)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
-
-"""
-    rgaplot(sys::LTISystem, w=nothing; kwargs...)
-    rgaplot(sys::Vector{<:LTISystem}, w=nothing; kwargs...)
+function rlocusplot end
+function rlocusplot! end
 
-Create a relative gain array plot using Makie.jl. Requires Makie to be loaded.
-
-See the main `rgaplot` documentation for available options.
-"""
-function rgaplot(args...; kwargs...)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
-
-"""
-    rlocusplot(P::LTISystem, K=500; kwargs...)
-
-Create a root locus plot using Makie.jl. Requires Makie to be loaded.
-
-See the main `rlocusplot` documentation for available options.
-"""
-function rlocusplot(args...; kwargs...)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
-
-"""
-    leadlinkcurve(start=1; kwargs...)
-
-Plot the phase advance curve for a lead link using Makie.jl. Requires Makie to be loaded.
-
-See the main `leadlinkcurve` documentation for available options.
-"""
-function leadlinkcurve(args...; kwargs...)
-    error(MAKIE_NOT_LOADED_ERROR)
-end
+function leadlinkcurve end
+function leadlinkcurve! end
 
-# Also export the mutating versions that will be defined by the extension
-export bodeplot, nyquistplot, sigmaplot, marginplot, pzmap, pzmap!,
-       nicholsplot, rgaplot, rlocusplot, leadlinkcurve
+# Export all functions and their mutating versions
+export bodeplot, bodeplot!, nyquistplot, nyquistplot!, sigmaplot, sigmaplot!, 
+       marginplot, marginplot!, pzmap, pzmap!, nicholsplot, nicholsplot!, 
+       rgaplot, rgaplot!, rlocusplot, rlocusplot!, leadlinkcurve, leadlinkcurve!
 
 end # module CSMakie
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/test/test_makie_plots.jl b/lib/ControlSystemsBase/test/test_makie_plots.jl
index 9651d66f0..9fbc096f5 100644
--- a/lib/ControlSystemsBase/test/test_makie_plots.jl
+++ b/lib/ControlSystemsBase/test/test_makie_plots.jl
@@ -17,6 +17,7 @@ using GLMakie
     @testset "pzmap" begin
         @test_nowarn begin
             fig = CSMakie.pzmap(P)
+            CSMakie.pzmap!(fig[1,2], P)
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
@@ -32,6 +33,7 @@ using GLMakie
     @testset "bodeplot" begin
         @test_nowarn begin
             fig = CSMakie.bodeplot(P)
+            CSMakie.bodeplot!(fig[1,2], P)
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
@@ -58,6 +60,7 @@ using GLMakie
     @testset "nyquistplot" begin
         @test_nowarn begin
             fig = CSMakie.nyquistplot(P)
+            CSMakie.nyquistplot!(fig[1,2], P)
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
@@ -73,6 +76,7 @@ using GLMakie
     @testset "sigmaplot" begin
         @test_nowarn begin
             fig = CSMakie.sigmaplot(P)
+            CSMakie.sigmaplot!(fig[1,2], P)
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
@@ -84,6 +88,7 @@ using GLMakie
     @testset "marginplot" begin
         @test_nowarn begin
             fig = CSMakie.marginplot(P)
+            CSMakie.marginplot!(fig[1,2], P)
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
@@ -95,6 +100,7 @@ using GLMakie
     @testset "rlocusplot" begin
         @test_nowarn begin
             fig = CSMakie.rlocusplot(P)
+            CSMakie.rlocusplot!(fig[1,2], P)
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
@@ -106,6 +112,7 @@ using GLMakie
     @testset "rgaplot" begin
         @test_nowarn begin
             fig = CSMakie.rgaplot(Pmimo)
+            CSMakie.rgaplot!(fig[2,1:2], Pmimo)
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
@@ -120,10 +127,12 @@ using GLMakie
         res = step(P, t)
         @test_nowarn begin
             fig = plot(res)
+            plot!(fig[1,2], res)
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
             fig = plot(res; plotu=true, plotx=true)
+            plot!(fig[1,2], res, plotu=true, plotx=true)
             @test fig isa Makie.Figure
         end
     end
@@ -133,6 +142,7 @@ using GLMakie
         si = stepinfo(res)
         @test_nowarn begin
             fig = plot(si)
+            plot!(fig[1,2], si)
             @test fig isa Makie.Figure
         end
     end
@@ -140,6 +150,7 @@ using GLMakie
     @testset "leadlinkcurve" begin
         @test_nowarn begin
             fig = CSMakie.leadlinkcurve()
+            CSMakie.leadlinkcurve!(fig[1,2])
             @test fig isa Makie.Figure
         end
         @test_nowarn begin

From 03cbb431e523d455891dad7344db6c47b89c263e Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 25 Aug 2025 11:31:58 +0200
Subject: [PATCH 07/10] Switch to CairoMakie in tests

---
 docs/src/lib/plotting.md                      | 90 ++++++++++---------
 lib/ControlSystemsBase/Project.toml           |  4 +-
 lib/ControlSystemsBase/src/CSMakie.jl         |  4 +
 .../test/test_makie_plots.jl                  |  2 +-
 4 files changed, 53 insertions(+), 47 deletions(-)

diff --git a/docs/src/lib/plotting.md b/docs/src/lib/plotting.md
index faab520c3..a29d6221e 100644
--- a/docs/src/lib/plotting.md
+++ b/docs/src/lib/plotting.md
@@ -8,50 +8,6 @@ Pages = ["plotting.md"]
 !!! note "Time-domain responses"
     There are no special functions to plot time-domain results, such as step and impulse responses, instead, simply call `plot` on the result structure (`ControlSystemsBase.SimResult`) returned by [`lsim`](@ref), [`step`](@ref), [`impulse`](@ref) etc.
 
-## Makie support
-!!! danger "Experimental"
-    
-    The support for plotting with Makie is currently experimental and at any time subject to breaking changes or removal **not** respecting semantic versioning.
-
-ControlSystemsBase provides experimental support for plotting with [Makie.jl](https://docs.makie.org/) through the `CSMakie` module. This support is loaded automatically when you load a Makie backend (GLMakie, CairoMakie, or WGLMakie).
-
-### Usage
-
-```julia
-using ControlSystemsBase, GLMakie  # or CairoMakie, WGLMakie
-
-# Create a system
-P = tf([1], [1, 2, 1])
-
-# Use CSMakie plotting functions
-CSMakie.bodeplot(P)
-CSMakie.nyquistplot(P)
-CSMakie.pzmap(P)
-# ... and more
-
-# Direct plotting of simulation results
-res = step(P, 10)
-plot(res)  # Creates a figure with time-domain response
-
-si = stepinfo(res)
-plot(si)  # Visualizes step response characteristics
-```
-
-### Available functions
-
-The `CSMakie` module provides Makie implementations of the following plotting functions:
-
-- `CSMakie.bodeplot` - Bode magnitude and phase plots
-- `CSMakie.nyquistplot` - Nyquist plots with optional M and Mt circles
-- `CSMakie.sigmaplot` - Singular value plots
-- `CSMakie.marginplot` - Gain and phase margin plots
-- `CSMakie.pzmap` - Pole-zero maps
-- `CSMakie.nicholsplot` - Nichols charts
-- `CSMakie.rgaplot` - Relative gain array plots
-- `CSMakie.rlocusplot` - Root locus plots
-- `CSMakie.leadlinkcurve` - Lead-link design curves
-
-Additionally, `SimResult` and `StepInfo` types can be plotted directly using Makie's `plot` function.
 
 # Plotting functions
 
@@ -160,3 +116,49 @@ res = step(sysd, 5)
 plot(res, l=(:dash, 4))
 # plot!(stepinfo(step(sysd[1,1], 5))) # adds extra info to the plot
 ```
+
+
+## Makie support
+!!! danger "Experimental"
+    
+    The support for plotting with Makie is currently experimental and at any time subject to breaking changes or removal **not** respecting semantic versioning.
+
+ControlSystemsBase provides experimental support for plotting with [Makie.jl](https://docs.makie.org/) through the `CSMakie` module. This support is loaded automatically when you load a Makie backend (GLMakie, CairoMakie, or WGLMakie).
+
+### Usage
+
+```julia
+using ControlSystemsBase, GLMakie  # or CairoMakie, WGLMakie
+
+# Create a system
+P = tf([1], [1, 2, 1])
+
+# Use CSMakie plotting functions
+CSMakie.bodeplot(P)
+CSMakie.nyquistplot(P)
+CSMakie.pzmap(P)
+# ... and more
+
+# Direct plotting of simulation results
+res = step(P, 10)
+plot(res)  # Creates a figure with time-domain response
+
+si = stepinfo(res)
+plot(si)  # Visualizes step response characteristics
+```
+
+### Available functions
+
+The `CSMakie` module provides Makie implementations of the following plotting functions:
+
+- `CSMakie.bodeplot` - Bode magnitude and phase plots
+- `CSMakie.nyquistplot` - Nyquist plots with optional M and Mt circles
+- `CSMakie.sigmaplot` - Singular value plots
+- `CSMakie.marginplot` - Gain and phase margin plots
+- `CSMakie.pzmap` - Pole-zero maps
+- `CSMakie.nicholsplot` - Nichols charts
+- `CSMakie.rgaplot` - Relative gain array plots
+- `CSMakie.rlocusplot` - Root locus plots
+- `CSMakie.leadlinkcurve` - Lead-link design curves
+
+Additionally, `SimResult` and `StepInfo` types can be plotted directly using Makie's `plot` function.
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml
index 912034afa..6eb6a70a1 100644
--- a/lib/ControlSystemsBase/Project.toml
+++ b/lib/ControlSystemsBase/Project.toml
@@ -53,11 +53,11 @@ julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
-GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
 ImplicitDifferentiation = "57b37032-215b-411a-8a7c-41a003a55207"
@@ -68,4 +68,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Aqua", "ComponentArrays", "Documenter", "DSP", "FiniteDifferences", "ImplicitDifferentiation", "Makie", "GLMakie", "GenericSchur", "GR", "Plots", "SparseArrays", "StaticArrays"]
+test = ["Test", "Aqua", "ComponentArrays", "Documenter", "DSP", "FiniteDifferences", "ImplicitDifferentiation", "Makie", "CairoMakie", "GenericSchur", "GR", "Plots", "SparseArrays", "StaticArrays"]
diff --git a/lib/ControlSystemsBase/src/CSMakie.jl b/lib/ControlSystemsBase/src/CSMakie.jl
index f22157c65..a681afde0 100644
--- a/lib/ControlSystemsBase/src/CSMakie.jl
+++ b/lib/ControlSystemsBase/src/CSMakie.jl
@@ -13,6 +13,10 @@ CSMakie.nyquistplot(sys)
 ```
 
 All functions will throw an informative error if called without Makie.jl loaded.
+
+!!! danger "Experimental"
+
+    The support for plotting with Makie is currently experimental and at any time subject to breaking changes or removal **not** respecting semantic versioning.
 """
 module CSMakie
 
diff --git a/lib/ControlSystemsBase/test/test_makie_plots.jl b/lib/ControlSystemsBase/test/test_makie_plots.jl
index 9fbc096f5..40f3f15dd 100644
--- a/lib/ControlSystemsBase/test/test_makie_plots.jl
+++ b/lib/ControlSystemsBase/test/test_makie_plots.jl
@@ -1,7 +1,7 @@
 using ControlSystemsBase
 using Test
 using LinearAlgebra
-using GLMakie
+using CairoMakie
 
 @testset "Makie Plot Tests" begin
     # Create test systems

From e05440bb5fdbcef302ca6f4cd1a0824524a12be1 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 25 Aug 2025 12:10:46 +0200
Subject: [PATCH 08/10] handle imports nicer

---
 .../ext/ControlSystemsBaseMakieExt.jl             |  8 ++++----
 lib/ControlSystemsBase/test/test_makie_plots.jl   | 15 ++++++++-------
 2 files changed, 12 insertions(+), 11 deletions(-)

diff --git a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
index e98076ae9..5a3f5933c 100644
--- a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
+++ b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
@@ -151,9 +151,9 @@ function CSMakie.bodeplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISy
                 
                 if adaptive && eltype(magdata) <: AbstractFloat
                     wsi, magi, _ = downsample(ws, magdata, _span(magdata)/500)
-                    lines!(ax_mag, wsi, magi, label=lab)
+                    lines!(ax_mag, wsi, magi; label=lab, kwargs...)
                 else
-                    lines!(ax_mag, ws, magdata, label=lab)
+                    lines!(ax_mag, ws, magdata; label=lab, kwargs...)
                 end
                 
                 # Plot phase
@@ -173,9 +173,9 @@ function CSMakie.bodeplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTISy
                     if adaptive && eltype(phasedata) <: AbstractFloat
                         wsp, phasep, _ = downsample(ws, phasedata, 
                                                    _span(phasedata)/500)
-                        lines!(ax_phase, wsp, phasep)
+                        lines!(ax_phase, wsp, phasep; kwargs...)
                     else
-                        lines!(ax_phase, ws, phasedata)
+                        lines!(ax_phase, ws, phasedata; kwargs...)
                     end
                 end
             end
diff --git a/lib/ControlSystemsBase/test/test_makie_plots.jl b/lib/ControlSystemsBase/test/test_makie_plots.jl
index 40f3f15dd..2c420cfcd 100644
--- a/lib/ControlSystemsBase/test/test_makie_plots.jl
+++ b/lib/ControlSystemsBase/test/test_makie_plots.jl
@@ -1,7 +1,8 @@
 using ControlSystemsBase
 using Test
 using LinearAlgebra
-using CairoMakie
+import CairoMakie
+import CairoMakie.Makie
 
 @testset "Makie Plot Tests" begin
     # Create test systems
@@ -126,13 +127,13 @@ using CairoMakie
         t = 0:0.01:5
         res = step(P, t)
         @test_nowarn begin
-            fig = plot(res)
-            plot!(fig[1,2], res)
+            fig = Makie.plot(res)
+            Makie.plot!(fig[1,2], res)
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
-            fig = plot(res; plotu=true, plotx=true)
-            plot!(fig[1,2], res, plotu=true, plotx=true)
+            fig = Makie.plot(res; plotu=true, plotx=true)
+            Makie.plot!(fig[1,2], res, plotu=true, plotx=true)
             @test fig isa Makie.Figure
         end
     end
@@ -141,8 +142,8 @@ using CairoMakie
         res = step(P, 100)
         si = stepinfo(res)
         @test_nowarn begin
-            fig = plot(si)
-            plot!(fig[1,2], si)
+            fig = Makie.plot(si)
+            Makie.plot!(fig[1,2], si)
             @test fig isa Makie.Figure
         end
     end

From 3d41fbce670b3ce000736aac0144baf3ea62f625 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 25 Aug 2025 14:39:46 +0200
Subject: [PATCH 09/10] improve test coverage slightly

---
 .../ext/ControlSystemsBaseMakieExt.jl                | 12 ------------
 lib/ControlSystemsBase/test/test_makie_plots.jl      |  8 +++++++-
 2 files changed, 7 insertions(+), 13 deletions(-)

diff --git a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
index 5a3f5933c..dbb0217e6 100644
--- a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
+++ b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
@@ -372,18 +372,6 @@ function CSMakie.marginplot!(fig, systems::Union{LTISystem, AbstractVector{<:LTI
                                                                      allMargins=true, 
                                                                      adjust_phase_start)
                 
-                # Limit number of margins shown
-                if length(gm) > 5
-                    idx = sortperm(gm)
-                    wgm = wgm[idx[1:5]]
-                    gm = gm[idx[1:5]]
-                end
-                if length(pm) > 5
-                    idx = sortperm(pm)
-                    wpm = wpm[idx[1:5]]
-                    pm = pm[idx[1:5]]
-                end
-                
                 # Magnitude plot
                 ax_mag = axes_mag[i, j]
                 magdata = vec(bmag[i, j, :])
diff --git a/lib/ControlSystemsBase/test/test_makie_plots.jl b/lib/ControlSystemsBase/test/test_makie_plots.jl
index 2c420cfcd..88064ed55 100644
--- a/lib/ControlSystemsBase/test/test_makie_plots.jl
+++ b/lib/ControlSystemsBase/test/test_makie_plots.jl
@@ -65,7 +65,7 @@ import CairoMakie.Makie
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
-            fig = CSMakie.nyquistplot(P, w; Ms_circles=[1.5, 2.0])
+            fig = CSMakie.nyquistplot(P, w; Ms_circles=[1.5, 2.0], Mt_circles=[1.5, 2.0])
             @test fig isa Makie.Figure
         end
         @test_nowarn begin
@@ -136,6 +136,12 @@ import CairoMakie.Makie
             Makie.plot!(fig[1,2], res, plotu=true, plotx=true)
             @test fig isa Makie.Figure
         end
+
+        res = step(Pmimo, t)
+        @test_nowarn begin
+            fig = Makie.plot(res; plotu=true)
+            @test fig isa Makie.Figure
+        end
     end
     
     @testset "StepInfo plot" begin

From bd153dd67913a7440d9b49adbb489c492eb6e27c Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Mon, 25 Aug 2025 15:02:48 +0200
Subject: [PATCH 10/10] add legend to step plot

---
 .../ext/ControlSystemsBaseMakieExt.jl                  | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
index dbb0217e6..203f33620 100644
--- a/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
+++ b/lib/ControlSystemsBase/ext/ControlSystemsBaseMakieExt.jl
@@ -631,6 +631,10 @@ function Makie.plot!(fig, r::SimResult; plotu=false, plotx=false, ploty=true)
                     stairs!(ax, t, r.y[i, :, ms], label=label, step=:post)
                 end
             end
+            # Add legend if there are multiple series
+            if n_series > 1
+                axislegend(ax, position=:rt)
+            end
             plotind += 1
         end
     end
@@ -662,12 +666,14 @@ function Makie.plot!(fig, r::SimResult; plotu=false, plotx=false, ploty=true)
             
             if n_series > 1
                 for ms in 1:n_series
+                    label = "From $(input_names(r.sys)[ms])"
                     if iscontinuous(r.sys)
-                        lines!(ax, t, r.x[i, :, ms])
+                        lines!(ax, t, r.x[i, :, ms], label=label)
                     else
-                        stairs!(ax, t, r.x[i, :, ms], step=:post)
+                        stairs!(ax, t, r.x[i, :, ms], label=label, step=:post)
                     end
                 end
+                axislegend(ax, position=:rt)
             else
                 if iscontinuous(r.sys)
                     lines!(ax, t, r.x[i, :])