docs: def predict (self, mjd, complain=True, n=20000):
          """Predict the object position at a given MJD.
  
          The return value is a tuple ``(ra, dec, major, minor, pa)``, all in
          radians. These are the predicted position of the object and its
          uncertainty at *mjd*. If *complain* is True, print out warnings for
          incomplete information. *n* is the number of Monte Carlo samples to
          draw for computing the positional uncertainty.
  
          The uncertainty ellipse parameters are sigmas, not FWHM. These may be
          converted with the :data:`S2F` constant.
  
          This function relies on the external :mod:`skyfield` package.
  
          """
          import sys
          from . import ellipses
  
          self.verify (complain=complain)
  
          planets, ts = load_skyfield_data() # might download stuff from the internet
          earth = planets['earth']
          t = ts.tdb(jd = mjd + 2400000.5)
  
          # "Best" position.
  
          args = {
              'ra_hours': self.ra * R2H,
              'dec_degrees': self.dec * R2D,
          }
  
          if self.pos_epoch is not None:
              args['jd_of_position'] = self.pos_epoch + 2400000.5
  
          if self.promo_ra is not None:
              args['ra_mas_per_year'] = self.promo_ra
              args['dec_mas_per_year'] = self.promo_dec
          if self.parallax is not None:
              args['parallax_mas'] = self.parallax
          if self.vradial is not None:
              args['radial_km_per_s'] = self.vradial
  
          bestra, bestdec, _ = earth.at(t).observe(PromoEpochStar(**args)).radec()
          bestra = bestra.radians
          bestdec = bestdec.radians
  
          # Monte Carlo to get an uncertainty. As always, astronomy position
          # angle convention requires that we treat declination as X and RA as
          # Y. First, we check sanity and generate randomized parameters:
  
          if self.pos_u_maj is None and self.promo_u_maj is None and self.u_parallax is None:
              if complain:
                  print_ ('AstrometryInfo.predict(): no uncertainties '
                          'available; cannot Monte Carlo!', file=sys.stderr)
              return (bestra, bestdec, 0., 0., 0.)
  
          if self.pos_u_maj is not None:
              sd, sr, cdr = ellipses.ellbiv (self.pos_u_maj, self.pos_u_min, self.pos_u_pa)
              decs, ras = ellipses.bivrandom (self.dec, self.ra, sd, sr, cdr, n).T
          else:
              ras = np.zeros (n) + self.ra
              decs = np.zeros (n) + self.dec
  
          if self.promo_ra is None:
              pmras = np.zeros (n)
              pmdecs = np.zeros (n)
          elif self.promo_u_maj is not None:
              sd, sr, cdr = ellipses.ellbiv (self.promo_u_maj, self.promo_u_min, self.promo_u_pa)
              pmdecs, pmras = ellipses.bivrandom (self.promo_dec, self.promo_ra, sd, sr, cdr, n).T
          else:
              pmras = np.zeros (n) + self.promo_ra
              pmdecs = np.zeros (n) + self.promo_dec
  
          if self.parallax is None:
              parallaxes = np.zeros (n)
          elif self.u_parallax is not None:
              parallaxes = np.random.normal (self.parallax, self.u_parallax, n)
          else:
              parallaxes = np.zeros (n) + self.parallax
  
          if self.vradial is None:
              vradials = np.zeros (n)
          elif self.u_vradial is not None:
              vradials = np.random.normal (self.vradial, self.u_vradial, n)
          else:
              vradials = np.zeros (n) + self.vradial
  
          # Now we compute the positions and summarize as an ellipse:
  
          results = np.empty ((n, 2))
  
          for i in range (n):
              args['ra_hours'] = ras[i] * R2H
              args['dec_degrees'] = decs[i] * R2D
              args['ra_mas_per_year'] = pmras[i]
              args['dec_mas_per_year'] = pmdecs[i]
              args['parallax_mas'] = parallaxes[i]
              args['radial_km_per_s'] = vradials[i]
              ara, adec, _ = earth.at(t).observe(PromoEpochStar(**args)).radec()
              results[i] = adec.radians, ara.radians
  
          maj, min, pa = ellipses.bivell (*ellipses.databiv (results))
  
          # All done.
  
          return bestra, bestdec, maj, min, pa