def propagate_astrometry(self, phi, theta, parallax, muphistar, mutheta, vrad, t0, t1):
        """
        Propagate the position of a source from the reference epoch t0 to the new epoch t1.

        Parameters
        ----------

        phi : float
            Longitude at reference epoch (radians).
        theta : float
            Latitude at reference epoch (radians).
        parallax : float
            Parallax at the reference epoch (mas).
        muphistar : float
            Proper motion in longitude (including cos(latitude) term) at reference epoch (mas/yr).
        mutheta : float
            Proper motion in latitude at reference epoch (mas/yr).
        vrad : float
            Radial velocity at reference epoch (km/s).
        t0 : float
            Reference epoch (Julian years).
        t1 : float
            New epoch (Julian years).

        Returns
        -------

        Astrometric parameters, including the "radial proper motion" (NOT the radial velocity), at the new epoch.
        phi1, theta1, parallax1, muphistar1, mutheta1, mur1 = epoch_prop_pos(..., t0, t1)
        """

        t = t1-t0
        p0, q0, r0 = normalTriad(phi, theta)

        # Convert input data to units of radians and Julian year. Use ICRS coordinate names internally to
        # avoid errors in translating the formulae to code.
        pmra0 = muphistar*self.mastorad
        pmdec0 = mutheta*self.mastorad
        plx0 = parallax*self.mastorad
        pmr0 = vrad*parallax/auKmYearPerSec*self.mastorad
        pmtot0sqr = (muphistar**2 + mutheta**2) * self.mastorad**2

        # Proper motion vector
        pmvec0 = pmra0*p0+pmdec0*q0

        f = (1 + 2*pmr0*t + (pmtot0sqr+pmr0**2)*t**2)**(-0.5)
        u = (r0*(1+pmr0*t) + pmvec0*t)*f

        _, phi1, theta1 = cartesianToSpherical(u[0], u[1], u[2])
        parallax1 = parallax*f
        pmr1 = (pmr0+(pmtot0sqr + pmr0**2)*t)*f**2
        pmvec1 = (pmvec0*(1+pmr0*t) - r0*pmr0**2*t)*f**3
        p1, q1, r1 = normalTriad(phi1, theta1)
        muphistar1 = sum(p1*pmvec1/self.mastorad, axis=0)
        mutheta1 = sum(q1*pmvec1/self.mastorad, axis =0)
        murad1 = pmr1/self.mastorad

        return phi1, theta1, parallax1, muphistar1, mutheta1, murad1
Propagate the position of a source from the reference epoch t0 to the new epoch t1.

        Parameters
        ----------

        phi : float
            Longitude at reference epoch (radians).
        theta : float
            Latitude at reference epoch (radians).
        parallax : float
            Parallax at the reference epoch (mas).
        muphistar : float
            Proper motion in longitude (including cos(latitude) term) at reference epoch (mas/yr).
        mutheta : float
            Proper motion in latitude at reference epoch (mas/yr).
        vrad : float
            Radial velocity at reference epoch (km/s).
        t0 : float
            Reference epoch (Julian years).
        t1 : float
            New epoch (Julian years).

        Returns
        -------

        Astrometric parameters, including the "radial proper motion" (NOT the radial velocity), at the new epoch.
        phi1, theta1, parallax1, muphistar1, mutheta1, mur1 = epoch_prop_pos(..., t0, t1)