from numpy import linspace
import matplotlib.pyplot as plt
from numpy import sqrt


def egg_radius(x, L, B, w, D):
    # Note that k does not involve x
    k = sqrt(5.5 * L**2 + 11 * L * w + 4 * w**2)
    k *= sqrt(3) * B * L - 2 * D * sqrt(L**2 + 2 * w * L + 4 * w**2)
    k /= sqrt(3) * B * L
    k /= sqrt(5.5 * L**2 + 11 * L * w + 4 * w * w) - 2 * sqrt(
        L**2 + 2 * w * L + 4 * w**2
    )

    # coefficients of a quadratic in x
    a = 2 * (L - 2 * w)
    b = L**2 + 8 * L * w - 4 * w**2
    c = 2 * L * w**2 + L**2 * w + L**3

    t = L * (L**2 + 8 * w * x + 4 * w**2)
    t /= a * x**2 + b * x + c

    y = 0.5 * B * sqrt(L**2 - 4 * x**2) / sqrt(L**2 + 8 * w * x + 4 * w**2)
    y *= 1 - k * (1 - sqrt(t))

    return y


L = 2
B = 1.3
w = 0.2
D = 0.9

plt.plot([-w, -w], [-B / 2, B / 2], "g--")
z = egg_radius(L / 4, L, B, w, D)
plt.plot([L / 4, L / 4], [-z, z], "r-.")
x = linspace(-L / 2, L / 2, 200)
y = egg_radius(x, L, B, w, D)
plt.plot(x, y, "b")
plt.plot(x, -y, "b")
plt.legend(["B", "D"])
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.show()