docs: def precesion_nutation(date):
      """Precession/nutation joint rotation matrix for the IAU2010 model
      """
  
      X, Y, s = _xys(date)
  
      d = np.arctan(np.sqrt((X**2 + Y**2) / (1 - X ** 2 - Y ** 2)))
      a = 1 / (1 + np.cos(d))
  
      return np.array([
          [1 - a * X ** 2, -a * X * Y, X],
          [-a * X * Y, 1 - a * Y ** 2, Y],
          [-X, -Y, 1 - a * (X**2 + Y**2)]
      ]) @ rot3(s)