0
votes

I am rotating a vector in 3D via two 2D rotations using the following code:

NOTE: L is

np.array([11.231303753070549, 9.27144871768164, 18.085790226916288])

a predefined vector shown in blue in the plot below.

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

def angle_between(p1, p2):
    ang1 = np.arctan2(*p1[::-1])
    ang2 = np.arctan2(*p2[::-1])
    return ((ang1 - ang2) % (2 * np.pi))

L = np.vstack([L,np.zeros(3)])
line_xy = [0.,1.]
line_L = [L[0,0],L[0,1]]
a = angle_between(line_xy, line_L)

def rotation(vector,theta):    
        v1_new = (vector[0]*np.cos(theta)) - (vector[1]*np.sin(theta))
        v2_new = (vector[1]*np.cos(theta)) + (vector[0]*np.sin(theta))        
        z_trans = [v1_new,v2_new,vector[2]]
        line_yz= [0.,1.]
        theta2 = angle_between(line_yz, [z_trans[1],z_trans[2]])
        v1_new = (z_trans[0]*np.cos(theta2)) - (z_trans[1]*np.sin(theta2))
        v2_new = (z_trans[1]*np.cos(theta2)) + (z_trans[0]*np.sin(theta2))
        y_trans = np.array([z_trans[0],v1_new,v2_new])        
        return z_trans,y_trans

L2,L3 = rotation(L[0,:],a)

L2 = np.vstack([L2,np.zeros(3)])
L3 = np.vstack([L3,np.zeros(3)])

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
#ax.scatter(x1*1000,y1*1000,z1*1000,c ='r',zorder=2)
ax.plot(L[:,0],L[:,1],L[:,2],color='b',zorder=1)
line = np.array([[0,0,0],[0,0,15]])
ax.plot(line[:,0],line[:,1],line[:,2],color = 'g')
ax.set_xlabel('X Kpc')
ax.set_ylabel('Y Kpc')
ax.set_zlabel('Z Kpc')

ax.plot(L2[:,0],L2[:,1],L2[:,2],color='g')
ax.plot(L3[:,0],L3[:,1],L3[:,2],color='y')

What I'm doing here is calculating the angle between x=0, y=1 (that's the line_xy part) and then rotating it around the z-axis using the first part of the rotation function:

v1_new = (vector[0]*np.cos(theta)) - (vector[1]*np.sin(theta))
v2_new = (vector[1]*np.cos(theta)) + (vector[0]*np.sin(theta))        
z_trans = [v1_new,v2_new,vector[2]]

then repeat the process but this time rotating around the x axis using the second part of the rotation function:

line_yz= [0.,1.]
theta2 = angle_between(line_yz, [z_trans[1],z_trans[2]])
v1_new = (z_trans[0]*np.cos(theta2)) - (z_trans[1]*np.sin(theta2))
v2_new = (z_trans[1]*np.cos(theta2)) + (z_trans[0]*np.sin(theta2))
y_trans = np.array([z_trans[0],v1_new,v2_new]) 

Rotations are done via the standard 2D rotation equations:

x' = x cos(theta) - y sin(theta) y' = y cos(theta) + x sin(theta)

But for some reason, after the second rotation, the line (in yellow) doesn't line up with the green line (the original target of rotating this vector).

enter image description here

I've tried checking the angles in both radians and degrees but it appears to only work with radians.

When checking the angle theta2, it comes out around 35 degrees which looks plausible.

2
Please include L and the desired resultant vector. Do you think you made a math error or a programmatic error?wwii
Hi, I have included L in the 'Note' at the top now. I think it's a programmatic error which could well be due to the way I'm approaching the problem mathematically.user8188120
With the L you specified, this line, line_L = [L[0,0],L[0,1]], throws a TypeError. Please provide a minimal reproducible examplewwii
I have updated L to be an array rather than a list so that should work now. This work is taken from a larger piece of work so apologies but I did provide to provide that.user8188120
Just a note that 'scipy.spatial.transform.Rotation` does all the rotations now without the need of going through all the theory.anishtain4

2 Answers

6
votes

I am not quite clear on your question, but hopefully this should help.

If you want to rotate a 3D vector around a particular axis, take advantage of matrix transformations instead of element wise (like you have written above). Below is code to rotate a 3-D vector around any axis:

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

def unit_vector(vector):
    """ Returns the unit vector of the vector."""
    return vector / np.linalg.norm(vector)

def angle_between(v1, v2):
    """Finds angle between two vectors"""
    v1_u = unit_vector(v1)
    v2_u = unit_vector(v2)
    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))

def x_rotation(vector,theta):
    """Rotates 3-D vector around x-axis"""
    R = np.array([[1,0,0],[0,np.cos(theta),-np.sin(theta)],[0, np.sin(theta), np.cos(theta)]])
    return np.dot(R,vector)

def y_rotation(vector,theta):
    """Rotates 3-D vector around y-axis"""
    R = np.array([[np.cos(theta),0,np.sin(theta)],[0,1,0],[-np.sin(theta), 0, np.cos(theta)]])
    return np.dot(R,vector)

def z_rotation(vector,theta):
    """Rotates 3-D vector around z-axis"""
    R = np.array([[np.cos(theta), -np.sin(theta),0],[np.sin(theta), np.cos(theta),0],[0,0,1]])
    return np.dot(R,vector)

Rotate Original Blue Vector 45 degrees (pi/2)

L_predef = np.array([11.231303753070549, 9.27144871768164, 18.085790226916288]) #blue vector
new_vect = z_rotation(L_predef, np.pi/2.0)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(np.linspace(0,L_predef[0]),np.linspace(0,L_predef[1]),np.linspace(0,L_predef[2]))
ax.plot(np.linspace(0,new_vect[0]),np.linspace(0,new_vect[1]),np.linspace(0,new_vect[2]))

plt.show()

Vector plot

0
votes

There is a general solution to this problem. Given a vector, a rotation axis and an anticlockwise angle, I wrote a simple code, which works of course also for the cases already mentioned. What it does is:

  1. projecting the vector onto the plane defined by the axis of rotation;
  2. rotating the component of the vector in the plane;
  3. finally reassembling all together to give the final result.


    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib    
    def rotve(v,erot,angle):
        rotmeasure=np.linalg.norm(erot)
        erot=erot/rotmeasure;
        norme=np.dot(v,erot)
        vplane=v-norme*erot
        plnorm=np.linalg.norm(vplane)
        ep=vplane/plnorm
        eo=np.cross(erot,ep)
        vrot=(np.cos(angle)*ep+np.sin(angle)*eo)*plnorm+norme*erot
        return(vrot)

If you want, you can check with an example which plots the "umbrella" made by the rotations:

axrot=np.array([1,0,1]); v=np.array([1.,1.,1.])
fig3 = plt.figure(3)
ax3d = fig3.add_subplot(111, projection='3d')
ax3d.quiver(0,0,0,axrot[0],axrot[1],axrot[2],length=.5, normalize=True, color='black')

angles=np.linspace(0,2,10)*np.pi
for i in range(len(angles)):
    vrot=rotve(v,axrot,angles[i]);
    ax3d.quiver(0,0,0,vrot[0],vrot[1],vrot[2],length=.1, normalize=True, color='red')
ax3d.quiver(0,0,0,v[0],v[1],v[2],length=.1, normalize=True, color='blue')
ax3d.set_title('rotations')
fig3.show()
plt.show()