code: def rotation_matrix(d):
      
      sin_angle = np.linalg.norm(d)
  
      if sin_angle == 0:
          return np.identity(3)
  
      d /= sin_angle
  
      eye = np.eye(3)
      ddt = np.outer(d, d)
      skew = np.array([[    0,  d[2],  -d[1]],
                       [-d[2],     0,   d[0]],
                       [ d[1], -d[0],     0]], dtype=np.float64)
  
      M = ddt + np.sqrt(1 - sin_angle**2) * (eye - ddt) + sin_angle * skew
      return M
docs: Calculates a rotation matrix given a vector d. The direction of d
      corresponds to the rotation axis. The length of d corresponds to
      the sin of the angle of rotation.
  
      Variant of: http://mail.scipy.org/pipermail/numpy-discussion/2009-March/040806.html