@@ -309,6 +309,9 @@ Bug Fixes
 
 - ``astropy.coordinates``
 
+  - Transformations between CIRS and AltAz now correctly account for the
+    location of the observer. [#5591]
+
   - GCRS frames representing a location on Earth with multiple obstimes are now
     allowed. This means that the solar system routines ``get_body``,
     ``get_moon`` and ``get_sun`` now work with non-scalar times and a
 
@@ -32,9 +32,22 @@ def cirs_to_altaz(cirs_coo, altaz_frame):
     # we use the same obstime everywhere now that we know they're the same
     obstime = cirs_coo.obstime
 
-    usrepr = cirs_coo.represent_as(UnitSphericalRepresentation)
-    cirs_ra = usrepr.lon.to(u.radian).value
-    cirs_dec = usrepr.lat.to(u.radian).value
+    # if the data are UnitSphericalRepresentation, we can skip the distance calculations
+    is_unitspherical = (isinstance(cirs_coo.data, UnitSphericalRepresentation) or
+                        cirs_coo.cartesian.x.unit == u.one)
+
+    if is_unitspherical:
+        usrepr = cirs_coo.represent_as(UnitSphericalRepresentation)
+        cirs_ra = usrepr.lon.to(u.radian).value
+        cirs_dec = usrepr.lat.to(u.radian).value
+    else:
+        # compute an "astrometric" ra/dec -i.e., the direction of the
+        # displacement vector from the observer to the target in CIRS
+        loccirs = altaz_frame.location.get_itrs(cirs_coo.obstime).transform_to(cirs_coo)
+        diffrepr = (cirs_coo.cartesian - loccirs.cartesian).represent_as(UnitSphericalRepresentation)
+
+        cirs_ra = diffrepr.lon.to(u.radian).value
+        cirs_dec = diffrepr.lat.to(u.radian).value
 
     lon, lat, height = altaz_frame.location.to_geodetic('WGS84')
     xp, yp = get_polar_motion(obstime)
@@ -52,10 +65,9 @@ def cirs_to_altaz(cirs_coo, altaz_frame):
                          altaz_frame.relative_humidity,
                          altaz_frame.obswl.value)
 
-    az, zen, ha, obs_dec, obs_ra = erfa.atioq(cirs_ra, cirs_dec, astrom)
+    az, zen, _, _, _ = erfa.atioq(cirs_ra, cirs_dec, astrom)
 
-    dat = cirs_coo.data
-    if dat.get_name() == 'unitspherical'  or dat.to_cartesian().x.unit == u.one:
+    if is_unitspherical:
         rep = UnitSphericalRepresentation(lat=u.Quantity(PIOVER2 - zen, u.radian, copy=False),
                                           lon=u.Quantity(az, u.radian, copy=False),
                                           copy=False)
@@ -95,32 +107,23 @@ def altaz_to_cirs(altaz_coo, cirs_frame):
 
     # the 'A' indicates zen/az inputs
     cirs_ra, cirs_dec = erfa.atoiq('A', az, zen, astrom)*u.radian
-
-    if isinstance(altaz_coo.data, UnitSphericalRepresentation):
-        distance = None
+    if isinstance(altaz_coo.data, UnitSphericalRepresentation) or altaz_coo.cartesian.x.unit == u.one:
+        cirs_at_aa_time = CIRS(ra=cirs_ra, dec=cirs_dec, distance=None,
+                               obstime=altaz_coo.obstime)
     else:
-        locitrs = altaz_coo.location.get_itrs(altaz_coo.obstime)
-
-        # To compute the distance in a way that is reversible with cirs_to_altaz
-        # we use basic trigonometry.  The altaz_coo's distance ("d") is one leg
-        # of a triangle, and the earth center -> EarthLocation distance ("r")
-        # is a neighboring leg.  We can also easily get the angle between the
-        # earth center->target (calculated above by apio13), and
-        # earth center->EarthLocation vectors.  This is a Side-Side-Angle
-        # situation, and the formula below is the trig formula to solve for
-        # the remaining side
-
-        ucirs = cirs_frame.realize_frame(UnitSphericalRepresentation(lon=cirs_ra, lat=cirs_dec))
-        delta = ucirs.separation(locitrs)
-        r = locitrs.spherical.distance
-        d = altaz_coo.distance
-
-        sindoverd = np.sin(delta) / d
-        distance = np.sin(delta + np.arcsin(r*sindoverd))/sindoverd
-
-    #the final transform may be a no-op if the obstimes are the same
-    return CIRS(ra=cirs_ra, dec=cirs_dec, distance=distance,
-                obstime=altaz_coo.obstime).transform_to(cirs_frame)
+        # treat the output of atoiq as an "astrometric" RA/DEC, so to get the
+        # actual RA/Dec from the observers vantage point, we have to reverse
+        # the vector operation of cirs_to_altaz (see there for more detail)
+
+        loccirs = altaz_coo.location.get_itrs(altaz_coo.obstime).transform_to(cirs_frame)
+
+        astrometric_rep = SphericalRepresentation(lon=cirs_ra, lat=cirs_dec,
+                                                  distance=altaz_coo.distance)
+        newrepr = astrometric_rep + loccirs.cartesian
+        cirs_at_aa_time = CIRS(newrepr, obstime=altaz_coo.obstime)
+
+    #this final transform may be a no-op if the obstimes are the same
+    return cirs_at_aa_time.transform_to(cirs_frame)
 
 
 @frame_transform_graph.transform(FunctionTransform, AltAz, AltAz)
 
@@ -14,7 +14,8 @@
 from ...time import Time
 from .. import (EarthLocation, get_sun, ICRS, GCRS, CIRS, ITRS, AltAz,
                 PrecessedGeocentric, CartesianRepresentation, SkyCoord,
-                SphericalRepresentation, HCRS, HeliocentricTrueEcliptic)
+                SphericalRepresentation, UnitSphericalRepresentation,
+                HCRS, HeliocentricTrueEcliptic)
 
 
 from ..._erfa import epv00
@@ -331,6 +332,10 @@ def test_gcrs_altaz_moonish(testframe):
     # now check that the distance change is similar to earth radius
     assert 1000*u.km < np.abs(moonaa.distance - moon.distance).to(u.au) < 7000*u.km
 
+    # now check that it round-trips
+    moon2 = moonaa.transform_to(moon)
+    assert_allclose(moon.cartesian.xyz, moon2.cartesian.xyz)
+
     # also should add checks that the alt/az are different for different earth locations
 
 
@@ -356,7 +361,7 @@ def test_gcrs_altaz_bothroutes(testframe):
 def test_cirs_altaz_moonish(testframe):
     """
     Sanity-check that an object resembling the moon goes to the right place with
-    a CIRS->AltAz transformation
+    a CIRS<->AltAz transformation
     """
     moon = CIRS(MOONDIST_CART, obstime=testframe.obstime)
 
@@ -368,6 +373,19 @@ def test_cirs_altaz_moonish(testframe):
     assert_allclose(moon.cartesian.xyz, moon2.cartesian.xyz)
 
 
+@pytest.mark.parametrize('testframe', totest_frames)
+def test_cirs_altaz_nodist(testframe):
+    """
+    Check that a UnitSphericalRepresentation coordinate round-trips for the
+    CIRS<->AltAz transformation.
+    """
+    coo0 = CIRS(UnitSphericalRepresentation(10*u.deg, 20*u.deg), obstime=testframe.obstime)
+
+    # check that it round-trips
+    coo1 = coo0.transform_to(testframe).transform_to(coo0)
+    assert_allclose(coo0.cartesian.xyz, coo1.cartesian.xyz)
+
+
 @pytest.mark.parametrize('testframe', totest_frames)
 def test_cirs_icrs_moonish(testframe):
     """
 
@@ -19,7 +19,7 @@
 from ...time import Time
 from ...utils import iers
 
-from ...tests.helper import pytest, assert_quantity_allclose, catch_warnings
+from ...tests.helper import pytest, assert_quantity_allclose, catch_warnings, quantity_allclose
 from .test_matching import HAS_SCIPY, OLDER_SCIPY
 
 
@@ -283,3 +283,68 @@ def test_regression_5209():
     moon = get_moon(time)
     new_coord = SkyCoord([moon])
     assert_quantity_allclose(new_coord[0].distance, moon.distance)
+
+
+def test_regression_5133():
+    N = 1000
+    np.random.seed(12345)
+    lon = np.random.uniform(-10, 10, N) * u.deg
+    lat = np.random.uniform(50, 52, N) * u.deg
+    alt = np.random.uniform(0, 10., N) * u.km
+
+    time = Time('2010-1-1')
+
+    objects = EarthLocation.from_geodetic(lon, lat, height=alt)
+    itrs_coo = objects.get_itrs(time)
+
+    homes = [EarthLocation.from_geodetic(lon=-1 * u.deg, lat=52 * u.deg, height=h)
+             for h in (0, 1000, 10000)*u.km]
+
+    altaz_frames = [AltAz(obstime=time, location=h) for h in homes]
+    altaz_coos = [itrs_coo.transform_to(f) for f in altaz_frames]
+
+    # they should all be different
+    for coo in altaz_coos[1:]:
+        assert not quantity_allclose(coo.az, coo.az[0])
+        assert not quantity_allclose(coo.alt, coo.alt[0])
+
+
+def test_itrs_vals_5133():
+    time = Time('2010-1-1')
+    el = EarthLocation.from_geodetic(lon=20*u.deg, lat=45*u.deg, height=0*u.km)
+
+    lons = [20, 30, 20]*u.deg
+    lats = [44, 45, 45]*u.deg
+    alts = [0, 0, 10]*u.km
+    coos = [EarthLocation.from_geodetic(lon, lat, height=alt).get_itrs(time)
+            for lon, lat, alt in zip(lons, lats, alts)]
+
+    aaf = AltAz(obstime=time, location=el)
+    aacs = [coo.transform_to(aaf) for coo in coos]
+
+    assert all([coo.isscalar for coo in aacs])
+
+    # the ~1 arcsec tolerance is b/c aberration makes it not exact
+    assert_quantity_allclose(aacs[0].az, 180*u.deg, atol=1*u.arcsec)
+    assert aacs[0].alt < 0*u.deg
+    assert aacs[0].distance > 50*u.km
+
+    # it should *not* actually be 90 degrees, b/c constant latitude is not
+    # straight east anywhere except the equator... but should be close-ish
+    assert_quantity_allclose(aacs[1].az, 90*u.deg, atol=5*u.deg)
+    assert aacs[1].alt < 0*u.deg
+    assert aacs[1].distance > 50*u.km
+
+    assert_quantity_allclose(aacs[2].alt, 90*u.deg, atol=1*u.arcsec)
+    assert_quantity_allclose(aacs[2].distance, 10*u.km)
+
+
+def test_regression_simple_5133():
+    t = Time('J2010')
+    obj = EarthLocation(-1*u.deg, 52*u.deg, height=[100., 0.]*u.km)
+    home = EarthLocation(-1*u.deg, 52*u.deg, height=10.*u.km)
+    aa = obj.get_itrs(t).transform_to(AltAz(obstime=t, location=home))
+
+    # az is more-or-less undefined for straight up or down
+    assert_quantity_allclose(aa.alt, [90, -90]*u.deg, rtol=1e-5)
+    assert_quantity_allclose(aa.distance, [90, 10]*u.km)