def propagate_astrometry(self, phi, theta, parallax, muphistar, mutheta, vrad, t0, t1):
t = t1-t0<EOL>p0, q0, r0 = normalTriad(phi, theta)<EOL>pmra0 = muphistar*self.mastorad<EOL>pmdec0 = mutheta*self.mastorad<EOL>plx0 = parallax*self.mastorad<EOL>pmr0 = vrad*parallax/auKmYearPerSec*self.mastorad<EOL>pmtot0sqr = (muphistar**<NUM_LIT:2> + mutheta**<NUM_LIT:2>) * self.mastorad**<NUM_LIT:2><EOL>pmvec0 = pmra0*p0+pmdec0*q0<EOL>f = (<NUM_LIT:1> + <NUM_LIT:2>*pmr0*t + (pmtot0sqr+pmr0**<NUM_LIT:2>)*t**<NUM_LIT:2>)**(-<NUM_LIT:0.5>)<EOL>u = (r0*(<NUM_LIT:1>+pmr0*t) + pmvec0*t)*f<EOL>_, phi1, theta1 = cartesianToSpherical(u[<NUM_LIT:0>], u[<NUM_LIT:1>], u[<NUM_LIT:2>])<EOL>parallax1 = parallax*f<EOL>pmr1 = (pmr0+(pmtot0sqr + pmr0**<NUM_LIT:2>)*t)*f**<NUM_LIT:2><EOL>pmvec1 = (pmvec0*(<NUM_LIT:1>+pmr0*t) - r0*pmr0**<NUM_LIT:2>*t)*f**<NUM_LIT:3><EOL>p1, q1, r1 = normalTriad(phi1, theta1)<EOL>muphistar1 = sum(p1*pmvec1/self.mastorad, axis=<NUM_LIT:0>)<EOL>mutheta1 = sum(q1*pmvec1/self.mastorad, axis =<NUM_LIT:0>)<EOL>murad1 = pmr1/self.mastorad<EOL>return phi1, theta1, parallax1, muphistar1, mutheta1, murad1<EOL>
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)