def rotate_van_dam(translation, obj, rotation_axis_dir, angle):
    #initial translation
    obj += translation

    e = rotation_axis_dir/sqrt(rotation_axis_dir*rotation_axis_dir)

    RxB = matrix([
            [1, 0, 0],
            [0, e[2], -e[1]/sqrt(e[1]^2 + e[2]^2)],
            [0, e[1], e[2]/sqrt(e[1]^2 + e[2]^2)]])

    RxBn = RxB.transpose()

    e_matrix = matrix([e[0], e[1], e[2]]).transpose()
    ee = RxB*e_matrix
    ee = vector([ee[0][0], ee[1][0], ee[2][0]])

    RyY = matrix([
            [ee[2], 0, ee[0]/sqrt(ee[0]^2 + ee[2]^2)],
            [0, 1, 0],
            [-ee[0], 0, ee[2]/sqrt(ee[0]^2 + ee[2]^2)]])

    RyYn = RyY.transpose()

    RzA = matrix([
            [cos(angle), -sin(angle), 0],
            [sin(angle), cos(angle), 0],
            [0, 0, 1]])

    M = matrix([obj[0], obj[1], obj[2]]).transpose()

    RA = RxBn*RyY*RzA*RyYn*RxB

    new_obj = RA*M
    new_obj = vector([new_obj[0][0], new_obj[1][0], new_obj[2][0]])

    #end translation
    new_obj -= translation

    return new_obj