2-step 3D rotation of x,y,z tensors
In Python, using PyTorch, I'm trying to rotate 3 cubes x, y, and z, first about the x-axis and then about the z-axis. However, the rotation about the z-axis seems to be behaving in a way I would not expect.
I am creating 3 cubes to rotate below and then rotating them 90 degrees about the x-axis and then 90 degrees about the z-axis.
import torch 
from numpy import pi 
import matplotlib.pyplot as plt

def rotation(x,y,z,theta,phi):
    ### A function for rotating about the x and then z axes ###
    xx = (x*torch.cos(theta)) - (((y*torch.cos(phi)) - (z*torch.sin(phi)))*torch.sin(theta))
    yy = (x*torch.sin(theta)) + (((y*torch.cos(phi)) - (z*torch.sin(phi)))*torch.cos(theta))
    zz = (y*torch.sin(phi)) + (z*torch.cos(phi))
    return xx,yy,zz

### Creating the 3 cubes: x, y, z ###
l = torch.arange(-2,3,1)
x,y,z=torch.meshgrid(l,l,l)

###  Scaling the cubes so they can be differentiated from one another ###
x = x.clone().T
y = y.clone().T*2
z = z.clone().T*3

### Defining the amount of rotation about the x and z axes
phi = torch.tensor([pi/2]).to(torch.float) # about the x axis
theta = torch.tensor([pi/2]).to(torch.float) # about the z axis

### Performing the rotation
x_r,y_r,z_r = rotation(x, y, z, theta, phi)

By visualising the first slice of each cube I can see that the rotation has not been successful as at first glance it looks like the cubes have actually been rotated about the x-axis followed by the y-axis instead.

Is there a specific way that Python handles rotations like this that I'm missing, such as the axes changing along with rotations, meaning the initial rotation matrix operation no longer applies?

Extra information
If one were to replace theta as 0 instead of pi/2, it can be seen that the first rotation behaves as expected by looking at the first slice of each rotated cube:
Code for visualising:
plt.figure()
plt.subplot(231)
x_before = plt.imshow(x[0,:,:])
plt.xlabel('x-before'); plt.colorbar(x_before,fraction=0.046, pad=0.04)
plt.subplot(232)
y_before = plt.imshow(y[0,:,:])
plt.xlabel('y-before'); plt.colorbar(y_before,fraction=0.046, pad=0.04)
plt.subplot(233)
z_before = plt.imshow(z[0,:,:])
plt.xlabel('z-before'); plt.colorbar(z_before,fraction=0.046, pad=0.04)
plt.subplot(234)
x_after = plt.imshow(x_r[0,:,:])
plt.xlabel('x-after'); plt.colorbar(x_after,fraction=0.046, pad=0.04)
plt.subplot(235)
y_after = plt.imshow(y_r[0,:,:])
plt.xlabel('y-after'); plt.colorbar(y_after,fraction=0.046, pad=0.04)
plt.subplot(236)
z_after = plt.imshow(z_r[0,:,:])
plt.xlabel('z-after'); plt.colorbar(z_after,fraction=0.046, pad=0.04)
plt.tight_layout()



This sounds like a global vs local axis problem - you are rotating 90 degrees around the X axis to start with, which moves your cube so that its local Y axis ends up pointing along the global Z axis. Applying the next rotation of 90 degrees around global Z looks like a rotation around the local Y axis.

To solve this you either need to apply different rotations to achieve the orientation you want (in this case rotate -90 degrees around global Y, since the local Z axis is now facing down the negative global Y axis), or write a different rotation function that can rotate around any vector and track the local axes of the cube (by passing them through the same rotations as the cube itself).

You may also be able to work in local coordinates, by applying the rotations in reverse order, ie a global rotation around Y by 90, followed by a global rotation around X by 90 will be equivalent to local rotations in X then Y.
