trouble with rotation in spherical coordinates / python implementation

Question

I am trying to implement a function for rotation in spherical coordinates. In the end I need the implementation in Javascript, but first I am trying to get it to work in python because I am more familiar with the language.

I want a function rotate which takes as parameters:

• rotation angle (alpha)
• spherical coordinates theta and phi of point to be rotated (theta, phi)
• spherical coordinates theta and phi of reference axis through the spheres origin around which the rotation should take place (THETA, PHI)

Using the approach from http://stla.github.io/stlapblog/posts/RotationSphericalCoordinates.html this is what I got so far:

Execution:

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
)

This should deliver a new coodinate theta = 90 and phi = 90. What I get is theta = 180 and phi = 90.

Here some ohter inputs and outputs/expected outputs:

The part I am really not sure about is the calculation of theta_ and psi_ in the rotate function. Within the article it says psi_ should be a 2x1 matrix, but what I get is a 2x2 matrix.

Here my implementation attempt:

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

Thanks for any help!

Answer 1

As mentioned in my comment, I don't think that using rotations around x, y, and z is the most clever solution if you actually want to rotate around an arbitrary axis. So I'd use Quaternions. This virtually uses x,y, z vectors, but the qubit solution uses all the sine, atan methods as well, so no advantage or disadvantage here

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

providing

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

and