<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 $\lambda$ variable for prolate ellipsoids
Step1:  Here, we follow the reasoning presented by Webster (1904) for analyzing the ellipsoidal coordinate  $\lambda$ describing a prolate ellipsoid.
Let's consider an ellipsoid with semi-axes $a$, $b$, $c$ oriented along the $x$-, $y$-, and $z$-axis, respectively, where $a > b = c > 0$. This ellipsoid is defined by the following equation
Step2:  By using these constants, we calculate the coefficients $p_{2}$ (equation 5), $p_{1}$ (equation 6) and $p_{0}$ (equation 7) as follows
Step3:  In the sequence, we define a set of values for the variable $\rho$ in an interval $\left[ \rho_{min} \, , \rho_{max} \right]$ and evaluate the quadratic equation $f(\rho)$ (equation 4).
Step4:  Finally, the cell below shows the quadratic equation $f(\rho)$ (equation 4) evaluated in the range $\left[ \rho_{min} \, , \rho_{max} \right]$ defined above.
Step5:  Remember that we are interested in a $\rho$ satisfying equation 3. Consequently, according to the figures shown above, we are interested in the largest root $\lambda$ of the quadratic equation $f(\rho)$ (equation 4).
The largest root $\lambda$ of $f(\rho)$ (equation 4) is given by
Step6:  By substituing $\lambda$ in equation 4, we can verify that it is a root of $f(\rho)$.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
Explanation: $\lambda$ variable for prolate ellipsoids
End of explanation
a = 20.
b = 13.
x = 21.
y = 23.
z = 30.
Explanation: Here, we follow the reasoning presented by Webster (1904) for analyzing the ellipsoidal coordinate  $\lambda$ describing a prolate ellipsoid.
Let's consider an ellipsoid with semi-axes $a$, $b$, $c$ oriented along the $x$-, $y$-, and $z$-axis, respectively, where $a > b = c > 0$. This ellipsoid is defined by the following equation:
<a id='eq1'></a>
$$
\frac{x^{2}}{a^{2}} + \frac{y^{2} + z^{2}}{b^{2}} = 1 \: . \tag{1}
$$
A quadric surface which is confocal with the ellipsoid defined in equation 1 can be described as follows:
<a id='eq2'></a>
$$
\frac{x^{2}}{a^{2} + \rho} + \frac{y^{2} + z^{2}}{b^{2} + \rho}= 1 \: , \tag{2}
$$
where $\rho$ is a real number. We know that equation 2 represents an ellipsoid for $\rho$ satisfying the condition
<a id='eq3'></a>
$$
\rho + b^{2} > 0 \: . \tag{3}
$$
Given $a$, $b$, and a $\rho$ satisfying equation 3, we may use equation 2 for determining a set of points $(x, y, z)$ lying on the surface of an ellipsoid confocal with that one defined in equation 1. Now, consider the problem of determining the ellipsoid which is confocal with that one defined in equation 1 and pass through a particular point $(x, y, z)$. This problem consists in determining the real number $\rho$ that, given $a$, $b$, $x$, $y$, and $z$, satisfies the equation 2.
By rearranging equation 2, we obtain the following quadratic equation for $\rho$:
$$
f(\rho) = (a^{2} + \rho)(b^{2} + \rho) - (b^{2} + \rho) \, x^{2}
- (a^{2} + \rho) \, (y^{2} + z^{2}) \: .
$$
This equation shows that:
$$
\rho = \begin{cases}
d \to \infty \: &, \quad f(\rho) > 0 \
-b^{2} \: &, \quad f(\rho) < 0 \
-a^{2} \: &, \quad f(\rho) > 0
\end{cases} \: .
$$
By rearanging this equation, we obtain a simpler one given by:
<a id='eq4'></a>
$$
f(\rho) = p_{2} \, \rho^{2} + p_{1} \, \rho + p_{0} \: , \tag{4}
$$
where
<a id='eq5'></a>
$$
p_{2} = 1 \: , \tag{5}
$$
<a id='eq6'></a>
$$
p_{1} =  a^{2} + b^{2} - x^{2} - y^{2} - z^{2} \tag{6}
$$
and
<a id='eq7'></a>
$$
p_{0} =  a^{2} \, b^{2} - b^{2} \, x^{2} - a^{2} \, y^{2} - a^{2} \, z^{2} \: . \tag{7}
$$
Note that a particular $\rho$ satisfying equation 2 results in $f(\rho) = 0$ (equation 4).
In order to illustrate the parameter $\rho$, consider the constants $a$, $b$, $x$, $y$, and $z$ given in the cell below:
End of explanation
p2 = 1.
p1 = a**2 + b**2 - (x**2) - (y**2) - (z**2)
p0 = (a*b)**2 - (b*x)**2 - (a*y)**2 - (a*z)**2
Explanation: By using these constants, we calculate the coefficients $p_{2}$ (equation 5), $p_{1}$ (equation 6) and $p_{0}$ (equation 7) as follows:
End of explanation
rho_min = -a**2 - 100.
rho_max = -b**2 + 2000.
rho = np.linspace(rho_min, rho_max, 100)
f = p2*(rho**2) + p1*rho + p0
Explanation: In the sequence, we define a set of values for the variable $\rho$ in an interval $\left[ \rho_{min} \, , \rho_{max} \right]$ and evaluate the quadratic equation $f(\rho)$ (equation 4).
End of explanation
ymin = np.min(f) - 0.1*(np.max(f) - np.min(f))
ymax = np.max(f) + 0.1*(np.max(f) - np.min(f))
plt.close('all')
plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.plot([rho_min, rho_max], [0., 0.], 'k-')
plt.plot([-a**2, -a**2], [ymin, ymax], 'r--', label = '$-a^{2}$')
plt.plot([-b**2, -b**2], [ymin, ymax], 'g--', label = '$-b^{2}$')
plt.plot(rho, f, 'k-', linewidth=2.)
plt.xlim(rho_min, rho_max)
plt.ylim(ymin, ymax)
plt.legend(loc = 'best')
plt.xlabel('$\\rho$', fontsize = 20)
plt.ylabel('$f(\\rho)$', fontsize = 20)
plt.subplot(1,2,2)
plt.plot([rho_min, rho_max], [0., 0.], 'k-')
plt.plot([-a**2, -a**2], [ymin, ymax], 'r--', label = '$-a^{2}$')
plt.plot([-b**2, -b**2], [ymin, ymax], 'g--', label = '$-b^{2}$')
plt.plot(rho, f, 'k-', linewidth=2.)
plt.xlim(rho_min, 0.)
plt.ylim(-0.5*10**6, 10**5)
plt.legend(loc = 'best')
plt.xlabel('$\\rho$', fontsize = 20)
#plt.ylabel('$f(\\rho)$', fontsize = 20)
plt.tight_layout()
plt.show()
Explanation: Finally, the cell below shows the quadratic equation $f(\rho)$ (equation 4) evaluated in the range $\left[ \rho_{min} \, , \rho_{max} \right]$ defined above.
End of explanation
delta = p1**2 - 4.*p2*p0
lamb = (-p1 + np.sqrt(delta))/(2.*p2)
print 'lambda = %.5f' % lamb
Explanation: Remember that we are interested in a $\rho$ satisfying equation 3. Consequently, according to the figures shown above, we are interested in the largest root $\lambda$ of the quadratic equation $f(\rho)$ (equation 4).
The largest root $\lambda$ of $f(\rho)$ (equation 4) is given by:
<a id='eq8'></a>
$$
\lambda = \frac{-p_{1} + \sqrt{\Delta}}{2} \: , \tag{8}
$$
where
<a id='eq9'></a>
$$
\Delta = p_{1}^{2} - 4 \, p_{0} \: . \tag{9}
$$
The cells below use the equations 8 and 9 to compute the root $\lambda$.
End of explanation
f_lamb = p2*(lamb**2) + p1*lamb + p0
print 'f(lambda) = %.5f' % f_lamb
Explanation: By substituing $\lambda$ in equation 4, we can verify that it is a root of $f(\rho)$.
End of explanation
<END_TASK>