@@ -685,6 +685,10 @@ astropy.convolution
 astropy.coordinates
 ^^^^^^^^^^^^^^^^^^^
 
+- Fix a bug where for light deflection by the Sun it was always assumed that the
+  source was at infinite distance, which in the (rare and) absolute worst-case
+  scenario could lead to errors up to 3 arcsec. [#10666]
+
 astropy.cosmology
 ^^^^^^^^^^^^^^^^^
 
 
@@ -35,10 +35,8 @@ def icrs_to_cirs(icrs_coo, cirs_frame):
 
     if icrs_coo.data.get_name() == 'unitspherical' or icrs_coo.data.to_cartesian().x.unit == u.one:
         # if no distance, just do the infinite-distance/no parallax calculation
-        usrepr = icrs_coo.represent_as(UnitSphericalRepresentation)
-        i_ra = usrepr.lon.to_value(u.radian)
-        i_dec = usrepr.lat.to_value(u.radian)
-        cirs_ra, cirs_dec = atciqz(i_ra, i_dec, astrom)
+        srepr = icrs_coo.spherical
+        cirs_ra, cirs_dec = atciqz(srepr.without_differentials(), astrom)
 
         newrep = UnitSphericalRepresentation(lat=u.Quantity(cirs_dec, u.radian, copy=False),
                                              lon=u.Quantity(cirs_ra, u.radian, copy=False),
@@ -51,11 +49,8 @@ def icrs_to_cirs(icrs_coo, cirs_frame):
         astrom_eb = CartesianRepresentation(astrom['eb'], unit=u.au,
                                             xyz_axis=-1, copy=False)
         newcart = icrs_coo.cartesian - astrom_eb
-
         srepr = newcart.represent_as(SphericalRepresentation)
-        i_ra = srepr.lon.to_value(u.radian)
-        i_dec = srepr.lat.to_value(u.radian)
-        cirs_ra, cirs_dec = atciqz(i_ra, i_dec, astrom)
+        cirs_ra, cirs_dec = atciqz(srepr.without_differentials(), astrom)
 
         newrep = SphericalRepresentation(lat=u.Quantity(cirs_dec, u.radian, copy=False),
                                          lon=u.Quantity(cirs_ra, u.radian, copy=False),
@@ -66,18 +61,11 @@ def icrs_to_cirs(icrs_coo, cirs_frame):
 
 @frame_transform_graph.transform(FunctionTransformWithFiniteDifference, CIRS, ICRS)
 def cirs_to_icrs(cirs_coo, icrs_frame):
-    srepr = cirs_coo.represent_as(SphericalRepresentation)
-    cirs_ra = srepr.lon.to_value(u.radian)
-    cirs_dec = srepr.lat.to_value(u.radian)
-
     # set up the astrometry context for ICRS<->cirs and then convert to
     # astrometric coordinate direction
-    jd1, jd2 = get_jd12(cirs_coo.obstime, 'tt')
-    x, y, s = get_cip(jd1, jd2)
-    earth_pv, earth_heliocentric = prepare_earth_position_vel(cirs_coo.obstime)
-    # erfa.apci requests TDB but TT can be used instead of TDB without any significant impact on accuracy
-    astrom = erfa.apci(jd1, jd2, earth_pv, earth_heliocentric, x, y, s)
-    i_ra, i_dec = aticq(cirs_ra, cirs_dec, astrom)
+    astrom = erfa_astrom.get().apci(cirs_coo)
+    srepr = cirs_coo.represent_as(SphericalRepresentation)
+    i_ra, i_dec = aticq(srepr.without_differentials(), astrom)
 
     if cirs_coo.data.get_name() == 'unitspherical' or cirs_coo.data.to_cartesian().x.unit == u.one:
         # if no distance, just use the coordinate direction to yield the
@@ -121,27 +109,13 @@ def cirs_to_cirs(from_coo, to_frame):
 
 @frame_transform_graph.transform(FunctionTransformWithFiniteDifference, ICRS, GCRS)
 def icrs_to_gcrs(icrs_coo, gcrs_frame):
-    # first set up the astrometry context for ICRS<->GCRS. There are a few steps...
-    # get the position and velocity arrays for the observatory.  Need to
-    # have xyz in last dimension, and pos/vel in one-but-last.
-    # (Note could use np.stack once our minimum numpy version is >=1.10.)
-    obs_pv = erfa.pav2pv(
-        gcrs_frame.obsgeoloc.get_xyz(xyz_axis=-1).to_value(u.m),
-        gcrs_frame.obsgeovel.get_xyz(xyz_axis=-1).to_value(u.m/u.s))
-
-    # find the position and velocity of earth
-    jd1, jd2 = get_jd12(gcrs_frame.obstime, 'tdb')
-    earth_pv, earth_heliocentric = prepare_earth_position_vel(gcrs_frame.obstime)
-
-    # get astrometry context object, astrom.
-    astrom = erfa.apcs(jd1, jd2, obs_pv, earth_pv, earth_heliocentric)
+    # first set up the astrometry context for ICRS<->GCRS.
+    astrom = erfa_astrom.get().apcs(gcrs_frame)
 
     if icrs_coo.data.get_name() == 'unitspherical' or icrs_coo.data.to_cartesian().x.unit == u.one:
         # if no distance, just do the infinite-distance/no parallax calculation
-        usrepr = icrs_coo.represent_as(UnitSphericalRepresentation)
-        i_ra = usrepr.lon.to_value(u.radian)
-        i_dec = usrepr.lat.to_value(u.radian)
-        gcrs_ra, gcrs_dec = atciqz(i_ra, i_dec, astrom)
+        srepr = icrs_coo.represent_as(SphericalRepresentation)
+        gcrs_ra, gcrs_dec = atciqz(srepr.without_differentials(), astrom)
 
         newrep = UnitSphericalRepresentation(lat=u.Quantity(gcrs_dec, u.radian, copy=False),
                                              lon=u.Quantity(gcrs_ra, u.radian, copy=False),
@@ -156,9 +130,7 @@ def icrs_to_gcrs(icrs_coo, gcrs_frame):
         newcart = icrs_coo.cartesian - astrom_eb
 
         srepr = newcart.represent_as(SphericalRepresentation)
-        i_ra = srepr.lon.to_value(u.radian)
-        i_dec = srepr.lat.to_value(u.radian)
-        gcrs_ra, gcrs_dec = atciqz(i_ra, i_dec, astrom)
+        gcrs_ra, gcrs_dec = atciqz(srepr.without_differentials(), astrom)
 
         newrep = SphericalRepresentation(lat=u.Quantity(gcrs_dec, u.radian, copy=False),
                                          lon=u.Quantity(gcrs_ra, u.radian, copy=False),
@@ -170,10 +142,6 @@ def icrs_to_gcrs(icrs_coo, gcrs_frame):
 @frame_transform_graph.transform(FunctionTransformWithFiniteDifference,
                                  GCRS, ICRS)
 def gcrs_to_icrs(gcrs_coo, icrs_frame):
-    srepr = gcrs_coo.represent_as(SphericalRepresentation)
-    gcrs_ra = srepr.lon.to_value(u.radian)
-    gcrs_dec = srepr.lat.to_value(u.radian)
-
     # set up the astrometry context for ICRS<->GCRS and then convert to BCRS
     # coordinate direction
     obs_pv = erfa.pav2pv(
@@ -185,7 +153,8 @@ def gcrs_to_icrs(gcrs_coo, icrs_frame):
     earth_pv, earth_heliocentric = prepare_earth_position_vel(gcrs_coo.obstime)
     astrom = erfa.apcs(jd1, jd2, obs_pv, earth_pv, earth_heliocentric)
 
-    i_ra, i_dec = aticq(gcrs_ra, gcrs_dec, astrom)
+    srepr = gcrs_coo.represent_as(SphericalRepresentation)
+    i_ra, i_dec = aticq(srepr.without_differentials(), astrom)
 
     if gcrs_coo.data.get_name() == 'unitspherical' or gcrs_coo.data.to_cartesian().x.unit == u.one:
         # if no distance, just use the coordinate direction to yield the
@@ -231,21 +200,11 @@ def gcrs_to_hcrs(gcrs_coo, hcrs_frame):
         frameattrs['obstime'] = hcrs_frame.obstime
         gcrs_coo = gcrs_coo.transform_to(GCRS(**frameattrs))
 
-    srepr = gcrs_coo.represent_as(SphericalRepresentation)
-    gcrs_ra = srepr.lon.to_value(u.radian)
-    gcrs_dec = srepr.lat.to_value(u.radian)
-
     # set up the astrometry context for ICRS<->GCRS and then convert to ICRS
     # coordinate direction
-    obs_pv = erfa.pav2pv(
-        gcrs_coo.obsgeoloc.get_xyz(xyz_axis=-1).to_value(u.m),
-        gcrs_coo.obsgeovel.get_xyz(xyz_axis=-1).to_value(u.m/u.s))
-
-    jd1, jd2 = get_jd12(hcrs_frame.obstime, 'tdb')
-    earth_pv, earth_heliocentric = prepare_earth_position_vel(gcrs_coo.obstime)
-    astrom = erfa.apcs(jd1, jd2, obs_pv, earth_pv, earth_heliocentric)
-
-    i_ra, i_dec = aticq(gcrs_ra, gcrs_dec, astrom)
+    astrom = erfa_astrom.get().apcs(gcrs_coo)
+    srepr = gcrs_coo.represent_as(SphericalRepresentation)
+    i_ra, i_dec = aticq(srepr.without_differentials(), astrom)
 
     # convert to Quantity objects
     i_ra = u.Quantity(i_ra, u.radian, copy=False)
 
@@ -14,6 +14,7 @@
 from astropy.time import Time
 from astropy.utils import iers
 from astropy.utils.exceptions import AstropyWarning
+from ..representation import CartesianRepresentation
 
 
 # We use tt as the time scale for this equinoxes, primarily because it is the
@@ -153,13 +154,16 @@ def get_cip(jd1, jd2):
     return x, y, s
 
 
-def aticq(ri, di, astrom):
+def aticq(srepr, astrom):
     """
     A slightly modified version of the ERFA function ``eraAticq``.
 
     ``eraAticq`` performs the transformations between two coordinate systems,
     with the details of the transformation being encoded into the ``astrom`` array.
 
+    There are two issues with the version of aticq in ERFA. Both are associated
+    with the handling of light deflection.
+
     The companion function ``eraAtciqz`` is meant to be its inverse. However, this
     is not true for directions close to the Solar centre, since the light deflection
     calculations are numerically unstable and therefore not reversible.
@@ -169,22 +173,30 @@ def aticq(ri, di, astrom):
     same approach used by the ERFA functions above, except that they use a threshold of
     9 arcseconds.
 
+    In addition, ERFA's aticq assumes a distant source, so there is no difference between
+    the object-Sun vector and the observer-Sun vector. This can lead to errors of up to a
+    few arcseconds in the worst case (e.g a Venus transit).
+
     Parameters
     ----------
-    ri : float or `~numpy.ndarray`
-        right ascension, radians
-    di : float or `~numpy.ndarray`
-        declination, radians
+    srepr : `~astropy.coordinates.SphericalRepresentation`
+        Astrometric GCRS or CIRS position of object from observer
     astrom : eraASTROM array
         ERFA astrometry context, as produced by, e.g. ``eraApci13`` or ``eraApcs13``
 
     Returns
     --------
     rc : float or `~numpy.ndarray`
+        Right Ascension in radians
     dc : float or `~numpy.ndarray`
+        Declination in radians
     """
+    # ignore parallax effects if no distance, or far away
+    srepr_distance = srepr.distance
+    ignore_distance = srepr_distance.unit == u.one
+
     # RA, Dec to cartesian unit vectors
-    pos = erfa.s2c(ri, di)
+    pos = erfa.s2c(srepr.lon.radian, srepr.lat.radian)
 
     # Bias-precession-nutation, giving GCRS proper direction.
     ppr = erfa.trxp(astrom['bpn'], pos)
@@ -201,7 +213,24 @@ def aticq(ri, di, astrom):
     d = np.zeros_like(pnat)
     for j in range(5):
         before = norm(pnat-d)
-        after = erfa.ld(1.0, before, before, astrom['eh'], astrom['em'], 5e-8)
+        if ignore_distance:
+            # No distance to object, assume a long way away
+            q = before
+        else:
+            # Find BCRS direction of Sun to object.
+            # astrom['eh'] and astrom['em'] contain Sun to observer unit vector,
+            # and distance, respectively.
+            eh = astrom['em'][..., np.newaxis] * astrom['eh']
+            # unit vector from Sun to object
+            q = eh + srepr_distance[..., np.newaxis].to_value(u.au) * before
+            sundist = np.linalg.norm(q, axis=-1, keepdims=True)
+            q /= sundist
+            # calculation above is extremely unstable very close to the sun
+            # in these situations, default back to ldsun-style behaviour,
+            # since this is reversible and drops to zero within stellar limb
+            q = np.where(sundist > 1.0e-10, q, before)
+
+        after = erfa.ld(1.0, before, q, astrom['eh'], astrom['em'], 1e-6)
         d = after - before
     pco = norm(pnat-d)
 
@@ -210,13 +239,16 @@ def aticq(ri, di, astrom):
     return erfa.anp(rc), dc
 
 
-def atciqz(rc, dc, astrom):
+def atciqz(srepr, astrom):
     """
     A slightly modified version of the ERFA function ``eraAtciqz``.
 
     ``eraAtciqz`` performs the transformations between two coordinate systems,
     with the details of the transformation being encoded into the ``astrom`` array.
 
+    There are two issues with the version of atciqz in ERFA. Both are associated
+    with the handling of light deflection.
+
     The companion function ``eraAticq`` is meant to be its inverse. However, this
     is not true for directions close to the Solar centre, since the light deflection
     calculations are numerically unstable and therefore not reversible.
@@ -226,25 +258,51 @@ def atciqz(rc, dc, astrom):
     same approach used by the ERFA functions above, except that they use a threshold of
     9 arcseconds.
 
+    In addition, ERFA's atciqz assumes a distant source, so there is no difference between
+    the object-Sun vector and the observer-Sun vector. This can lead to errors of up to a
+    few arcseconds in the worst case (e.g a Venus transit).
+
     Parameters
     ----------
-    rc : float or `~numpy.ndarray`
-        right ascension, radians
-    dc : float or `~numpy.ndarray`
-        declination, radians
+    srepr : `~astropy.coordinates.SphericalRepresentation`
+        Astrometric ICRS position of object from observer
     astrom : eraASTROM array
         ERFA astrometry context, as produced by, e.g. ``eraApci13`` or ``eraApcs13``
 
     Returns
     --------
     ri : float or `~numpy.ndarray`
+        Right Ascension in radians
     di : float or `~numpy.ndarray`
+        Declination in radians
     """
+    # ignore parallax effects if no distance, or far away
+    srepr_distance = srepr.distance
+    ignore_distance = srepr_distance.unit == u.one
+
     # BCRS coordinate direction (unit vector).
-    pco = erfa.s2c(rc, dc)
+    pco = erfa.s2c(srepr.lon.radian, srepr.lat.radian)
+
+    # Find BCRS direction of Sun to object
+    if ignore_distance:
+        # No distance to object, assume a long way away
+        q = pco
+    else:
+        # Find BCRS direction of Sun to object.
+        # astrom['eh'] and astrom['em'] contain Sun to observer unit vector,
+        # and distance, respectively.
+        eh = astrom['em'][..., np.newaxis] * astrom['eh']
+        # unit vector from Sun to object
+        q = eh + srepr_distance[..., np.newaxis].to_value(u.au) * pco
+        sundist = np.linalg.norm(q, axis=-1, keepdims=True)
+        q /= sundist
+        # calculation above is extremely unstable very close to the sun
+        # in these situations, default back to ldsun-style behaviour,
+        # since this is reversible and drops to zero within stellar limb
+        q = np.where(sundist > 1.0e-10, q, pco)
 
     # Light deflection by the Sun, giving BCRS natural direction.
-    pnat = erfa.ld(1.0, pco, pco, astrom['eh'], astrom['em'], 5e-8)
+    pnat = erfa.ld(1.0, pco, q, astrom['eh'], astrom['em'], 1e-6)
 
     # Aberration, giving GCRS proper direction.
     ppr = erfa.ab(pnat, astrom['v'], astrom['em'], astrom['bm1'])
 
@@ -12,6 +12,7 @@
 from astropy import _erfa as erfa
 from .utils import randomly_sample_sphere
 from astropy.coordinates.builtin_frames.utils import get_jd12, atciqz, aticq
+from astropy.coordinates import SphericalRepresentation
 
 times = [Time("2014-06-25T00:00"), Time(["2014-06-25T00:00", "2014-09-24"])]
 ra, dec, _ = randomly_sample_sphere(2)
@@ -27,7 +28,8 @@ def test_atciqz_aticq(st):
     astrom, _ = erfa.apci13(jd1, jd2)
 
     ra, dec = pos
+    srepr = SphericalRepresentation(ra, dec, 1)
     ra = ra.value
     dec = dec.value
-    assert_allclose(erfa.atciqz(ra, dec, astrom), atciqz(ra, dec, astrom))
-    assert_allclose(erfa.aticq(ra, dec, astrom), aticq(ra, dec, astrom))
+    assert_allclose(erfa.atciqz(ra, dec, astrom), atciqz(srepr, astrom))
+    assert_allclose(erfa.aticq(ra, dec, astrom), aticq(srepr, astrom))
@@ -14,14 +14,15 @@
 
 
 from astropy import units as u
-from astropy.coordinates import (AltAz, EarthLocation, SkyCoord, get_sun, ICRS,
-                GeocentricMeanEcliptic, Longitude, Latitude, GCRS, HCRS, CIRS,
-                get_moon, FK4, FK4NoETerms, BaseCoordinateFrame, ITRS,
-                QuantityAttribute, UnitSphericalRepresentation,
-                SphericalRepresentation, CartesianRepresentation,
-                FunctionTransform,
-                CylindricalRepresentation, CylindricalDifferential,
-                CartesianDifferential)
+from astropy.coordinates import (
+    AltAz, EarthLocation, SkyCoord, get_sun, ICRS,
+    GeocentricMeanEcliptic, Longitude, Latitude, GCRS, HCRS, CIRS,
+    get_moon, FK4, FK4NoETerms, BaseCoordinateFrame, ITRS,
+    QuantityAttribute, UnitSphericalRepresentation,
+    SphericalRepresentation, CartesianRepresentation,
+    FunctionTransform, get_body,
+    CylindricalRepresentation, CylindricalDifferential,
+    CartesianDifferential)
 from astropy.coordinates.sites import get_builtin_sites
 from astropy.utils.exceptions import ErfaWarning
 from astropy.time import Time
@@ -672,3 +673,20 @@ def test_regression_10422(mjd):
         loc = EarthLocation(88258.0 * u.m, -4924882.2 * u.m, 3943729.0 * u.m)
         p, v = loc.get_gcrs_posvel(obstime=t)
         assert p.shape == v.shape == t.shape
+
+
+@pytest.mark.remote_data
+def test_regression_10291():
+    """
+    According to https://eclipse.gsfc.nasa.gov/OH/transit12.html,
+    the minimum separation between Venus and the Sun during the 2012
+    transit is 554 arcseconds for an observer at the Geocenter.
+
+    If light deflection from the Sun is incorrectly applied, this increases
+    to 557 arcseconds.
+    """
+    t = Time('2012-06-06 01:29:36')
+    sun = get_body('sun', t)
+    venus = get_body('venus', t)
+    assert_quantity_allclose(venus.separation(sun),
+                             554.427*u.arcsecond, atol=0.001*u.arcsecond)