import numpy as np
import matplotlib.pyplot as plt

M_values = np.linspace(1.01, 10, 20)

fig, ax = plt.subplots(1, 1, figsize=(10, 10))

for M in M_values:
    x = np.linspace((1-M)/(1+M)+1e-6, 1-1e-6, 1000)

    y = -(M**4 * (x**2 - 1) - 2 * M**2 * ((x - 1) * x + 1) + 2 * np.sqrt(M**4 * (-2 * M**2 * (x - 1) + 2 * x - 1)) + (x - 1)**2)/(M**2 - 1)**2
    y = np.sqrt(y)
    ax.plot(x, y, color='blue', alpha=0.5)
    ax.plot(x, -y, color='blue', alpha=0.5)
    xmin = 0.4
    ax.plot([xmin, 1], [(x-1) * np.tan(np.arccos(1/M)) for x in [xmin, 1]], 'r-', alpha=0.2)
    ax.plot([xmin, 1], [-(x-1) * np.tan(np.arccos(1/M)) for x in [xmin, 1]], 'r-', alpha=0.2)

plt.gca().add_patch(plt.Circle((0, 0), 1, fill=False))

plt.xlim(-1.01,1.01)
plt.ylim(-1.01,1.01)
plt.show()
plt.savefig('plot.svg', format='svg')