def METHOD_NAME(axis, r, p, t):
    """
    Function that rotates the point about the axis with given angle
    # 1) Translate space so that the reference point locate at the origin (T)
    # 2) Rotate space so that the rotation axis lies in the xz plane (R_x)
    # 3) Rotate space so that the rotation axis lies in the z axis (R_y)
    # 4) Rotate angle t about the z axis (R_z)
    # 5) Rotate and translate space back to original space (T_inv, R_xinv, R_yinv)
    Parameters
    ----------
    axis : The rotation axis vector
    r : The reference point of the axis vector
    p : The point to rotate
    t : Rotation angle
    Returns
    -------
    new : new coordinates
    """
    unit_axis = axis / np.linalg.norm(axis) # normalize axis vector
    a = unit_axis[0]
    b = unit_axis[1]
    c = unit_axis[2]
    d = np.sqrt(unit_axis[1]**2 + unit_axis[2]**2)
    T = [[1, 0, 0, -r[0]], [0, 1, 0, -r[1]], [0, 0, 1, -r[2]], [0, 0, 0, 1]]
    T_inv = np.linalg.inv(T)
    old = [p[0], p[1], p[2], 1]
    if d != 0:
        R_x = [[1, 0, 0, 0], [0, c/d, -b/d, 0], [0, b/d, c/d, 0], [0, 0, 0, 1]]
        R_xinv = np.linalg.inv(R_x)
        R_y = [[d, 0, -a, 0], [0, 1, 0, 0], [a, 0, d, 0], [0, 0, 0, 1]]
        R_yinv = np.linalg.inv(R_y)
        R_z = [[np.cos(t), np.sin(t), 0, 0], [-np.sin(t), np.cos(t), 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
        # new = T_inv * R_xinv * R_yinv * R_z * R_y * R_x * T * old
        new = T_inv.dot(R_xinv).dot(R_yinv).dot(R_z).dot(R_y).dot(R_x).dot(T).dot(old)
    else: # if d ==0, rotation axis is along the x axis -> no rotation along y/z axis
        R_x = [[1, 0, 0, 0], [0, np.cos(t), np.sin(t), 0], [0, -np.sin(t), np.cos(t), 0], [0, 0, 0, 1]]
        # new = T_inv * R_x * T * old
        new = T_inv.dot(R_x).dot(T).dot(old)
    return (new[0:3])