$\lambda$ variable for triaxial ellipsoids



In [1]:

    
%matplotlib inline



In [2]:

    
import matplotlib.pyplot as plt
import numpy as np

Here, we follow the reasoning presented by Webster (1904) for analyzing the ellipsoidal coordinate $\lambda$ describing a triaxial 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:

$$ \frac{x^{2}}{a^{2}} + \frac{y^{2}}{b^{2}} + \frac{z^{2}}{c^{2}} = 1 \: . \tag{1} $$

A quadric surface which is confocal with the ellipsoid defined in equation 1 can be described as follows:

$$ \frac{x^{2}}{a^{2} + \rho} + \frac{y^{2}}{b^{2} + \rho} + \frac{z^{2}}{c^{2} + \rho} = 1 \: , \tag{2} $$

where $\rho$ is a real number. We know that equation 2 represents an ellipsoid for $\rho$ satisfying the condition

$$ \rho + c^{2} > 0 \: . \tag{3} $$

Given $a$, $b$, $c$, 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$, $c$, $x$, $y$, and $z$, satisfies the equation 2.

By rearranging equation 2, we obtain the following cubic equation for $\rho$:

$$ f(\rho) = (a^{2} + \rho)(b^{2} + \rho)(c^{2} + \rho) - (b^{2} + \rho)(c^{2} + \rho) \, x^{2} - (a^{2} + \rho)(c^{2} + \rho) \, y^{2} - (a^{2} + \rho)(b^{2} + \rho) \, z^{2} \: . $$

This equation shows that:

$$ \rho = \begin{cases} d \to \infty \: &, \quad f(\rho) > 0 \\ -c^{2} \: &, \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:

$$ f(\rho) = \rho^{3} + p_{2} \, \rho^{2} + p_{1} \, \rho + p_{0} \: , \tag{4} $$

where

$$ p_{2} = a^{2} + b^{2} + c^{2} - x^{2} - y^{2} - z^{2} \: , \tag{5} $$

$$ p_{1} = b^{2} \, c^{2} + a^{2} \, c^{2} + a^{2} \, b^{2} - (b^{2} + c^{2}) \, x^{2} \ - (a^{2} + c^{2}) \, y^{2} - (a^{2} + b^{2}) \, z^{2} \tag{6} $$

and

$$ p_{0} = a^{2} \, b^{2} \, c^{2} - b^{2} \, c^{2} \, x^{2} - a^{2} \, c^{2} \, y^{2} - a^{2} \, b^{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$, $c$, $x$, $y$, and $z$ given in the cell below:



In [3]:

    
a = 200.
b = 180.
c = 150.
x = 210.
y = 230.
z = 300.

By using these constants, we calculate the coefficients $p_{2}$ (equation 5), $p_{1}$ (equation 6) and $p_{0}$ (equation 7) as follows:



In [4]:

    
p2 = a**2 + b**2 + c**2 - x**2 - y**2 - z**2
p1 = (b*c)**2 + (a*c)**2 + (a*b)**2 - (b**2 + c**2)*(x**2) - (a**2 + c**2)*(y**2) - (a**2 + b**2)*(z**2)
p0 = (a*b*c)**2 - (b*c*x)**2 - (a*c*y)**2 - (a*b*z)**2

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 cubic equation $f(\rho)$ (equation 4).



In [5]:

    
rho_min = -a**2 - 10.
rho_max = -c**2 + 2000.
rho = np.linspace(rho_min, rho_max, 100)



In [6]:

    
f = rho**3 + p2*(rho**2) + p1*rho + p0

Finally, the cell below shows the cubic equation $f(\rho)$ (equation 4) evaluated in the range $\left[ \rho_{min} \, , \rho_{max} \right]$ defined above.



In [7]:

    
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([-c**2, -c**2], [ymin, ymax], 'b--', label = '$-c^{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([-c**2, -c**2], [ymin, ymax], 'b--', label = '$-c^{2}$')
plt.plot(rho, f, 'k-', linewidth=2.)
plt.xlim(rho_min, -200.)
plt.ylim(-0.3*10**8, 10**7)
plt.legend(loc = 'best')
plt.xlabel('$\\rho$', fontsize = 20)
#plt.ylabel('$f(\\rho)$', fontsize = 20)

plt.tight_layout()

plt.show()

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 cubic equation $f(\rho)$ (equation 4).

The largest root $\lambda$ of $f(\rho)$ (equation 4) can be calculated as follows (Weisstein):

$$ \lambda = 2 \, \sqrt{-Q} \, \cos \left( \frac{\theta}{3}\right) - \frac{p_{2}}{3} \: , \tag{8} $$

where

$$ \theta = \cos^{-1} \left( \frac{R}{\sqrt{Q^{3}}} \right) \: , \tag{9} $$

$$ Q = \frac{3 \, p_{1} - p_{2}^{2}}{9} \tag{10a} $$

and

$$ R = \frac{9 \, p_{1} \, p_{2} - 27 \, p_{0} - 2 \, p_{2}^{3}}{54} \: , \tag{10b} $$

for the case in which the quantity $D < 0$, where $D = Q^{3} + R^{2}$. The cells below use the equations 8, 9, and 10 to compute the root $\lambda$.



In [8]:

    
p = (3.*p1 - p2**2)/3.
q = (9.*p1*p2 - 27.*p0 - 2.*p2**3)/27.

Q = p/3.
R = q/2.

D = Q**3 + R**2

print 'D = %.3f' % D









    



D = -1034963516680816233524756480.000



In [9]:

    
theta = np.arccos(R/np.sqrt(-Q*Q*Q))

lamb = 2.*np.sqrt(-Q)*np.cos(theta/3.) - p2/3.



In [10]:

    
print 'lambda = %.5f' % lamb









    



lambda = 157846.44992



In [14]:

    
np.roots([1.0, p2, p1, p0])









    Out[14]:





array([ 157846.44991877,  -37471.61650327,  -28274.8334155 ])

By substituing $\lambda$ in equation 4, we can verify that it is a root of $f(\rho)$.



In [11]:

    
f_lamb = lamb**3 + p2*(lamb**2) + p1*lamb + p0



In [12]:

    
print 'f(lambda) = %.5f' % f_lamb









    



f(lambda) = 0.00000

References

Webster, A. G. 1904. The Dynamics of Particles and of Rigid, Elastic and Fluid Bodies. Universidade de Michigan.

Weisstein, Eric W. Cubic Formula. From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/CubicFormula.html



In [ ]: