沿着椭圆的主轴旋转

Question

我的目标是绘制一个3D的旋转椭球体。我有一段代码可以用来创建一个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()

这段代码生成了这样的椭圆：

现在我需要把这个2D的表示扩展到3D。两个焦点仍然在x-y平面上，椭圆的属性保持不变。我的目标是围绕椭圆的长轴旋转它，以得到一个旋转椭球体。

我的方法是生成360个椭圆，每个椭圆旋转1度。但我总是无法做到只旋转而不改变它的属性。

另一种方法是直接将其绘制为椭球体。这是我的代码：

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()

由于我对轴的计算是完全相同的，我原本期待结果是正确的。但生成的椭球体并没有包含两个焦点，因此它的属性与2D的椭圆不一样。

我需要帮助找到正确的方法来围绕椭圆的长轴旋转它。如果在Python中没有合适的方法，我也希望能找到其他编程语言的解决方案。

祝好

编辑：使用旋转矩阵绘制椭球体生成了以下图片：

正如你所看到的，当围绕两个焦点之间的轴对齐时，似乎没有正确定位，并且在焦点1到椭球体再到焦点2的路径上，有些点的距离比其他点要长。

Answer 1

注意：已编辑，因为提问者现在希望飞机在空中（因此一个焦点）！

我建议你先把椭球体按照标准方式放置（中心在原点，半轴 a、b、c 分别沿 x、y、z 方向），然后再围绕 z 轴旋转一个合适的角度，同时移动到正确的中心位置。

对于标准对齐：

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

然后你可以像现在这样添加中心的平移，同时进行旋转。旋转的过程是先围绕旧的 z 轴旋转，以对齐 xy 平面，然后再围绕新的 y 轴旋转，以获得飞机的仰角。你可以通过两次连续的旋转来实现，但更简单的方法是通过一个线性变换来表达，注意它的列就是变换后的单位笛卡尔向量。因此，新的 x 轴（ex）指向飞机的方向，新的 y 轴（ey）仅在 xy 平面内旋转，而新的 z 轴（ez）则形成一个右手正交坐标系。你（或者我）需要记住，numpy 的一维数组基本上是行向量，这在创建矩阵或进行矩阵乘法时会有影响。

在代码中看起来像这样：

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()

输出：