def propagate_astrometry_and_covariance_matrix(self, a0, c0, t0, t1):
zero, one, two, three = <NUM_LIT:0>, <NUM_LIT:1>, <NUM_LIT:2>, <NUM_LIT:3><EOL>tau = t1-t0<EOL>p0, q0, r0 = normalTriad(a0[<NUM_LIT:0>], a0[<NUM_LIT:1>])<EOL>par0 = a0[<NUM_LIT:2>]*self.mastorad<EOL>pma0 = a0[<NUM_LIT:3>]*self.mastorad<EOL>pmd0 = a0[<NUM_LIT:4>]*self.mastorad<EOL>pmr0 = a0[<NUM_LIT:5>]*a0[<NUM_LIT:2>]/auKmYearPerSec*self.mastorad<EOL>pmvec0 = pma0*p0+pmd0*q0<EOL>tau2 = tau*tau<EOL>pm02 = pma0**<NUM_LIT:2> + pmd0**<NUM_LIT:2><EOL>w = one + pmr0*tau<EOL>f2 = one/(one + two*pmr0*tau + (pm02+pmr0**<NUM_LIT:2>)*tau2)<EOL>f = sqrt(f2)<EOL>f3 = f2*f<EOL>f4 = f2*f2<EOL>u = (r0*w + pmvec0*tau)*f<EOL>_, ra, dec = cartesianToSpherical(u[<NUM_LIT:0>], u[<NUM_LIT:1>], u[<NUM_LIT:2>])<EOL>par = par0*f<EOL>pmvec = (pmvec0*(one+pmr0*tau) - r0*pmr0**<NUM_LIT:2>*tau)*f3<EOL>pmr = (pmr0+(pm02 + pmr0**<NUM_LIT:2>)*tau)*f2<EOL>p, q, r = normalTriad(ra, dec)<EOL>pma = sum(p*pmvec, axis=<NUM_LIT:0>)<EOL>pmd = sum(q*pmvec, axis =<NUM_LIT:0>)<EOL>a = zeros_like(a0)<EOL>a[<NUM_LIT:0>] = ra<EOL>a[<NUM_LIT:1>] = dec<EOL>a[<NUM_LIT:2>] = par/self.mastorad<EOL>a[<NUM_LIT:3>] = pma/self.mastorad<EOL>a[<NUM_LIT:4>] = pmd/self.mastorad<EOL>a[<NUM_LIT:5>] = pmr/self.mastorad<EOL>pmz = pmvec0*f - three*pmvec*w<EOL>pp0 = sum(p*p0, axis=<NUM_LIT:0>)<EOL>pq0 = sum(p*q0, axis=<NUM_LIT:0>)<EOL>pr0 = sum(p*r0, axis=<NUM_LIT:0>)<EOL>qp0 = sum(q*p0, axis=<NUM_LIT:0>)<EOL>qq0 = sum(q*q0, axis=<NUM_LIT:0>)<EOL>qr0 = sum(q*r0, axis=<NUM_LIT:0>)<EOL>ppmz = sum(p*pmz, axis=<NUM_LIT:0>)<EOL>qpmz = sum(q*pmz, axis=<NUM_LIT:0>)<EOL>J = zeros_like(c0)<EOL>if (c0.ndim==<NUM_LIT:2>):<EOL><INDENT>J = J[newaxis,:,:]<EOL><DEDENT>J[:,<NUM_LIT:0>,<NUM_LIT:0>] = pp0*w*f - pr0*pma0*tau*f<EOL>J[:,<NUM_LIT:0>,<NUM_LIT:1>] = pq0*w*f - pr0*pmd0*tau*f<EOL>J[:,<NUM_LIT:0>,<NUM_LIT:2>] = zero<EOL>J[:,<NUM_LIT:0>,<NUM_LIT:3>] = pp0*tau*f<EOL>J[:,<NUM_LIT:0>,<NUM_LIT:4>] = pq0*tau*f<EOL>J[:,<NUM_LIT:0>,<NUM_LIT:5>] = -pma*tau2<EOL>J[:,<NUM_LIT:1>,<NUM_LIT:0>] = qp0*w*f - qr0*pma0*tau*f<EOL>J[:,<NUM_LIT:1>,<NUM_LIT:1>] = qq0*w*f - qr0*pmd0*tau*f<EOL>J[:,<NUM_LIT:1>,<NUM_LIT:2>] = zero<EOL>J[:,<NUM_LIT:1>,<NUM_LIT:3>] = qp0*tau*f<EOL>J[:,<NUM_LIT:1>,<NUM_LIT:4>] = qq0*tau*f<EOL>J[:,<NUM_LIT:1>,<NUM_LIT:5>] = -pmd*tau2<EOL>J[:,<NUM_LIT:2>,<NUM_LIT:0>] = zero<EOL>J[:,<NUM_LIT:2>,<NUM_LIT:1>] = zero<EOL>J[:,<NUM_LIT:2>,<NUM_LIT:2>] = f<EOL>J[:,<NUM_LIT:2>,<NUM_LIT:3>] = -par*pma0*tau2*f2<EOL>J[:,<NUM_LIT:2>,<NUM_LIT:4>] = -par*pmd0*tau2*f2<EOL>J[:,<NUM_LIT:2>,<NUM_LIT:5>] = -par*w*tau*f2<EOL>J[:,<NUM_LIT:3>,<NUM_LIT:0>] = -pp0*pm02*tau*f3 - pr0*pma0*w*f3<EOL>J[:,<NUM_LIT:3>,<NUM_LIT:1>] = -pq0*pm02*tau*f3 - pr0*pmd0*w*f3<EOL>J[:,<NUM_LIT:3>,<NUM_LIT:2>] = zero<EOL>J[:,<NUM_LIT:3>,<NUM_LIT:3>] = pp0*w*f3 - two*pr0*pma0*tau*f3 - three*pma*pma0*tau2*f2<EOL>J[:,<NUM_LIT:3>,<NUM_LIT:4>] = pq0*w*f3 - two*pr0*pmd0*tau*f3 - three*pma*pmd0*tau2*f2<EOL>J[:,<NUM_LIT:3>,<NUM_LIT:5>] = ppmz*tau*f2<EOL>J[:,<NUM_LIT:4>,<NUM_LIT:0>] = -qp0*pm02*tau*f3 - qr0*pma0*w*f3<EOL>J[:,<NUM_LIT:4>,<NUM_LIT:1>] = -qq0*pm02*tau*f3 - qr0*pmd0*w*f3<EOL>J[:,<NUM_LIT:4>,<NUM_LIT:2>] = zero<EOL>J[:,<NUM_LIT:4>,<NUM_LIT:3>] = qp0*w*f3 - two*qr0*pma0*tau*f3 - three*pmd*pma0*tau2*f2<EOL>J[:,<NUM_LIT:4>,<NUM_LIT:4>] = qq0*w*f3 - two*qr0*pmd0*tau*f3 - three*pmd*pmd0*tau2*f2<EOL>J[:,<NUM_LIT:4>,<NUM_LIT:5>] = qpmz*tau*f2<EOL>J[:,<NUM_LIT:5>,<NUM_LIT:0>] = zero<EOL>J[:,<NUM_LIT:5>,<NUM_LIT:1>] = zero<EOL>J[:,<NUM_LIT:5>,<NUM_LIT:2>] = zero<EOL>J[:,<NUM_LIT:5>,<NUM_LIT:3>] = two*pma0*w*tau*f4<EOL>J[:,<NUM_LIT:5>,<NUM_LIT:4>] = two*pmd0*w*tau*f4<EOL>J[:,<NUM_LIT:5>,<NUM_LIT:5>] = (w**<NUM_LIT:2> - pm02*tau2)*f4<EOL>JT = zeros_like(J)<EOL>for i in range(J.shape[<NUM_LIT:0>]):<EOL><INDENT>JT[i] = J[i].T<EOL><DEDENT>if (c0.ndim==<NUM_LIT:2>):<EOL><INDENT>c = matmul(J,matmul(c0[newaxis,:,:],JT))<EOL><DEDENT>else:<EOL><INDENT>c = matmul(J,matmul(c0,JT))<EOL><DEDENT>return a, squeeze(c)<EOL>
Propagate the covariance matrix of the astrometric parameters and radial proper motion of a
source from epoch t0 to epoch t1.

Code based on the Hipparcos Fortran implementation by Lennart Lindegren.

Parameters
----------

a0 : array_like
    6-element vector: (phi, theta, parallax, muphistar, mutheta, vrad) in units of (radians,
    radians, mas, mas/yr, mas/yr, km/s). Shape of a should be (6,) or (6,N), with N the number of
    sources for which the astrometric parameters are provided.

c0 : array_like
    Covariance matrix stored in a 6x6 element array. This can be constructed from the columns
    listed in the Gaia catalogue. The units are [mas^2, mas^2/yr, mas^2/yr^2] for the various
    elements. Note that the elements in the 6th row and column should be:
    c[6,i]=c[i,6]=c[i,3]*vrad/auKmYearPerSec for i=1,..,5 and
    c[6,6]=c[3,3]*(vrad^2+vrad_error^2)/auKmYearPerSec^2+(parallax*vrad_error/auKmYearPerSec)^2

    Shape of c0 should be (6,6) or (N,6,6).

t0 : float
    Reference epoch (Julian years).

t1 : float
    New epoch (Julian years).

Returns
-------

Astrometric parameters, including the "radial proper motion" (NOT the radial velocity), and
covariance matrix at the new epoch as a 2D matrix with the new variances on the diagional and the
covariance in the off-diagonal elements.