python-control · murrayrm · Feb 16, 2025 · Feb 1, 2025 · Feb 2, 2025 · Feb 4, 2025
diff --git a/control/phaseplot.py b/control/phaseplot.py
diff --git a/control/tests/ctrlplot_test.py b/control/tests/ctrlplot_test.py
@@ -116,6 +116,11 @@ def setup_plot_arguments(resp_fcn, plot_fcn, compute_time_response=True):
                 args2 = (sys2, )
                 argsc = ([sys1, sys2], )
 
+        case (None, ct.phase_plane_plot):
+            args1 = (sys1, )
+            args2 = (sys2, )
+            plot_kwargs = {'plot_streamlines': True}
+
         case _, _:
             args1 = (sys1, )
             args2 = (sys2, )

diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
@@ -172,6 +172,7 @@ def test_unrecognized_kwargs(function, nsssys, ntfsys, moreargs, kwargs,
      (control.phase_plane_plot, 1, ([-1, 1, -1, 1], 1), {}),
      (control.phaseplot.streamlines, 1, ([-1, 1, -1, 1], 1), {}),
      (control.phaseplot.vectorfield, 1, ([-1, 1, -1, 1], ), {}),
+     (control.phaseplot.streamplot, 1, ([-1, 1, -1, 1], ), {}),
      (control.phaseplot.equilpoints, 1, ([-1, 1, -1, 1], ), {}),
      (control.phaseplot.separatrices, 1, ([-1, 1, -1, 1], ), {}),
      (control.singular_values_plot, 1, (), {})]
@@ -360,6 +361,7 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
         optimal_test.test_oep_argument_errors,
     'phaseplot.streamlines': test_matplotlib_kwargs,
     'phaseplot.vectorfield': test_matplotlib_kwargs,
+    'phaseplot.streamplot': test_matplotlib_kwargs,
     'phaseplot.equilpoints': test_matplotlib_kwargs,
     'phaseplot.separatrices': test_matplotlib_kwargs,
 }

diff --git a/control/tests/phaseplot_test.py b/control/tests/phaseplot_test.py
@@ -12,6 +12,7 @@
 import warnings
 from math import pi
 
+import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 import pytest
@@ -138,29 +139,39 @@ def invpend_ode(t, x, m=0, l=0, b=0, g=0):
 
     # Use callable form, with parameters (if not correct, will get /0 error)
     ct.phase_plane_plot(
-        invpend_ode, [-5, 5, -2, 2], params={'args': (1, 1, 0.2, 1)})
+        invpend_ode, [-5, 5, -2, 2], params={'args': (1, 1, 0.2, 1)},
+        plot_streamlines=True)
 
     # Linear I/O system
     ct.phase_plane_plot(
-        ct.ss([[0, 1], [-1, -1]], [[0], [1]], [[1, 0]], 0))
+        ct.ss([[0, 1], [-1, -1]], [[0], [1]], [[1, 0]], 0),
+        plot_streamlines=True)
 
 
 @pytest.mark.usefixtures('mplcleanup')
 def test_phaseplane_errors():
     with pytest.raises(ValueError, match="invalid grid specification"):
-        ct.phase_plane_plot(ct.rss(2, 1, 1), gridspec='bad')
+        ct.phase_plane_plot(ct.rss(2, 1, 1), gridspec='bad',
+                            plot_streamlines=True)
 
     with pytest.raises(ValueError, match="unknown grid type"):
-        ct.phase_plane_plot(ct.rss(2, 1, 1), gridtype='bad')
+        ct.phase_plane_plot(ct.rss(2, 1, 1), gridtype='bad',
+                            plot_streamlines=True)
 
     with pytest.raises(ValueError, match="system must be planar"):
-        ct.phase_plane_plot(ct.rss(3, 1, 1))
+        ct.phase_plane_plot(ct.rss(3, 1, 1),
+                            plot_streamlines=True)
 
     with pytest.raises(ValueError, match="params must be dict with key"):
         def invpend_ode(t, x, m=0, l=0, b=0, g=0):
             return (x[1], -b/m*x[1] + (g*l/m) * np.sin(x[0]))
         ct.phase_plane_plot(
-            invpend_ode, [-5, 5, 2, 2], params={'stuff': (1, 1, 0.2, 1)})
+            invpend_ode, [-5, 5, 2, 2], params={'stuff': (1, 1, 0.2, 1)},
+                            plot_streamlines=True)
+
+    with pytest.raises(ValueError, match="gridtype must be 'meshgrid' when using streamplot"):
+        ct.phase_plane_plot(ct.rss(2, 1, 1), plot_streamlines=False,
+                            plot_streamplot=True, gridtype='boxgrid')
 
     # Warning messages for invalid solutions: nonlinear spring mass system
     sys = ct.nlsys(
@@ -170,14 +181,87 @@ def invpend_ode(t, x, m=0, l=0, b=0, g=0):
     with pytest.warns(UserWarning, match=r"X0=array\(.*\), solve_ivp failed"):
         ct.phase_plane_plot(
             sys, [-12, 12, -10, 10], 15, gridspec=[2, 9],
-            plot_separatrices=False)
+            plot_separatrices=False, plot_streamlines=True)
 
     # Turn warnings off
     with warnings.catch_warnings():
         warnings.simplefilter("error")
         ct.phase_plane_plot(
             sys, [-12, 12, -10, 10], 15, gridspec=[2, 9],
-            plot_separatrices=False, suppress_warnings=True)
+            plot_streamlines=True, plot_separatrices=False,
+            suppress_warnings=True)
+
+@pytest.mark.usefixtures('mplcleanup')
+def test_phase_plot_zorder():
+    # some of these tests are a bit akward since the streamlines and separatrices
+    # are stored in the same list, so we separate them by color
+    key_color = "tab:blue" # must not be 'k', 'r', 'b' since they are used by separatrices
+
+    def get_zorders(cplt):
+        max_zorder = lambda items: max([line.get_zorder() for line in items])
+        assert isinstance(cplt.lines[0], list)
+        streamline_lines = [line for line in cplt.lines[0] if line.get_color() == key_color]
+        separatrice_lines = [line for line in cplt.lines[0] if line.get_color() != key_color]
+        streamlines = max_zorder(streamline_lines) if streamline_lines else None
+        separatrices = max_zorder(separatrice_lines) if separatrice_lines else None
+        assert cplt.lines[1] == None or isinstance(cplt.lines[1], mpl.quiver.Quiver)
+        quiver = cplt.lines[1].get_zorder() if cplt.lines[1] else None
+        assert cplt.lines[2] == None or isinstance(cplt.lines[2], list)
+        equilpoints = max_zorder(cplt.lines[2]) if cplt.lines[2] else None
+        assert cplt.lines[3] == None or isinstance(cplt.lines[3], mpl.streamplot.StreamplotSet)
+        streamplot = max(cplt.lines[3].lines.get_zorder(), cplt.lines[3].arrows.get_zorder()) if cplt.lines[3] else None
+        return streamlines, quiver, streamplot, separatrices, equilpoints
+
+    def assert_orders(streamlines, quiver, streamplot, separatrices, equilpoints):
+        print(streamlines, quiver, streamplot, separatrices, equilpoints)
+        if streamlines is not None:
+            assert streamlines < separatrices < equilpoints
+        if quiver is not None:
+            assert quiver < separatrices < equilpoints
+        if streamplot is not None:
+            assert streamplot < separatrices < equilpoints
+
+    def sys(t, x):
+        return np.array([4*x[1], -np.sin(4*x[0])])
+
+    # ensure correct zordering for all three flow types
+    res_streamlines = ct.phase_plane_plot(sys, plot_streamlines=dict(color=key_color))
+    assert_orders(*get_zorders(res_streamlines))
+    res_vectorfield = ct.phase_plane_plot(sys, plot_vectorfield=True)
+    assert_orders(*get_zorders(res_vectorfield))
+    res_streamplot = ct.phase_plane_plot(sys, plot_streamplot=True)
+    assert_orders(*get_zorders(res_streamplot))
+
+    # ensure that zorder can still be overwritten
+    res_reversed = ct.phase_plane_plot(sys, plot_streamlines=dict(color=key_color, zorder=50), plot_vectorfield=dict(zorder=40),
+                                       plot_streamplot=dict(zorder=30), plot_separatrices=dict(zorder=20), plot_equilpoints=dict(zorder=10))
+    streamlines, quiver, streamplot, separatrices, equilpoints = get_zorders(res_reversed)
+    assert streamlines > quiver > streamplot > separatrices > equilpoints
+
+
+@pytest.mark.usefixtures('mplcleanup')
+def test_stream_plot_magnitude():
+    def sys(t, x):
+        return np.array([4*x[1], -np.sin(4*x[0])])
+
+    # plt context with linewidth
+    with plt.rc_context({'lines.linewidth': 4}):
+        res = ct.phase_plane_plot(sys, plot_streamplot=dict(vary_linewidth=True))
+    linewidths = res.lines[3].lines.get_linewidths()
+    # linewidths are scaled to be between 0.25 and 2 times default linewidth
+    # but the extremes may not exist if there is no line at that point
+    assert min(linewidths) < 2 and max(linewidths) > 7
+
+    # make sure changing the colormap works
+    res = ct.phase_plane_plot(sys, plot_streamplot=dict(vary_color=True, cmap='viridis'))
+    assert res.lines[3].lines.get_cmap().name == 'viridis'
+    res = ct.phase_plane_plot(sys, plot_streamplot=dict(vary_color=True, cmap='turbo'))
+    assert res.lines[3].lines.get_cmap().name == 'turbo'
+
+    # make sure changing the norm at least doesn't throw an error
+    ct.phase_plane_plot(sys, plot_streamplot=dict(vary_color=True, norm=mpl.colors.LogNorm()))
+
+
 
 
 @pytest.mark.usefixtures('mplcleanup')
@@ -189,7 +273,7 @@ def test_basic_phase_plots(savefigs=False):
     plt.figure()
     axis_limits = [-1, 1, -1, 1]
     T = 8
-    ct.phase_plane_plot(sys, axis_limits, T)
+    ct.phase_plane_plot(sys, axis_limits, T, plot_streamlines=True)
     if savefigs:
         plt.savefig('phaseplot-dampedosc-default.png')
 
@@ -202,7 +286,7 @@ def invpend_update(t, x, u, params):
     ct.phase_plane_plot(
         invpend, [-2*pi, 2*pi, -2, 2], 5,
         gridtype='meshgrid', gridspec=[5, 8], arrows=3,
-        plot_separatrices={'gridspec': [12, 9]},
+        plot_separatrices={'gridspec': [12, 9]}, plot_streamlines=True,
         params={'m': 1, 'l': 1, 'b': 0.2, 'g': 1})
     plt.xlabel(r"$\theta$ [rad]")
     plt.ylabel(r"$\dot\theta$ [rad/sec]")
@@ -217,7 +301,8 @@ def oscillator_update(t, x, u, params):
         oscillator_update, states=2, inputs=0, name='nonlinear oscillator')
 
     plt.figure()
-    ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9)
+    ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9,
+                        plot_streamlines=True)
     pp.streamlines(
     oscillator, np.array([[0, 0]]), 1.5,
     gridtype='circlegrid', gridspec=[0.5, 6], dir='both')
@@ -227,6 +312,18 @@ def oscillator_update(t, x, u, params):
     if savefigs:
         plt.savefig('phaseplot-oscillator-helpers.png')
 
+    plt.figure()
+    ct.phase_plane_plot(
+        invpend, [-2*pi, 2*pi, -2, 2],
+        plot_streamplot=dict(vary_color=True, vary_density=True),
+        gridspec=[60, 20], params={'m': 1, 'l': 1, 'b': 0.2, 'g': 1}
+    )
+    plt.xlabel(r"$\theta$ [rad]")
+    plt.ylabel(r"$\dot\theta$ [rad/sec]")
+
+    if savefigs:
+        plt.savefig('phaseplot-invpend-streamplot.png')
+
 
 if __name__ == "__main__":
     #

diff --git a/doc/figures/phaseplot-dampedosc-default.png b/doc/figures/phaseplot-dampedosc-default.png
diff --git a/doc/figures/phaseplot-invpend-meshgrid.png b/doc/figures/phaseplot-invpend-meshgrid.png
diff --git a/doc/figures/phaseplot-oscillator-helpers.png b/doc/figures/phaseplot-oscillator-helpers.png
diff --git a/doc/functions.rst b/doc/functions.rst
@@ -103,6 +103,7 @@ Phase plane plots
     phaseplot.separatrices
     phaseplot.streamlines
     phaseplot.vectorfield
+    phaseplot.streamplot
 
 
 Frequency Response

diff --git a/doc/phaseplot.rst b/doc/phaseplot.rst
@@ -12,7 +12,7 @@ functionality is supported by a set of mapping functions that are part of
 the `phaseplot` module.
 
 The default method for generating a phase plane plot is to provide a
-2D dynamical system along with a range of coordinates and time limit:
+2D dynamical system along with a range of coordinates in phase space:
 
 .. testsetup:: phaseplot
 
@@ -27,8 +27,7 @@ The default method for generating a phase plane plot is to provide a
         sys_update, states=['position', 'velocity'],
 	inputs=0, name='damped oscillator')
     axis_limits = [-1, 1, -1, 1]
-    T = 8
-    ct.phase_plane_plot(sys, axis_limits, T)
+    ct.phase_plane_plot(sys, axis_limits)
 
 .. testcode:: phaseplot
     :hide:
@@ -39,12 +38,12 @@ The default method for generating a phase plane plot is to provide a
 .. image:: figures/phaseplot-dampedosc-default.png
    :align: center
 
-By default, the plot includes streamlines generated from starting
-points on limits of the plot, with arrows showing the flow of the
-system, as well as any equilibrium points for the system.  A variety
+By default the plot includes streamlines infered from function values
+on a grid, equilibrium points and separatrices if they exist. A variety
 of options are available to modify the information that is plotted,
-including plotting a grid of vectors instead of streamlines and
-turning on and off various features of the plot.
+including plotting a grid of vectors instead of streamlines, plotting
+streamlines from arbitrary starting points and turning on and off
+various features of the plot. 
 
 To illustrate some of these possibilities, consider a phase plane plot for
 an inverted pendulum system, which is created using a mesh grid:
@@ -62,9 +61,7 @@ an inverted pendulum system, which is created using a mesh grid:
     invpend = ct.nlsys(invpend_update, states=2, inputs=1, name='invpend')
 
     ct.phase_plane_plot(
-        invpend, [-2 * np.pi, 2 * np.pi, -2, 2], 5,
-        gridtype='meshgrid', gridspec=[5, 8], arrows=3,
-        plot_equilpoints={'gridspec': [12, 9]},
+        invpend, [-2 * np.pi, 2 * np.pi, -2, 2],
         params={'m': 1, 'l': 1, 'b': 0.2, 'g': 1})
     plt.xlabel(r"$\theta$ [rad]")
     plt.ylabel(r"$\dot\theta$ [rad/sec]")
@@ -79,16 +76,17 @@ an inverted pendulum system, which is created using a mesh grid:
 
 This figure shows several features of more complex phase plane plots:
 multiple equilibrium points are shown, with saddle points showing
-separatrices, and streamlines generated along a 5x8 mesh of initial
-conditions.  At each mesh point, a streamline is created that goes 5 time
-units forward and backward in time.  A separate grid specification is used
-to find equilibrium points and separatrices (since the course grid spacing
-of 5x8 does not find all possible equilibrium points).  Together, the
-multiple features in the phase plane plot give a good global picture of the
+separatrices, and streamlines generated generated from a rectangular
+25x25 grid (default) of function evaluations. Together, the multiple
+features in the phase plane plot give a good global picture of the
 topological structure of solutions of the dynamical system.
 
-Phase plots can be built up by hand using a variety of helper functions that
-are part of the :mod:`phaseplot` (pp) module:
+Phase plots can be built up by hand using a variety of helper
+functions that are part of the :mod:`phaseplot` (pp) module. For more
+precise control, the streamlines can also generated by integrating the
+system forwards or backwards in time from a set of initial
+conditions. The initial conditions can be chosen on a rectangular
+grid, rectangual boundary, circle or from an arbitrary set of points.
 
 .. testcode:: phaseplot
     :hide:
@@ -105,7 +103,8 @@ are part of the :mod:`phaseplot` (pp) module:
     oscillator = ct.nlsys(
         oscillator_update, states=2, inputs=0, name='nonlinear oscillator')
 
-    ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9)
+    ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9,
+        plot_streamlines=True)
     pp.streamlines(
         oscillator, np.array([[0, 0]]), 1.5,
         gridtype='circlegrid', gridspec=[0.5, 6], dir='both')
@@ -128,6 +127,7 @@ The following helper functions are available:
    phaseplot.equilpoints
    phaseplot.separatrices
    phaseplot.streamlines
+   phaseplot.streamplot
    phaseplot.vectorfield
 
 The :func:`phase_plane_plot` function calls these helper functions