This is what I have now. It would be more difficulty to hunt down bugs stemming from this method, so it's time to run it through the python interpreter and fix things. X = np.array([1.0,0,0,0]) Y = np.array([0,1.0,0,0]) Z = np.array([0,0,1.0,0]) W = np.array([0,0,0,1.0]) class CoordFrame: def __init__(self, mat = None): if mat is None: mat = np.identity(4) self.mat = mat def apply(self, vec): return self.mat @ vec def inverted(self): return CoordFrame(self.mat.inverse()) @classmethod def fromaxisangle(cls, axis, angle): # from https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_... cos_theta = np.cos(angle) sin_theta = np.sin(angle) axis = axis / np.linalg.norm(axis) return CoordFrame( cos_theta * identity(4) + sin_theta * np.cross(axis, -np.identity(4)) + (1 - cos_theta) * np.outer(axis) )