docs: def rotation_matrix_3D(u, th):
      """
      rotation_matrix_3D(u, t) yields a 3D numpy matrix that rotates any vector about the axis u
      t radians counter-clockwise.
      """
      # normalize the axis:
      u = normalize(u)
      # We use the Euler-Rodrigues formula;
      # see https://en.wikipedia.org/wiki/Euler-Rodrigues_formula
      a = math.cos(0.5 * th)
      s = math.sin(0.5 * th)
      (b, c, d) = -s * u
      (a2, b2, c2, d2) = (a*a, b*b, c*c, d*d)
      (bc, ad, ac, ab, bd, cd) = (b*c, a*d, a*c, a*b, b*d, c*d)
      return np.array([[a2 + b2 - c2 - d2, 2*(bc + ad),         2*(bd - ac)],
                       [2*(bc - ad),       a2 + c2 - b2 - d2,   2*(cd + ab)],
                       [2*(bd + ac),       2*(cd - ab),         a2 + d2 - b2 - c2]])