From d34361266b91555e541becf8d3c5e82d87bca6f8 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 30 May 2023 18:20:49 +0200
Subject: [PATCH] add option to draw an M circle in the Nyquist plot

---
 lib/ControlSystemsBase/src/plotting.jl    | 18 ++++++++++++++++--
 lib/ControlSystemsBase/test/test_plots.jl |  2 +-
 2 files changed, 17 insertions(+), 3 deletions(-)

diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl
index 4cf8c0ec1..78914e7bb 100644
--- a/lib/ControlSystemsBase/src/plotting.jl
+++ b/lib/ControlSystemsBase/src/plotting.jl
@@ -388,7 +388,7 @@ optionally provided.
 `kwargs` is sent as argument to plot.
 """
 nyquistplot
-@recipe function nyquistplot(p::Nyquistplot; Ms_circles=Float64[], Mt_circles=Float64[], unit_circle=false, hz=false, critical_point=-1, balance=true)
+@recipe function nyquistplot(p::Nyquistplot; Ms_circles=Float64[], Mt_circles=Float64[], M_circles=[], unit_circle=false, hz=false, critical_point=-1, balance=true)
     systems, w = _processfreqplot(Val{:nyquist}(), p.args...)
     ny, nu = size(systems[1])
     nw = length(w)
@@ -445,7 +445,21 @@ nyquistplot
                             rt = Mt/(Mt^2-1)    # Mt radius
                             ct.+rt.*C, rt.*S
                         end
-                    end                
+                    end  
+                    for M in M_circles
+                        @series begin
+                            subplot --> s2i(i,j)
+                            primary := false
+                            linestyle := :dash
+                            linecolor := :gray
+                            seriestype := :path
+                            markershape := :none
+                            label := "M = $(round(M, digits=2))"
+                            ct = -(2M^2 - 2M + 1)/(2M*(M-1)) # M center
+                            rt = (2M - 1)/(2M*(M-1)) # M radius
+                            ct.+rt.*C, rt.*S
+                        end
+                    end        
                     if unit_circle 
                         @series begin
                             subplot --> s2i(i,j)
diff --git a/lib/ControlSystemsBase/test/test_plots.jl b/lib/ControlSystemsBase/test/test_plots.jl
index 5b17b06f8..1d5274ffe 100644
--- a/lib/ControlSystemsBase/test/test_plots.jl
+++ b/lib/ControlSystemsBase/test/test_plots.jl
@@ -18,7 +18,7 @@ function getexamples()
       setPlotScale("dB")
       bodeplot(sys,ws)
     end
-    nyquistgen() = nyquistplot(sysss,ws, Ms_circles=1.2, Mt_circles=1.2)
+    nyquistgen() = nyquistplot(sysss,ws, Ms_circles=1.2, Mt_circles=1.2, M_circles=1.2)
     sigmagen() = sigmaplot(sysss,ws)
     #Only siso for now
     nicholsgen() = nicholsplot(tf1,ws)