球坐标系：笛卡尔坐标系和后坐标系中的旋转

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中看到的，曲面是碎片，无法绘制。在笛卡尔空间中，旋转是成功的（如图9和图10所示，相当于3D中的x、y和z绘图）

编辑：问题似乎出在arctan2()上，正如小组2和小组4所建议的那样。我不明白为什么phi的所有值的结果theta都变为0，并且phi的值是-180的两倍。我想我应该处理arctan2()的例外情况，但我不理解逻辑。。正确的方法是什么

致以最良好的祝愿

Answer 1

简而言之：

1. plt中几乎所有的绘图实例都假定x和y数据是单调级数和/或在网格上。旋转显然破坏了它
2. 要解决此问题，可以：

• 创建网格并将原始分布投影到该网格上

  从scipy.interpolate导入网格数据 r_temp=griddata（列表（zip（phi_测试整形（-1）、ctheta_测试整形（-1））、e0mom\u LF.整形（-1），（phi.T，THETA.T），方法='linear'）

最大的缺点：对于接近零的cos（θ）值，将出现扭曲区域

• 使用b样条曲线并在常规栅格上进行渲染

  f=插值bisplrep（ctheta_检验，phi_检验，e0mom_LF，w=e0mom_LF**-0.5，s=3.5） r_new=插值.bisplev（θ，φ，f）.T

它真的很好用！球坐标上的b样条曲线有一个特定的类，但输入被限制为θ=[0，pi]phi[0,2*pi]

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