docs: def R_rot_3d(th):
      """Return a 3-dimensional rotation matrix.
  
      Parameters
      ----------
      th: array, shape (n, 3)
          Angles about which to rotate along each axis.
  
      Returns
      -------
      R: array, shape (n, 3, 3)
      """
      sx, sy, sz = np.sin(th).T
      cx, cy, cz = np.cos(th).T
      R = np.empty((len(th), 3, 3), dtype=np.float)
  
      R[:, 0, 0] = cy * cz
      R[:, 0, 1] = -cy * sz
      R[:, 0, 2] = sy
  
      R[:, 1, 0] = sx * sy * cz + cx * sz
      R[:, 1, 1] = -sx * sy * sz + cx * cz
      R[:, 1, 2] = -sx * cy
  
      R[:, 2, 0] = -cx * sy * cz + sx * sz
      R[:, 2, 1] = cx * sy * sz + sx * cz
      R[:, 2, 2] = cx * cy
      return R