NumPy mgrid二次旋转异常问题求助（Python）

我尝试按**R_z(a)R_y(b)R_x(c)**的顺序旋转两个3D对象，但多次旋转后两者结果无法匹配：

• 长椭球（含主轴）：任意旋转均正常
• 包含该椭球密度值的NumPy mgrid：单次旋转正常，多次旋转后结果偏离

形状与密度生成函数 #

import numpy as np

def cart2sph(x,y,z):
    return [
        np.sqrt(x**2 + y**2 + z**2),
        np.arctan2(y,x),
        np.arctan(np.sqrt(x**2+y**2)/z)
    ]
def Y20(theta):
    Y = 3*np.cos(theta)**2 - 1
    Y *= (1/4) * np.sqrt(5/np.pi)
    return Y

def Y22(theta, phi):
    Y = 2*np.cos(2*phi)*np.sin(theta)**2
    Y *= (1./4.) * np.sqrt((15.0/2.0)/np.pi)
    return Y

def Rp(theta, phi):
        RPrime = 1. + beta2*(np.cos(gamma)*Y20(theta) 
                    + 1./np.sqrt(2)*np.sin(gamma)*Y22(theta, phi))
        return RPrime 

def Density(x, y, z):
        r, phi, theta = cart2sph(x,y,z)
        density = 1 + np.exp((r-Rp(theta,phi))/a0)
        return 1./density

旋转实现代码 #

def EulerXYZ(matrix, alpha, beta, gamma): 
    X3 = np.copy(matrix[0])
    Y3 = np.copy(matrix[1])
    Z3 = np.copy(matrix[2])
    if(len(X3.shape)<3):
        X2, Y2, Z2 = RotateX(X3, Y3, Z3, alpha)
        X1, Y1, Z1 = RotateY(X2, Y2, Z2, beta)
        X,  Y,  Z  = RotateZ(X1, Y1, Z1, gamma)
    else:
        X2, Y2, Z2 = RotateX(X3, Y3, Z3, -alpha)
        X1, Y1, Z1 = RotateY(X2, Y2, Z2, -beta)
        X,  Y,  Z  = RotateZ(X1, Y1, Z1, -gamma)
    return [X, Y, Z]

def RotateX(x1,y1,z1, gamma):
    if(gamma==0): return [x1,y1,z1]
    x = x1
    y = y1*np.cos(gamma) - z1*np.sin(gamma)
    z = y1*np.sin(gamma) + z1*np.cos(gamma)
    return [x,y,z]

def RotateY(x1,y1,z1, beta):
    if(beta==0): return [x1,y1,z1]
    x = x1*np.cos(beta) + z1*np.sin(beta)
    y = y1
    z = -x1*np.sin(beta) + z1*np.cos(beta)
    return [x,y,z]

def RotateZ(x1,y1,z1, alpha):
    if(alpha==0): return [x1,y1,z1]
    x = x1*np.cos(alpha) - y1*np.sin(alpha)
    y = x1*np.sin(alpha) + y1*np.cos(alpha)
    z = z1
    return [x,y,z]

初始化与绘图代码 #

from plotly.subplots import make_subplots
import plotly.graph_objects as go

# Make surface
theta, phi = np.mgrid[0:np.pi:50j, 0:2*np.pi:50j]
xyz = np.array([np.sin(theta) * np.cos(phi),
                        np.sin(theta) * np.sin(phi),
                        np.cos(theta)])
Rx, Ry, Rz = Rp(theta,phi)*xyz

# Make denisty grid
Xx, Yy, Zz = np.mgrid[-2:2:50j, -2:2:50j, -2:2:50j]
rho = Density(Xx, Yy, Zz)

# Make rotation
#r1,r2,r3 = np.pi/3, 0, 0  # or 0, np.pi/3, 0 produces valid output
r1,r2,r3 = np.pi/3,np.pi/3,np.pi/3

Rx, Ry, Rz = EulerXYZ( [Rx,Ry,Rz], r1,r2,r3)
x,y,z = EulerXYZ( [Xx, Yy, Zz], r1,r2,r3 )
rho = Density(x,y,z) # make density grid in original space

# Plot
fig = make_subplots(rows=1, cols=1, specs=[[{'is_3d': True}]], subplot_titles=[r'$Rotating$'])
fig.add_trace(go.Volume(x=Xx.flatten(), y=Yy.flatten(), z=Zz.flatten(), value=rho.flatten(), isomin=0.05, isomax=0.4, opacity=0.3, surface_count=4, showscale=False), 1, 1)
fig.add_trace(go.Surface(x=Rx, y=Ry, z=Rz, surfacecolor=Rx**2 + Ry**2 + Rz**2, showscale=False, colorscale='Plasma'), 1, 1)
fig.show()

已尝试的解决方案 #

• 测试随机角度与固定角度（π/2旋转无问题）
• 调整旋转相位（发现mgrid旋转方向相反，已在EulerXYZ()中添加负号）
• 尝试不同的旋转矩阵实现方式
• 参考使用np.einsum的单次旋转方案（多次旋转仍出错）
• 添加主轴直线测试，确认主轴与椭球旋转一致，仅mgrid变换结果偏离

内容的提问来源于stack exchange，提问作者DaKehn