1

我正在尝试实现一个在球坐标中旋转的函数。最后我需要用 Javascript 实现，但首先我试图让它在 python 中工作，因为我更熟悉这种语言。

我想要一个以参数为参数的函数旋转：

  • 旋转角度（阿尔法）
  • 要旋转的点的球坐标 theta 和 phi (theta, phi)
  • 通过球体原点的参考轴的球坐标 theta 和 phi，旋转应围绕该原点进行 (THETA, PHI)

使用来自http://stla.github.io/stlapblog/posts/RotationSphericalCoordinates.html的方法，这是我到目前为止得到的：

执行：

rotate(
    alpha = 90,        # roation angle
    theta = 90,        # point theta
    phi = 0,           # point phi
    THETA = 0,         # axis theta
    PHI = 0,           # axis phi
    degrees = True     # angles in degrees
)

这应该提供一个新的坐标 theta = 90 和 phi = 90。我得到的是 theta = 180 和 phi = 90。

这里有一些其他的输入和输出/预期输出： 旋转结果（）

我真的不确定的部分是旋转函数中 theta_ 和 psi_ 的计算。在文章中它说 psi_ 应该是一个 2x1 矩阵，但我得到的是一个 2x2 矩阵。

这是我的实现尝试：

import numpy as np
from math import cos, sin, atan, pi
from cmath import exp, phase

#####################################################

def rotate(alpha, theta, phi, THETA, PHI, degrees=True):

    ## DEGREES TO RAD
    if degrees:
        alpha = pi/180 * alpha
        theta = pi/180 * theta
        phi = pi/180 * phi
        THETA = pi/180 * THETA
        PHI = pi/180 * PHI

    psi_ = Psi_(alpha, theta, phi, THETA, PHI)

    theta_  = 2 * atan(abs(psi_[1][1])/abs(psi_[0][0]))
    phi_ = phase(psi_[1][1]) - phase(psi_[0][0])

    ## RAD TO DEGREES
    if degrees:
        return theta_ * 180/pi, phi_ * 180/pi

    return theta_, phi_

#####################################################

def Psi_(alpha, theta, phi, THETA, PHI):

    return Rn(THETA, PHI, alpha) * \
           Psi(alpha, theta, phi)

#####################################################

def Psi(alpha, theta, phi):

    return np.array([
        [cos(theta)/2], 
        [exp(1j*phi) * sin(theta/2)]
    ])

#####################################################

def Rn(THETA, PHI, alpha):

    return Rz(PHI) * \
           Ry(THETA) * \
           Rz(alpha) * \
           Ry(THETA).conj().T * \
           Rz(PHI).conj().T

#####################################################

def Rx(alpha):

    return np.array([
        [cos(alpha/2), -1j * sin(alpha/2)], 
        [-1j * sin(alpha/2), cos(alpha/2)]
    ])

#####################################################

def Ry(alpha):

    return np.array([
        [cos(alpha/2), -sin(alpha/2)], 
        [sin(alpha/2), cos(alpha/2)]
    ])

#####################################################

def Rz(alpha):

    return np.array([
        [exp(-1j * alpha/2), 0], 
        [0, exp(1j * alpha/2)]
    ])

#####################################################

if __name__ == "__main__":

    print(rotate(
        alpha = 90,        # roation angle
        theta = 90,        # point theta
        phi = 0,           # point phi
        THETA = 0,         # axis theta
        PHI = 0,           # axis phi
        degrees = True     # angles in degrees
    ))

谢谢你的帮助！

4

1 回答 1

3

正如我在评论中提到的，如果您实际上想围绕任意轴旋转，我认为使用围绕x、y和的旋转不是最聪明的解决方案。z所以我会使用四元数。这实际上使用x, y,z向量，但量子位解决方案也使用了所有sine,atan方法，因此这里没有优点或缺点

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

class Quaternion( object ):
    """ 
    Simplified Quaternion class for rotation of normalized vectors only!
    """

    def __init__( self, q0, qx, qy, qz ):
        """ 
        Internally uses floats to avoid integer division issues.

        @param q0: int or float
        @param qx: int or float
        @param qy: int or float
        @param qz: int or float
        """
        self._q0 = float( q0 )
        self._qx = float( qx )
        self._qy = float( qy )
        self._qz = float( qz )
        """
        Note if interpreted as rotation q0 -> -q0 doesn't make a difference
        q0 = cos( w ) so -cos( w ) = cos( w + pi ) and as the rotation
        is by twice the angle it is either 2w or 2w + 2pi, the latter being equivalent to the former.
        """

    def conjugate(q):
        """
        @return Quaternion
        """
        conjq = Quaternion( q._q0, -q._qx, -q._qy, -q._qz )
        return conjq

    def __mul__(q, r):
        """ 
        Non commutative quaternion multiplication.
        @return Quaternion
        """
        if isinstance(r, Quaternion):
            mq0 = q._q0 * r._q0 - q._qx * r._qx - q._qy * r._qy - q._qz * r._qz
            mqx = q._q0 * r._qx + q._qx * r._q0 + q._qy * r._qz - q._qz * r._qy
            mqy = q._q0 * r._qy - q._qx * r._qz + q._qy * r._q0 + q._qz * r._qx
            mqz = q._q0 * r._qz + q._qx * r._qy - q._qy * r._qx + q._qz * r._q0
            out = Quaternion(mq0, mqx, mqy, mqz)
        else:
            raise TypeError
        return out

    def __getitem__( q, idx ):
        """
        @return float
        """
        if idx < 0:
            idx = 4 + idx
        if idx in [ 0, 1, 2, 3 ]:
            out = (q._q0, q._qx, q._qy, q._qz)[idx]
        else:
            raise IndexError
        return out

theta, phi = .4, .89
xA, yA, zA = np.sin( theta ) * np.cos( phi ), np.sin( theta ) * np.sin( phi ), np.cos( theta ) 
Theta, Phi = .5, 1.13
xB, yB, zB = np.sin( Theta ) * np.cos( Phi ), np.sin( Theta ) * np.sin( Phi ), np.cos( Theta ) 

qB = Quaternion( 0, xB, yB, zB  )

cX = [ xB ]
cY = [ yB ]
cZ = [ zB ]

for alpha in np.linspace( 0.1, 6, 20 ):
    qA = Quaternion( np.cos( 0.5 * alpha ), xA * np.sin( 0.5 * alpha ), yA * np.sin( 0.5 * alpha ), zA * np.sin( 0.5 * alpha )  )
    qAi = qA.conjugate()
    qBr = qA * ( qB * qAi )
    cX += [ qBr[1] ]
    cY += [ qBr[2] ]
    cZ += [ qBr[3] ]
    print np.arccos( qBr[3] ), np.arctan2( qBr[2], qBr[1] )


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

u = np.linspace( 0, 2 * np.pi, 50 )
v = np.linspace( 0, np.pi, 25 )
x = .9 * np.outer( np.cos( u ), np.sin( v ) )
y = .9 * np.outer( np.sin( u ), np.sin( v ) )
z = .9 * np.outer( np.ones( np.size( u ) ), np.cos( v ) )

ax.plot_wireframe( x, y, z, color='g', alpha=.3 )
ax.plot( [ 0, xA ], [ 0, yA ],[ 0, zA ], color='r' )
ax.plot( [ 0, xB ], [ 0, yB ],[ 0, zB ], color='b' )
ax.plot( cX, cY, cZ , color='b' )


plt.show()

提供

>> 0.49031916121373825 1.1522714737763464
>> 0.45533365052161895 1.2122741888530444
>> 0.41447110732929837 1.2534150991034823
>> 0.3704040237686721 1.2671812656784103
>> 0.32685242086086375 1.2421569673912964
>> 0.28897751220432055 1.1656787444306542
>> 0.26337170669521853 1.0325160977992986
>> 0.2562961184275642 0.8617797986161756
>> 0.26983294601232743 0.6990291355811976
>> 0.30014342513738007 0.5835103693125616
>> 0.3405035923275427 0.5247781593073798
>> 0.38470682535027323 0.5136174978518265
>> 0.42809208202393517 0.5372807783495164
>> 0.4673177317395864 0.5852787001209924
>> 0.49997646587457817 0.6499418738891971
>> 0.5243409810228178 0.7256665899898235
>> 0.5392333590629659 0.8081372118739611
>> 0.5439681824890205 0.8937546559885136
>> 0.5383320845959003 0.9792306451808166
>> 0.5225792805478816 1.0612632858722035

和

围绕 A 旋转 B

于 2019-05-03T08:14:27.317 回答