0

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.

6
  • Please include L and the desired resultant vector. Do you think you made a math error or a programmatic error? Commented Jan 15, 2018 at 15:29
  • 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. Commented Jan 15, 2018 at 15:34
  • With the L you specified, this line, line_L = [L[0,0],L[0,1]], throws a TypeError. Please provide a minimal reproducible example Commented Jan 15, 2018 at 15:51
  • 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. Commented Jan 15, 2018 at 15:55
  • 1
    Just a note that 'scipy.spatial.transform.Rotation` does all the rotations now without the need of going through all the theory. Commented Jun 11, 2019 at 19:19

2 Answers 2

8

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

Sign up to request clarification or add additional context in comments.

2 Comments

This looks like it will be a perfect fix! My only issue is that when I pass: array([ 11.23130375, 9.27144872, 18.08579023]) and array([ 0., 0., 1.]) to angle_between, I get the error: # ang1 = np.arctan2(*p1[::-1]) TypeError: return arrays must be of ArrayType
I was trying to rotate the blue vector to lie along the z axis only. I've fixed the ArrayType error now.
1

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

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.