c*********************************************************************** SUBROUTINE pos_prop(t0, a0, c0, t, a, c) c*********************************************************************** c c Propagates the 6-dimensional vector of barycentric astrometric c parameters and the associated covariance matrix from epoch t0 to t c assuming uniform space motion c c This is an implementation of the equations in The Hipparcos c and Tycho Catalogues (ESA SP-1200), Volume 1, Section 1.5.5, c 'Epoch Transformation: Rigorous Treatment'. c c c INPUT (all real variables are DOUBLE PRECISION): c c t0 = original epoch [yr] - see Note 1 c a0(1) = right ascension at t0 [deg] c a0(2) = declination at t0 [deg] c a0(3) = parallax at t0 [mas] c a0(4) = proper motion in R.A., mult by cos(Dec), at t0 [mas/yr] c a0(5) = proper motion in Dec at t0 [mas/yr] c a0(6) = normalised radial velocity at t0 [mas/yr] - see Note 2 c c0(i,j) = cov[a0(i),a0(j)] in [mas^2, mas^2/yr, mas^2/yr^2] c (i, j = 1, 2, .., 6) - see Note 3 c t = new epoch [yr] - see Note 1 c c c OUTPUT (all real variables are DOUBLE PRECISION): c c a(1) = right ascension at t [deg] c a(2) = declination at t [deg] c a(3) = parallax at t [mas] c a(4) = proper motion in R.A., mult by cos(Dec), at t [mas/yr] c a(5) = proper motion in Dec at t [mas/yr] c a(6) = normalised radial velocity at t [mas/yr] - see Note 2 c c(i,j) = cov[a(i),a(j)] in [mas^2, mas^2/yr, mas^2/yr^2] c (i, j = 1, 2, .., 6) - see Note 3 c c c FUNCTIONS/SUBROUTINES CALLED: c c DOUBLE PRECISION FUNCTION scalar = scalar product of two vectors c c c NOTES: c c 1. Only t-t0 is used; the origin of the time scale is c therefore irrelevant but must be the same for t0 and t. c c 2. The normalised radial velocity at epoch t0 is given by c a0(6) = vr0*a0(3)/4.740470446 c where vr0 is the barycentric radial velocity in [km/s] at c epoch t0; similarly a(6) = vr*a(3)/4.740470446 at epoch t. c c 3. The matrices c0() and c() are dimensioned as c0(6,6), c(6,6) c in this subroutine and must have compatible dimensions in the c calling program unit. Although c0() and c() are symmetric c matrices, all elements in c0() must be defined on entry, and c all elements in c() are defined on return. c c 4. This routine adheres to ANSI Fortran 77, except for the use c of lower-case letters. All variables are declared. c c c WRITTEN BY: c c L Lindegren, Lund Observatory, 15 March 1995 c c------------------------- c Input/output variables: c------------------------- DOUBLE PRECISION t0, a0(6), c0(6,6), t, a(6), c(6,6) c------------------ c Function called: c------------------ DOUBLE PRECISION scalar c------------ c Constants: c------------ DOUBLE PRECISION zero, one, two, three PARAMETER (zero = 0d0, one = 1d0, two = 2d0, three = 3d0) DOUBLE PRECISION pi, rdpdeg, rdpmas, eps PARAMETER (pi = 3.14159265358979324d0) PARAMETER (rdpdeg = (pi/180d0)) PARAMETER (rdpmas = (pi/(180d0*3600d3))) PARAMETER (eps = 1d-9) c---------------------- c Auxiliary variables: c---------------------- INTEGER i, j, k, l DOUBLE PRECISION tau, alpha0, delta0, par0, pma0, pmd0, zeta0, 1 ca0, sa0, cd0, sd0, xy, alpha, delta, par, pma, pmd, zeta, 2 w, f, f2, f3, f4, tau2, pm02, pp0, pq0, pr0, qp0, qq0, qr0, 3 ppmz, qpmz, sum, r0(3), p0(3), q0(3), pmv0(3), d(6,6), 4 pmv(3), p(3), q(3), r(3), pmz(3) c--------------------------------------------------- c Convert input data to internal units (rad, year): c--------------------------------------------------- tau = t - t0 alpha0 = a0(1)*rdpdeg delta0 = a0(2)*rdpdeg par0 = a0(3)*rdpmas pma0 = a0(4)*rdpmas pmd0 = a0(5)*rdpmas zeta0 = a0(6)*rdpmas c------------------------------------------------ c Calculate normal triad [p0 q0 r0] at t0; r0 is c also the unit vector to the star at epoch t0: c------------------------------------------------ ca0 = dcos(alpha0) sa0 = dsin(alpha0) cd0 = dcos(delta0) sd0 = dsin(delta0) c p0(1) = -sa0 p0(2) = ca0 p0(3) = zero c q0(1) = -sd0*ca0 q0(2) = -sd0*sa0 q0(3) = cd0 c r0(1) = cd0*ca0 r0(2) = cd0*sa0 r0(3) = sd0 c----------------------- c Proper motion vector: c----------------------- pmv0(1) = p0(1)*pma0 + q0(1)*pmd0 pmv0(2) = p0(2)*pma0 + q0(2)*pmd0 pmv0(3) = p0(3)*pma0 + q0(3)*pmd0 c------------------------------- c Various auxiliary quantities: c------------------------------- tau2 = tau**2 pm02 = pma0**2 + pmd0**2 w = one + zeta0*tau f2 = one/(one + two*zeta0*tau + (pm02 + zeta0**2)*tau2) f = dsqrt(f2) f3 = f2*f f4 = f2**2 c---------------------------------------- c The position vector and parallax at t: c---------------------------------------- r(1) = (r0(1)*w + pmv0(1)*tau)*f r(2) = (r0(2)*w + pmv0(2)*tau)*f r(3) = (r0(3)*w + pmv0(3)*tau)*f c par = par0*f c--------------------------------------------------------------- c The proper motion vector and normalised radial velocity at t: c--------------------------------------------------------------- pmv(1) = (pmv0(1)*w - r0(1)*pm02*tau)*f3 pmv(2) = (pmv0(2)*w - r0(2)*pm02*tau)*f3 pmv(3) = (pmv0(3)*w - r0(3)*pm02*tau)*f3 c zeta = (zeta0 + (pm02 + zeta0**2)*tau)*f2 c------------------------------------------------ c The normal triad [p q r] at t; if r is very c close to the pole, select p towards RA=90 deg: c------------------------------------------------ xy = dsqrt(r(1)**2 + r(2)**2) IF (xy .LT. eps) THEN p(1) = zero p(2) = one p(3) = zero ELSE p(1) = -r(2)/xy p(2) = r(1)/xy p(3) = zero ENDIF c q(1) = r(2)*p(3) - r(3)*p(2) q(2) = r(3)*p(1) - r(1)*p(3) q(3) = r(1)*p(2) - r(2)*p(1) c-------------------------------------------- c Convert parameters at t to external units: c-------------------------------------------- alpha = datan2(-p(1), p(2)) IF (alpha .LT. zero) alpha = alpha + two*pi delta = datan2(r(3), xy) pma = scalar(p, pmv) pmd = scalar(q, pmv) c a(1) = alpha/rdpdeg a(2) = delta/rdpdeg a(3) = par/rdpmas a(4) = pma/rdpmas a(5) = pmd/rdpmas a(6) = zeta/rdpmas c--------------------------------------------------- c Auxiliary quantities for the partial derivatives: c--------------------------------------------------- pmz(1) = pmv0(1)*f - three*pmv(1)*w pmz(2) = pmv0(2)*f - three*pmv(2)*w pmz(3) = pmv0(3)*f - three*pmv(3)*w c pp0 = scalar(p, p0) pq0 = scalar(p, q0) pr0 = scalar(p, r0) qp0 = scalar(q, p0) qq0 = scalar(q, q0) qr0 = scalar(q, r0) ppmz = scalar(p, pmz) qpmz = scalar(q, pmz) c---------------------- c Partial derivatives: c---------------------- d(1,1) = pp0*w*f - pr0*pma0*tau*f d(1,2) = pq0*w*f - pr0*pmd0*tau*f d(1,3) = zero d(1,4) = pp0*tau*f d(1,5) = pq0*tau*f d(1,6) = -pma*tau2 c d(2,1) = qp0*w*f - qr0*pma0*tau*f d(2,2) = qq0*w*f - qr0*pmd0*tau*f d(2,3) = zero d(2,4) = qp0*tau*f d(2,5) = qq0*tau*f d(2,6) = -pmd*tau2 c d(3,1) = zero d(3,2) = zero d(3,3) = f d(3,4) = -par*pma0*tau2*f2 d(3,5) = -par*pmd0*tau2*f2 d(3,6) = -par*w*tau*f2 c d(4,1) = -pp0*pm02*tau*f3 - pr0*pma0*w*f3 d(4,2) = -pq0*pm02*tau*f3 - pr0*pmd0*w*f3 d(4,3) = zero d(4,4) = pp0*w*f3 - two*pr0*pma0*tau*f3 - three*pma*pma0*tau2*f2 d(4,5) = pq0*w*f3 - two*pr0*pmd0*tau*f3 - three*pma*pmd0*tau2*f2 d(4,6) = ppmz*tau*f2 c d(5,1) = -qp0*pm02*tau*f3 - qr0*pma0*w*f3 d(5,2) = -qq0*pm02*tau*f3 - qr0*pmd0*w*f3 d(5,3) = zero d(5,4) = qp0*w*f3 - two*qr0*pma0*tau*f3 - three*pmd*pma0*tau2*f2 d(5,5) = qq0*w*f3 - two*qr0*pmd0*tau*f3 - three*pmd*pmd0*tau2*f2 d(5,6) = qpmz*tau*f2 c d(6,1) = zero d(6,2) = zero d(6,3) = zero d(6,4) = two*pma0*w*tau*f4 d(6,5) = two*pmd0*w*tau*f4 d(6,6) = (w**2 - pm02*tau2)*f4 c------------------------------------------ c Propagate the covariance as c = d*c0*d': c------------------------------------------ DO 140 i = 1, 6 DO 130 j = 1, 6 sum = zero DO 120 k = 1, 6 DO 110 l = 1, 6 sum = sum + d(i,k)*c0(k,l)*d(j,l) 110 CONTINUE 120 CONTINUE c(i,j) = sum 130 CONTINUE 140 CONTINUE c RETURN END c*********************************************************************** DOUBLE PRECISION FUNCTION scalar(x, y) c*********************************************************************** c c Scalar product of x(1:3) and y(1:3) c DOUBLE PRECISION x(3), y(3) scalar = x(1)*y(1) + x(2)*y(2) + x(3)*y(3) RETURN END