I try to find points on the surface of an ellipsoid that are also within a conical frustum.
The ellipsoid and frustum are centred at ce
and cf
, respectively.
The ellipsoid has semi-axes lengths defined in l
.
The frustum had a radius r1
at its base, an opening div
, and a height h
.
To solve the issue at hand, I use python
, together with numpy
and scipy
.
More specifically, I use the fsolve
function from scipy.optimize
to minimise the difference between the ellipsoid and frustum equations.
Based on that function, I have a Jacobian.
My code is
import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def get_ellispoid(c, l):
X_ell = c[0] + l[0] * np.outer(np.sin(Phi), np.cos(Theta))
Y_ell = c[1] + l[1] * np.outer(np.sin(Phi), np.sin(Theta))
Z_ell = c[2] + l[2] * np.outer(np.cos(Phi), np.ones_like(Theta))
return X_ell, Y_ell, Z_ell
def get_conical_frustum(c: tuple[float, float, float],
r1: float,
alpha: float,
z: np.ndarray,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
Theta, Z = np.meshgrid(theta, z)
R = r1 + Z * np.tan(alpha)
X_fru = c[0] + R * np.cos(Theta)
Y_fru = c[1] + R * np.sin(Theta)
Z_fru = c[2] + Z
return X_fru, Y_fru, Z_fru
def ellipsoid(theta, phi, ce, l):
return np.array(
[ce[0] + l[0] * np.cos(theta) * np.cos(phi),
ce[1] + l[1] * np.sin(theta) * np.cos(phi),
ce[2] + l[2] * np.sin(phi)]
)
def frustum(theta, z, cf, r1, div):
return np.array(
[cf[0] + (r1 + z * np.tan(div)) * np.cos(theta),
cf[1] + (r1 + z * np.tan(div)) * np.sin(theta),
cf[2] + z]
)
# Define the function for which we want to find the roots
def func(x, theta, phi, z, ce, cf, l, r1, div):
return ellipsoid(theta, phi, ce, l) - frustum(theta, z, cf, r1, div)
# Define the Jacobian matrix of the function
def jacobian(x, theta, phi, z, ce, cf, l, r1, div):
return np.array(
[
[-(l[0] * np.cos(phi) - r1 - z * np.tan(div)) * np.sin(theta), -l[0] * np.cos(theta) * np.sin(phi), -np.tan(div) * np.cos(theta)],
[-(l[1] * np.cos(phi) - r1 - z * np.tan(div)) * np.cos(theta), -l[1] * np.cos(theta) * np.sin(phi), -np.tan(div) * np.sin(theta)],
[0.0, l[2] * np.cos(phi), -1],
]
)
num_points = 20
theta = np.linspace(0, 2*np.pi, num_points)
phi = np.linspace(0, np.pi, num_points)
Theta, Phi = np.meshgrid(theta, phi)
# Ellipsoid
ce = [0.0, 0.0, 0.0]
l = [1.5, 1.1, 0.45]
# Conical frustum
cf = [0.0, 0.0, -1.5]
h = np.linalg.norm(np.array(cf) - np.array(ce))
z = np.linspace(cf[-1], h, num_points)
r1 = 0.15
div = np.pi / 18
# Initial guess
x0 = np.array([0.0, 0.0, ce[-1]])
# Iterate over theta, phi, and z values
sol = []
for it, t in enumerate(theta):
for ip, p in enumerate(phi):
for iz, z_val in enumerate(z):
try:
root = fsolve(func, [t, p, z_val], args=(t, p, z_val, ce, cf, l, r1, div), fprime=jacobian)
sol.append(root)
except RuntimeError:
continue
sol = np.array(sol)
# Create a 3D plot
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection = '3d')
# Plot the volumes
xell, yell, zell = get_ellispoid(ce, l)
ax.plot_surface(xell, yell, zell, color = 'r', alpha = 0.005)
xfru, yfru, zfru = get_conical_frustum(cf, r1, div, z)
ax.plot_surface(xfru, yfru, zfru, color = 'b', alpha = 0.35)
# Plot the intersections
ax.scatter(sol[:, 0], sol[:, 1], sol[:, 2], color = 'g', s = 1.5)
# Set labels and title
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
# Show the plot
plt.show()
When the code is run, the points to which the code converges are the initial condition [t, p, z_val]
or, if I try to put x0
as initial guess, the code converges to that value (which is not a valid solution), repeated over and over again.
The solutions should be all the points of the frustum lying on the ellipsoid surface.
I do not get why my code does not behave as expected?
Can someone help me find out?
There are a number of problems with your approach.
The intersection of the two surfaces is a curve. A curve in space has just one parameter, so you should not be using a triple-nested loop to do anything.
You appear to be under the impression that the polar angle theta in the ellipsoid and the frustrum is the same. It is not. Theta and phi in the ellipsoid are just parameters: they are not actually the spherical polar angles (unless the ellipsoid happens to reduce to a sphere).
There may be no intersections or there may in fact be two intersection curves, depending on the height of the frustrum. fsolve would never be able to find both (or even indicate which it had found).
You can solve the problem exactly. For each polar angle phi for the generators of the frustrum there may be 0, 1 or 2 values of z, according to the number of intersections. These values are simply the real roots of a quadratic equation, which may be set up for each value of phi - see the bottom of this post.
I give a code solution below and a sketch proof of the setting up of the quadratic equations at the bottom of this post. Note that I do not know whether you are referring to r1 as the larger or smaller radius of the frustrum. I have taken it as the larger one, but had to increase it to get any intersection. I do not use fsolve at all: there is no need and it could not solve the problem.
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Ellipsoid
xe, ye, ze = 0.0, 0.0, 0.0
a, b, c = 1.5, 1.1, 0.45
# Conical frustum
xf, yf, zf = 0.0, 0.0, -1.5
h = 1.5
R = 0.45
alpha = math.pi / 18
tanalpha = math.tan( alpha )
angles = np.linspace( 0.0, 2 * np.pi, 40 )
upper = []
lower = []
for phi in angles:
# Set up coefficients in quadratic equation
cosphi, sinphi = math.cos( phi ), math.sin( phi )
x1, x2 = xf - xe + ( R + zf * tanalpha ) * cosphi, tanalpha * cosphi
y1, y2 = yf - ye + ( R + zf * tanalpha ) * sinphi, tanalpha * sinphi
z1, z2 = ze, 1
A = ( x2 / a ) ** 2 + ( y2 / b ) ** 2 + ( z2 / c ) ** 2
B = -2 * ( x1 * x2 / a ** 2 + y1 * y2 / b ** 2 + z1 * z2 / c ** 2 )
C = ( x1 / a ) ** 2 + ( y1 / a ) ** 2 + ( z1 / c ) ** 2 - 1
discriminant = B ** 2 - 4 * A * C
if discriminant >= 0:
root1 = ( -B - math.sqrt( discriminant ) ) / ( 2 * A )
root2 = -B / A - root1
if zf <= root1 <= zf + h: lower.append( [ xf + ( R - ( root1 - zf ) * tanalpha ) * cosphi, yf + ( R - ( root1 - zf ) * tanalpha ) * sinphi, root1 ] )
if zf <= root2 <= zf + h: upper.append( [ xf + ( R - ( root2 - zf ) * tanalpha ) * cosphi, yf + ( R - ( root2 - zf ) * tanalpha ) * sinphi, root2 ] )
lower = np.array( lower )
upper = np.array( upper )
################################################################################################
# Create a 3D plot
def get_ellipsoid( c, l ):
num_points = 20
th = np.linspace( 0, 2*np.pi, num_points )
ph = np.linspace( 0, np.pi, num_points )
Theta, Phi = np.meshgrid(th, ph)
X_ell = c[0] + l[0] * np.outer(np.sin(Phi), np.cos(Theta))
Y_ell = c[1] + l[1] * np.outer(np.sin(Phi), np.sin(Theta))
Z_ell = c[2] + l[2] * np.outer(np.cos(Phi), np.ones_like(Theta))
return X_ell, Y_ell, Z_ell
def get_frustum( c, r1, alpha, h ):
num_points = 20
th = np.linspace( 0, 2*np.pi, num_points )
z = np.linspace( 0, h, num_points )
Theta, Z = np.meshgrid(th, z)
R = r1 - Z * np.tan(alpha)
X_fru = c[0] + R * np.cos(Theta)
Y_fru = c[1] + R * np.sin(Theta)
Z_fru = c[2] + Z
return X_fru, Y_fru, Z_fru
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection = '3d')
# Plot the volumes
xell, yell, zell = get_ellipsoid( [xe,ye,ze], [a,b,c] )
ax.plot_surface(xell, yell, zell, color = 'r', alpha = 0.005)
xfru, yfru, zfru = get_frustum( [xf,yf,zf], R, alpha, h )
ax.plot_surface( xfru, yfru, zfru, color = 'b', alpha = 0.35)
# Plot the intersections
if lower.size > 0: ax.plot( lower[:,0], lower[:, 1], lower[:, 2], color = 'g' )
if upper.size > 0: ax.plot( upper[:,0], upper[:, 1], upper[:, 2], color = 'g' )
# Set labels and title
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()
Sketch summary of the solution. You will have to complete the lengthy algebra to find the coefficients of the quadratic: they are given in the code.