Rotating an ellipse alongside its major axes

Question

My goal is to plot a rotational ellipsoid in 3D. I have a code for creating a simulation of the respective ellipse in 2D:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

# Define ground station and airplane positions
ground_station = np.array([2, 2])
airplane = np.array([4, 5])

# Actual distance between ground station and airplane
true_distance = np.linalg.norm(ground_station - airplane)

# Distance measured by DME
dme_distance = 1.25 * true_distance

# Calculate the center of the ellipse (midpoint between ground station and airplane)
center = (ground_station + airplane) / 2

# Calculate the semi-major and semi-minor axes of the ellipse
a = dme_distance / 2
b = np.sqrt(a**2 - (true_distance / 2)**2)

# Calculate the angle of rotation for the ellipse
angle = np.arctan2(airplane[1] - ground_station[1], airplane[0] - ground_station[0])

# Visualize the ellipse
ellipse = Ellipse(xy=center, width=2 * a, height=2 * b, angle=np.degrees(angle), edgecolor='r', linestyle='dotted', fill=False)

# Visualize simulation
plt.figure(figsize=(8, 6))
plt.plot(*ground_station, 'ro', label='Ground Station')
plt.plot(*airplane, 'bo', label='Airplane')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('2D Simulation')
plt.legend()
plt.grid(True)

plt.gca().add_patch(ellipse)
plt.axis('equal')
plt.show()

This creates such an ellpse:

I now need to extend this 2D representation into 3D. The two foci stay on the x-y-plane and the properties of the ellipse stay the same. The goal is to rotate the ellipse alongside its major axis to get a rotational ellipsoid.

My approach was to generate 360 ellipses where every ellipse is rotated for 1 degree. Somehow I cant get it to just rotate without changing its properties.

Another approach was to just plot it directly as an ellipsoid. Here's my code for this:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.patches import Ellipse

ground_station_3d = np.array([2, 2, 0])  
airplane_3d = np.array([4, 5, 0])        

true_distance = np.linalg.norm(ground_station_3d[:2] - airplane_3d[:2])

dme_distance = 1.25 * true_distance

center = (ground_station_3d + airplane_3d) / 2

a = dme_distance / 2
b = np.sqrt(a**2 - (true_distance / 2)**2)

angle = np.arctan2(airplane_3d[1] - ground_station_3d[1], airplane_3d[0] - ground_station_3d[0])

ellipse_height = 5

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = center[0] + a * np.outer(np.cos(u), np.sin(v))
y = center[1] + b * np.outer(np.sin(u), np.sin(v))
z = center[2] + ellipse_height * np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, color='r', alpha=0.5)

ax.scatter(*ground_station_3d, color='b', label='Bodenstation')
ax.scatter(*airplane_3d, color='g', label='Flugzeug')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D-Simulation')
ax.legend()

plt.show()

Since my calculations of the axes are exactly the same I expected the result to be somehow correct. But the resulting ellipsoid does not include any of the 2 focis and therefore doesn't have the same properties as the 2D one.

I need help finding the correct way of rotating the ellipse alongside its major axis. If there is no proper way to do so in python I would appreciate help on finding ways in other programming languages.

Regards

EDIT: Using the rotational matrix and plotting the ellipsoid generates the following picture:

As you can see when aligning it around the axis between both foci it seems that it isn't oriented properly and that there are points that have a longer distance on the way foci1->ellipsoid->foci2 than others.

Answer 1

Note: edited, as the OP now wants the aircraft (and hence one focus) in the air!

I suggest that your first orient your ellipsoid in the standard sense (centred at the origin with half axes a, b, c along x, y, z directions), then rotate by the appropriate angle about the z axis, together with translating to the correct centre.

For the canonical alignment:

xp = a cos(u)
yp = b sin(u) cos(v)
zp = c sin(u) sin(v)

Then you add the translation to the centre (as you are currently doing) and simultaneously rotate. The rotation consists of a rotation about the old z-axis to align the xy plane, then a rotation about the NEW y-axis to get the aircraft elevation angle. You could do this by two successive rotations, but its easier just to formulate one linear transformation by noting that its columns are just the transformed unit cartesian vectors. Thus, the new x axis (ex) stares in the direction of the aircraft, the new y axis (ey) is just rotated in the xy plane, and the new z axis (ez) is formed to make a right-handed orthogonal set. You (well, I) have to remember that numpy 1-d arrays are basically row vectors, which has consequences when creating matrices or doing matrix multiplication.

In code it looks like

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

ground_station_3d = np.array([2, 2, 0])  
airplane_3d = np.array([4, 5, 5])

true_distance = np.linalg.norm( airplane_3d - ground_station_3d )
center = ( ground_station_3d + airplane_3d ) / 2

dme_distance = 1.25 * true_distance
a = dme_distance / 2
b = np.sqrt( a ** 2 - ( true_distance / 2 ) ** 2 )
c = b

# construct rotation matrix R (columns are the transforms of the cartesian basis vectors)
diffx = airplane_3d[0] -  ground_station_3d[0]
diffy = airplane_3d[1] -  ground_station_3d[1]
ex = ( airplane_3d - ground_station_3d ) / true_distance     # new x axis stares at aircraft
ey = np.array( [ -diffy, diffx, 0 ] ) / math.sqrt( diffx * diffx + diffy * diffy )
ez = np.cross( ex, ey )
R = ( np.vstack( ( ex, ey, ez ) ) ).T

N = 100
u = np.linspace( 0, 2 * np.pi, N )
v = np.linspace( 0, np.pi, N )

# First form an ellipsoid conventionally aligned with the axes (semi-axes are a, b, c)
xp = a * np.outer( np.cos(u), np.ones(N) )
yp = b * np.outer( np.sin(u), np.cos(v) )
zp = c * np.outer( np.sin(u), np.sin(v))

# Now rotate (with R) and translate (by center)
x = np.zeros( ( N, N ) )
y = np.zeros( ( N, N ) )
z = np.zeros( ( N, N ) )
for i in range( N ):
    for j in range( N ):
        xyzp = np.array( ( xp[i,j], yp[i,j], zp[i,j] ) )   # row vector (unfortunately)
        xyz = center + xyzp @ ( R.T )                      # ditto
        x[i,j] = xyz[0]
        y[i,j] = xyz[1]
        z[i,j] = xyz[2]

# Plot
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, color='r', alpha=0.3)
ax.scatter(*ground_station_3d, color='b', label='Bodenstation')
ax.scatter(*airplane_3d, color='g', label='Flugzeug')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D-Simulation')
ax.legend()
plt.show()

Output: