import numpy as np

def make_rotation_matrix(phi, n):
    """
    Make the rotation matrix to rotate points through an angle `phi`
    counter-clockwise around the axis `n` (when looking towards +infinity).

    Source: Weissten, Eric W. "Rotation Formula." Mathworld.
    """
    n = np.asarray(n)/np.linalg.norm(n)
    
    return np.cos(phi)*np.identity(3) + (1-np.cos(phi))*np.outer(n,n) + \
        np.sin(phi)*np.array([[0,n[2],-n[1]],[-n[2],0,n[0]],[n[1],-n[0],0]])

def rotate(x, phi, n):
    """
    Rotate an array of points `x` through an angle phi counter-clockwise
    around the axis `n` (when looking towards +infinity).
    """
    return np.inner(np.asarray(x),make_rotation_matrix(phi, n))