I'm setting out a rigid body in python (cuboid) with a defined height, width, depth and centre position. I want to find the normal to each face of the cuboid, and then the bound plane equations for each. I've set out some code that looks to be returning what I want in terms of the normal vector, is there a more efficient way to do this? I know I've been a bit loose with how I've used the matrices but is the general idea of this okay?

My code is as follows. I've set out a couple of rotation functions at the start, one for rotation around y and one for rotation around x. The next step will be combining these into 1, just haven't got around to it yet. I then set out a class for the rigid body, which contains functions that determine the centrepoint of each face, the vertices of each face and then the normal equation of each face. Currently it does work, when the angles of rotation are zero I get -x for front and back, -y for left and right and -z for top and bottom as the normal vectors. Is this is a good way to find the normals? I want to be using this for some monte carlo modelling in time where the angles and position are changing for each t

import numpy as np

def rotation_xy(theta, vector):
    theta = np.radians(theta)
    R = np.array([[np.cos(theta), -np.sin(theta), 0], [np.sin(theta), np.cos(theta), 0], [0, 0, 1]])
    vector_rotated = np.dot(R, vector)

    return vector_rotated

def rotation_xz(phi, vector):
    phi = np.radians(phi)
    R = np.array([[np.cos(phi), 0,np.sin(phi)], [0, 1, 0], [-np.sin(phi), 0, np.cos(phi)]])
    #vector_rotated = np.dot(R, vector)
    vector_rotated = np.dot(R, vector)

    return vector_rotated
    
class paddle():
    def __init__(self, height, width, depth, centre, theta = 0, phi = 0):
        self.height = height
        self.width = width
        self.depth = depth
        self.centre = centre
        self.theta = theta
        self.phi = phi

        self.faces = self.calculate_faces()
        
    def calculate_faces(self):
        """ 
        Calculate the plane equations for each face
        """
        half_width = self.width/2
        half_height = self.height/2
        half_depth = self.depth/2

        face_offsets = {
            "front": np.array([half_depth, 0, 0]),
            "back": np.array([-half_depth, 0, 0]),
            "left": np.array([0, -half_width, 0]),
            "right": np.array([0, half_width, 0]),
            "top": np.array([0, 0, half_height]),
            "bottom": np.array([0, 0, -half_height]),
        }
        
        plane_equations = {}
        for face_name, offset in face_offsets.items():
            # Get the four vertices for the current face dynamically
            corners = self.get_face_vertices(offset, half_width, half_height, half_depth)
            # Calculate the normal vector and D value in Ax + By + Cz + D = 0
            normal_vector, D = self.get_normal_and_D(corners, self.theta, self.phi)
            plane_equations[face_name] = (normal_vector, D)
        return plane_equations

    
    def get_face_vertices(self, offset, half_width, half_height, half_depth):
        """
        Calculate the 4 vertices of the selected face
        Uses the face offset to determine these
        If offset is non-zeo in x, then vertices are at constant x, changing y and z 
        -> Done dynamically through matrices
        for example, for the front face the offset matrix is
            [0 ,0    ,0]
            [0 ,+/-1 ,0]
            [0 ,0    ,+/-1]
        """
        face_array = np.array([half_width, half_height, half_depth])
        non_zero_indices = np.where(offset !=0)[0] # offset indices that are non-zero
        zero_indices = np.where(offset == 0)[0] # offset indices that are zero
        offset_matrix = np.eye(3)
        offset_matrix[non_zero_indices,: ] = 0
        start = self.centre 
        corners = []
        corners.append(start + np.dot(face_array, (offset_matrix * -np.eye(3))))
        corners.append(start + np.dot(face_array, offset_matrix))
        for a in zero_indices:
            offset_matrix = np.eye(3)
            offset_matrix[non_zero_indices,: ] = 0
            offset_matrix[a,a] = -1
            corners.append(start + np.dot(face_array, offset_matrix))
        return corners
    
    def get_normal_and_D(self, vertices, theta, phi):
        """
        Calculate the normal vector and D value for the plane equation.
        Uses the first 3 vertices to compute the normal.
        """
        # Get two vectors from the vertices
        v1 = vertices[1] - vertices[0]  # Vector from corner 1 to corner 2
        v2 = vertices[3] - vertices[0]  # Vector from corner 1 to corner 4
        v1, v2 = rotation_xy(theta, v1), rotation_xz(theta, v2) # rotate around z axis
        v1, v2 = rotation_xz(phi, v1), rotation_xz(phi, v2) # rotate around y axis 
        normal_vector = np.cross(v1, v2)  # Cross product to get normal vector
        normal_vector = normal_vector / np.linalg.norm(normal_vector)  # Normalize the vector

        # Calculate D for the plane equation: Ax + By + Cz + D = 0
        D = -np.dot(normal_vector, self.centre)

        return normal_vector, D

paddle = paddle(height=1, width=0.1, depth=0.04, centre=np.array([0, 0, 0]))
plane_equations = paddle.calculate_faces()
for face, (normal, D) in plane_equations.items():
    print(f"Plane equation for {face} face: {normal[0]}x + {normal[1]}y + {normal[2]}z + {D}")
    print(normal)