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.

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