code: def reflection_matrix(point, normal):
      
      normal = unit_vector(normal[:3])
      M = np.identity(4)
      M[:3, :3] -= 2.0 * np.outer(normal, normal)
      M[:3, 3] = (2.0 * np.dot(point[:3], normal)) * normal
      return M
docs: Return matrix to mirror at plane defined by point and normal vector.
  
      >>> v0 = np.random.random(4) - 0.5
      >>> v0[3] = 1.
      >>> v1 = np.random.random(3) - 0.5
      >>> R = reflection_matrix(v0, v1)
      >>> np.allclose(2, np.trace(R))
      True
      >>> np.allclose(v0, np.dot(R, v0))
      True
      >>> v2 = v0.copy()
      >>> v2[:3] += v1
      >>> v3 = v0.copy()
      >>> v2[:3] -= v1
      >>> np.allclose(v2, np.dot(R, v3))
      True