docs: def _xysxy2(date):
      """Here we deviate from what has been done everywhere else. Instead of taking the formulas
      available in the Vallado, we take those described in the files tab5.2{a,b,d}.txt.
  
      The result should be equivalent, but they are the last iteration of the IAU2000A as of June 2016
  
      Args:
          date (Date)
      Return:
          3-tuple of float: Values of X, Y, s + XY/2 in arcsecond
      """
  
      planets = _planets(date)
      x_tab, y_tab, s_tab = _tab('X'), _tab('Y'), _tab('s')
  
      ttt = date.change_scale('TT').julian_century
  
      # Units: micro-arcsecond
      X = -16616.99 + 2004191742.88 * ttt - 427219.05 * ttt ** 2 - 198620.54 * ttt ** 3\
          - 46.05 * ttt ** 4 + 5.98 * ttt ** 5
  
      Y = -6950.78 - 25381.99 * ttt - 22407250.99 * ttt ** 2 + 1842.28 * ttt ** 3\
          + 1113.06 * ttt ** 4 + 0.99 * ttt ** 5
  
      s_xy2 = 94.0 + 3808.65 * ttt - 122.68 * ttt ** 2 - 72574.11 * ttt ** 3\
          + 27.98 * ttt ** 4 + 15.62 * ttt ** 5
  
      for j in range(5):
  
          _x, _y, _s = 0, 0, 0
          for i in range(len(x_tab[j])):
              Axs, Axc, *p_coefs = x_tab[j][i]
              ax_p = np.dot(p_coefs, planets)
              _x += Axs * np.sin(ax_p) + Axc * np.cos(ax_p)
  
          for i in range(len(y_tab[j])):
              Ays, Ayc, *p_coefs = y_tab[j][i]
              ay_p = np.dot(p_coefs, planets)
              _y += Ays * np.sin(ay_p) + Ayc * np.cos(ay_p)
  
          for i in range(len(s_tab[j])):
              Ass, Asc, *p_coefs = s_tab[j][i]
              as_p = np.dot(p_coefs, planets)
              _s += Ass * np.sin(as_p) + Asc * np.cos(as_p)
  
          X += _x * ttt ** j
          Y += _y * ttt ** j
          s_xy2 += _s * ttt ** j
  
      # Conversion to arcsecond
      return X * 1e-6, Y * 1e-6, s_xy2 * 1e-6