import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

A = np.array([[8, 16], [20, 23j]])

# The following code plots 500 random points in the numerical range of a matrix formed from the numbers in today’s date. The eigenvalues are plotted in black.

for _ in range(5000):
    v = norm(0, 1).rvs(size=4)
    v /= np.linalg.norm(v)
    x = np.array([v[0] + 1j * v[1], v[2] + 1j * v[3]])
    z = np.conj(x).dot(A.dot(x))
    plt.plot(z.real, z.imag, "o", color="0.9")

eigenvalues, eigenvectors = np.linalg.eig(A)
plt.plot(eigenvalues[0].real, eigenvalues[0].imag, "ko")
plt.plot(eigenvalues[1].real, eigenvalues[1].imag, "ko")


# The following code draws the ellipse guaranteed by the elliptical range theorem.
A_star_A = np.dot(np.conj(A).T, A)
b = (np.trace(A_star_A) - abs(eigenvalues[0]) ** 2 - abs(eigenvalues[1]) ** 2) ** 0.5
c = 0.5 * abs(eigenvalues[0] - eigenvalues[1])
a = (2 / 3) * (2 * (2 * b**2 / 4 + c**2) ** 0.5 + c)


t = np.linspace(0, np.pi * 2)
z = a * np.cos(t) / 2 + 1j * b * np.sin(t) / 2
u = (eigenvalues[1] - eigenvalues[0]) / 2
u /= np.linalg.norm(u)
m = eigenvalues.mean()
z = z * u + m

plt.plot(z.real, z.imag, "b-")
plt.show()


#  Ulrich Daepp, Pamela Gorkin, Andrew Shaffer, and Karl Voss. Finding Ellipses: What Blaschke Products, Poncelet’s Theorem, and the Numerical Range Know about Each Other. MAA Press 2018.