1

I am trying to rotate two 3D objects with rotation matrices but one of them does not rotate correctly after the first rotation (rotation axis does not effect this outcome). The objects are

  • A prolate ellipsoid (3D object)
  • A mgrid containing the density values of the 3D object

I am using rotation matricec applied in the following sequence Rz(a)Ry(b)Rx(c). The ellipsoid can be rotated (together with its principal axis) in any way, but when trying to rotate the mgrid more then ónce, the objects does not follow each other. I have the following function defined to create shape and density values. (and sorry for the long post but could not think of an more quick way set up my problem). My attempt at solution can be found in the bottom of the page

Figure of current state. Left: single rotation around y. Right: rotation around all three axis

def cart2sph(x,y,z):
    return [
        np.sqrt(x**2 + y**2 + z**2),
        np.arctan2(y,x),
        np.arctan(np.sqrt(x**2+y**2)/z)
    ]
def Y20(theta):
    Y = 3*np.cos(theta)**2 - 1
    Y *= (1/4) * np.sqrt(5/np.pi)
    return Y

def Y22(theta, phi):
    Y = 2*np.cos(2*phi)*np.sin(theta)**2
    Y *= (1./4.) * np.sqrt((15.0/2.0)/np.pi)
    return Y

def Rp(theta, phi):
        RPrime = 1. + beta2*(np.cos(gamma)*Y20(theta) 
                    + 1./np.sqrt(2)*np.sin(gamma)*Y22(theta, phi))
        return RPrime 

def Density(x, y, z):
        r, phi, theta = cart2sph(x,y,z)
        density = 1 + np.exp((r-Rp(theta,phi))/a0)
        return 1./density

My rotation are made as follows, but have also tried with a specific numpy array implementation (same output)

def EulerXYZ(matrix, alpha, beta, gamma): 
    X3 = np.copy(matrix[0])
    Y3 = np.copy(matrix[1])
    Z3 = np.copy(matrix[2])
    if(len(X3.shape)<3):
        X2, Y2, Z2 = RotateX(X3, Y3, Z3, alpha)
        X1, Y1, Z1 = RotateY(X2, Y2, Z2, beta)
        X,  Y,  Z  = RotateZ(X1, Y1, Z1, gamma)
    else:
        X2, Y2, Z2 = RotateX(X3, Y3, Z3, -alpha)
        X1, Y1, Z1 = RotateY(X2, Y2, Z2, -beta)
        X,  Y,  Z  = RotateZ(X1, Y1, Z1, -gamma)
    return [X, Y, Z]

def RotateX(x1,y1,z1, gamma):
    if(gamma==0): return [x1,y1,z1]
    x = x1
    y = y1*np.cos(gamma) - z1*np.sin(gamma)
    z = y1*np.sin(gamma) + z1*np.cos(gamma)
    return [x,y,z]

def RotateY(x1,y1,z1, beta):
    if(beta==0): return [x1,y1,z1]
    x = x1*np.cos(beta) + z1*np.sin(beta)
    y = y1
    z = -x1*np.sin(beta) + z1*np.cos(beta)
    return [x,y,z]

def RotateZ(x1,y1,z1, alpha):
    if(alpha==0): return [x1,y1,z1]
    x = x1*np.cos(alpha) - y1*np.sin(alpha)
    y = x1*np.sin(alpha) + y1*np.cos(alpha)
    z = z1
    return [x,y,z]

I then initiate my object as

# Make surface
theta, phi = np.mgrid[0:np.pi:50j, 0:2*np.pi:50j]
xyz = np.array([np.sin(theta) * np.cos(phi),
                        np.sin(theta) * np.sin(phi),
                        np.cos(theta)])
Rx, Ry, Rz = Rp(theta,phi)*xyz

# Make denisty grid
Xx, Yy, Zz = np.mgrid[-2:2:50j, -2:2:50j, -2:2:50j]
rho = Density(Xx, Yy, Zz)

# Make rotation
#r1,r2,r3 = np.pi/3, 0, 0  # or 0, np.pi/3, 0 produces valid output
r1,r2,r3 = np.pi/3,np.pi/3,np.pi/3

Rx, Ry, Rz = EulerXYZ( [Rx,Ry,Rz], r1,r2,r3)
x,y,z = EulerXYZ( [Xx, Yy, Zz], r1,r2,r3 )
rho = Density(x,y,z) # make density grid in original space

From where I draw them with plotly as

fig = make_subplots(rows=1, cols=1, specs=[[{'is_3d': True}]], subplot_titles=[r'$Rotating$'])
fig.add_trace(go.Volume(x=Xx.flatten(), y=Yy.flatten(), z=Zz.flatten(), value=rho.flatten(), isomin=0.05, isomax=0.4, opacity=0.3, surface_count=4, showscale=False), 1, 1)
fig.add_trace(go.Surface(x=Rx, y=Ry, z=Rz, surfacecolor=Rx**2 + Ry**2 + Rz**2, showscale=False, colorscale='Plasma'), 1, 1)

As can bee seen from plotting the figure and looking over the data the two object does not rotate equivilant (under two rotation)

Attempt at solution

Currently I have tried

  • Generating random angles and fixed angles (problems does not arise for pi/2 rotations)
  • Applying different phases (found that the mgrid rotates in opesite direction, hence the minus sign in function decleration of EulerXYZ())
  • Different implementation of rotation matrix
  • A np.einsum solution as references in Numpy einsum() for rotation of meshgrid (but the problem only involves óne rotation, same error occours when doing two rotations.)
  • Atached a principal axis (line). The object rotates under same conditions as the surface, meaning only the mgrid object deviates from the transformation.

I am at a loss for solution and the fact that óne single rotation does not give a problem gives me a headeche. Hopefully it is just some basic fact that I am missing, and hopefulle some of your guys can help me see it.

DaKehn
  • 11
  • 3

0 Answers0