球面坐标：笛卡尔和反向旋转

Question

我试图在笛卡尔坐标旋转后将分布重新转换为球坐标。形状是一个椭圆体，其公式取自 this paper。这是一个小的工作示例：

import math
import numpy as np
import pandas as pd

import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource
%matplotlib inline

import itertools

def rot3d(alpha,beta,gamma):
    """
    It computes the intrinsic rotation according to the convention z,x',y''.
    The angles are α around -z, β around y, and γ around x. 
    It returns 3X3 Rtot=Rz(α)Ry(β)Rx(γ) matrix
    ref: https://de.wikipedia.org/wiki/Eulersche_Winkel
    """
    rot_matrix = np.array([[np.cos(beta)*np.cos(alpha), np.sin(gamma)*np.sin(beta)*np.cos(alpha)-np.cos(gamma)*np.sin(alpha), np.cos(gamma)*np.sin(beta)*np.cos(alpha)+np.sin(gamma)*np.sin(alpha)], 
                           [np.cos(beta)*np.sin(alpha), np.sin(gamma)*np.sin(beta)*np.sin(alpha)+np.cos(gamma)*np.cos(alpha), np.cos(gamma)*np.sin(beta)*np.sin(alpha)-np.sin(gamma)*np.cos(alpha)],
                           [-np.sin(beta)             , np.sin(gamma)*np.cos(beta)                                          , np.cos(gamma)*np.cos(beta)                                         ]])
    return rot_matrix

nth=5;nph=10;
theta = np.linspace(0, np.pi, nth)
phi = np.linspace(-np.pi, np.pi, nph)

THETA, PHI = np.meshgrid(theta, phi)

THETA_deg=[el*180/np.pi for el in THETA]
PHI_deg=[el*180/np.pi for el in PHI]

#Definiton of the ellipsoid
# from https://arxiv.org/pdf/1104.5145.pdf
a=1; b=2; c=3
R = (a*b*c) / np.sqrt(b**2*c**2*np.cos(THETA)**2 + c**2*a**2*np.sin(THETA)**2*np.cos(PHI)**2 + a**2*b**2*np.sin(THETA)**2*np.sin(PHI)**2)

#Conversion to cartesian coordiantes
e0mom_x_MF = R.reshape(-1) * np.sin(THETA.reshape(-1)) * np.cos(PHI.reshape(-1))
e0mom_y_MF = R.reshape(-1) * np.sin(THETA.reshape(-1)) * np.sin(PHI.reshape(-1))
e0mom_z_MF = R.reshape(-1) * np.cos(THETA.reshape(-1))

xyzm = np.stack((e0mom_x_MF, e0mom_y_MF, e0mom_z_MF))

#Ratation in cartesia
#arbitrary rotation of (30°,15°,0.)
e0mom_x_LF, e0mom_y_LF, e0mom_z_LF = np.einsum('ik, kj -> ij', rot3d(10.*np.pi/180.,0.*np.pi/180.,0.*np.pi/180.), xyzm)
e0mom_LF = np.sqrt(e0mom_x_LF**2+e0mom_y_LF**2+e0mom_z_LF**2)

#Riconversion in spherical coordiantes
r_test=(e0mom_LF)
ctheta_test=np.arccos(e0mom_z_LF/e0mom_LF)*180/np.pi
phi_test=np.arctan2(e0mom_y_LF,e0mom_x_LF)*180./np.pi

#Plotting
fig1, ax = plt.subplots(5,2, figsize=(10,25), constrained_layout=True)

ax[0,0].scatter(PHI_deg, THETA_deg)
ax[0,0].set_xlabel('phi [DEG]')
ax[0,0].set_ylabel('theta [DEG]')
ax[0,0].set_title("initial coordinates")

ax[0,1].scatter(phi_test,ctheta_test)
ax[0,1].set_xlabel('phi [DEG]')
ax[0,1].set_ylabel('theta [DEG]')
ax[0,1].set_title("rotated coordiantes")

ax[1,0].plot(ctheta_test, "o")
ax[1,0].set_xlabel('# of bin')
ax[1,0].set_ylabel('theta [DEG]')
ax[1,0].set_title("rotated distrib")

ax[1,1].plot(phi_test, "o")
ax[1,1].set_xlabel('# of bin')
ax[1,1].set_ylabel('phi [DEG]')
ax[1,1].set_title("rotated distrib")

ax[2,0].plot(np.array(THETA_deg).reshape(-1), "o")
ax[2,0].set_xlabel('# of bin')
ax[2,0].set_ylabel('theta [DEG]')
ax[2,0].set_title("initial distrib")

ax[2,1].plot(np.array(PHI_deg).reshape(-1), "o")
ax[2,1].set_xlabel('# of bin')
ax[2,1].set_ylabel('phi [DEG]')
ax[2,1].set_title("initial distrib")

ax[3,0].contourf(PHI_deg, THETA_deg, np.array(R).reshape(nph,nth))
ax[3,0].set_xlabel('phi [DEG]')
ax[3,0].set_ylabel('theta [DEG]')
ax[3,0].set_title("Original ellipsoind in spherical coordinates")

ax[3,1].contourf(phi_test.reshape(nth,nph), ctheta_test.reshape(nth,nph), r_test.reshape(nth,nph))
ax[3,1].set_xlabel('phi [DEG]')
ax[3,1].set_ylabel('theta [DEG]')
ax[3,1].set_title("Rotated ellipsoind in spherical coordinates")

ax[4,0].contourf(np.array(R).reshape(nth,nph))
ax[4,0].set_xlabel('shape[0]')
ax[4,0].set_ylabel('shape[1]')
ax[4,0].set_title("original R array (no theta/phi)")

ax[4,1].contourf(r_test.reshape(nth,nph))
ax[4,1].set_xlabel('shape[0]')
ax[4,1].set_ylabel('shape[1]')
ax[4,1].set_title("rotated R array (no theta/phi)")

正如您在附加输出的面板 8 中看到的那样，表面支离破碎，无法绘制。在笛卡尔 space 中，旋转成功（如面板 9 和 10 所示，相当于绘制 3D 中的 x、y 和 z）。

编辑：问题似乎出在 arctan2() 面板 2 和面板 4 上。我不明白为什么对于 theta 的所有值，结果 phi 都变为 0，并且 phi 的 -180 值有两倍。我想我应该处理 arctan2() 的例外情况，但我无法理解其中的逻辑。应该遵循的正确方法是什么？

此致，

Answer 1

简要说明：

plt 中的几乎所有绘图实例都假定 x 和 y 数据在网格上呈单调递增 and/or。旋转显然打破了它。
要克服这个问题，可以：

创建一个网格并将原始分布投影到它上面

从 scipy.interpolate 导入网格数据 r_temp = griddata(列表(zip(phi_test.reshape(-1), ctheta_test.reshape(-1))),e0mom_LF.reshape(-1), (PHI.T ,THETA.T),方法='linear')

大缺点：cos(theta) 的值接近于零时会出现扭曲区域。

使用 b 样条并在规则网格上渲染它

f = interpolate.bisplrep(ctheta_test, phi_test, e0mom_LF, w=e0mom_LF**-0.5, s=3.5) r_new = interpolate.bisplev(THETA,PHI,f).T

效果真的很好！球坐标上的b样条有一个特定的class，但输入被限制为theta=[0,pi] phi[0,2*pi].

最好使用 plt 中的 pcolormesh 来可视化数据。

球面坐标：笛卡尔和反向旋转

Spherical coordiantes: rotation in cartesian and back

python

rotation