Repository: astropy/astropy

Issue:
Erfa function apco13 is unsupported?
### Description

In https://github.com/astropy/astropy/pull/17729#discussion_r1946538474, it was shown that if one passes in `Quantity` to `erfa.apco13`, it raises
```
TypeError: unknown ufunc apco13. If you believe this ufunc should be
supported, please raise an i...
```

Indeed, that function is not listed in `astropy/units/quantity_helper/erfa.py`, so it should be added.


Discussion:
github-actions[bot]: Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

  - [ ] Do the proposed changes actually accomplish desired goals?
  - [ ] Do the proposed changes follow the [Astropy coding guidelines](https://docs.astropy.org/en/latest/development/codeguide.html)?
  - [ ] Are tests added/updated as required? If so, do they follow the [Astropy testing guidelines](https://docs.astropy.org/en/latest/development/testguide.html)?
  - [ ] Are docs added/updated as required? If so, do they follow the [Astropy documentation guidelines](https://docs.astropy.org/en/latest/development/docguide.html)?
  - [ ] Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for [rebase](https://docs.astropy.org/en/latest/development/development_details.html#rebase-if-necessary) and [squash](https://docs.astropy.org/en/latest/development/development_details.html#squash-if-necessary).
  - [ ] Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the [bot](https://docs.astropy.org/en/latest/development/development_details.html#pre-commit).
  - [ ] Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
  - [ ] Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
  - [ ] At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) *before* merge.
mhvk: Note that the `s390x` and `ppc64le` failures are unrelated to this PR, while the -dev problems are due to timeouts. So, I think this should be ready...

Patch:
diff --git a/astropy/units/quantity_helper/erfa.py b/astropy/units/quantity_helper/erfa.py
index e88307d6da..75a0ec2210 100644
--- a/astropy/units/quantity_helper/erfa.py
+++ b/astropy/units/quantity_helper/erfa.py
@@ -22,8 +22,8 @@ erfa_ufuncs = (
     "s2c", "s2p", "c2s", "p2s", "pm", "pdp", "pxp", "rxp", "cpv", "p2pv", "pv2p",
     "pv2s", "pvdpv", "pvm", "pvmpv", "pvppv", "pvstar", "pvtob", "pvu", "pvup",
     "pvxpv", "rxpv", "s2pv", "s2xpv", "starpv", "sxpv", "trxpv", "gd2gc", "gd2gce",
-    "gc2gd", "gc2gde", "ldn", "aper", "apio", "atciq", "atciqn", "atciqz", "aticq",
-    "atioq", "atoiq",
+    "gc2gd", "gc2gde", "ldn", "apco13", "aper", "apio",
+    "atciq", "atciqn", "atciqz", "aticq", "atioq", "atoiq",
 )  # fmt: skip
 
 
@@ -50,6 +50,11 @@ def check_structured_unit(unit, dtype):
         raise UnitTypeError(f"{msg} input needs unit matching dtype={dtype}.")
 
 
+def check_no_time_units(unit1, unit2, f):
+    if unit1 is not None or unit2 is not None:
+        raise TypeError(f"cannot pass in units for 2-part time in {f.__name__}.")
+
+
 def helper_s2c(f, unit1, unit2):
     from astropy.units.si import radian
 
@@ -319,6 +324,52 @@ def helper_aper(f, unit_theta, unit_astrom):
     return [get_converter(unit_theta, unit_along), None], result_unit
 
 
+def helper_apco13(
+    f,
+    unit_utc1,
+    unit_utc2,
+    unit_dut1,
+    unit_elong,
+    unit_phi,
+    unit_hm,
+    unit_xp,
+    unit_yp,
+    unit_phpa,
+    unit_tc,
+    unit_rh,
+    unit_wl,
+):
+    from astropy.units import (
+        add_enabled_equivalencies,
+        deg_C,
+        hPa,
+        m,
+        micron,
+        one,
+        radian,
+        second,
+        temperature,
+    )
+
+    check_no_time_units(unit_utc1, unit_utc2, f)
+
+    with add_enabled_equivalencies(temperature()):
+        return [
+            None,
+            None,
+            get_converter(unit_dut1, second),
+            get_converter(unit_elong, radian),
+            get_converter(unit_phi, radian),
+            get_converter(unit_hm, m),
+            get_converter(unit_xp, radian),
+            get_converter(unit_yp, radian),
+            get_converter(unit_phpa, hPa),
+            get_converter(unit_tc, deg_C),
+            get_converter(_d(unit_rh), one),
+            get_converter(unit_wl, micron),
+        ], (astrom_unit(), radian, None)
+
+
 def helper_apio(
     f,
     unit_sp,
@@ -423,49 +474,50 @@ def helper_atoiq(f, unit_type, unit_ri, unit_di, unit_astrom):
 
 
 def get_erfa_helpers():
-    ERFA_HELPERS = {}
-    ERFA_HELPERS[erfa_ufunc.s2c] = helper_s2c
-    ERFA_HELPERS[erfa_ufunc.s2p] = helper_s2p
-    ERFA_HELPERS[erfa_ufunc.c2s] = helper_c2s
-    ERFA_HELPERS[erfa_ufunc.p2s] = helper_p2s
-    ERFA_HELPERS[erfa_ufunc.pm] = helper_invariant
-    ERFA_HELPERS[erfa_ufunc.cpv] = helper_invariant
-    ERFA_HELPERS[erfa_ufunc.p2pv] = helper_p2pv
-    ERFA_HELPERS[erfa_ufunc.pv2p] = helper_pv2p
-    ERFA_HELPERS[erfa_ufunc.pv2s] = helper_pv2s
-    ERFA_HELPERS[erfa_ufunc.pvdpv] = helper_pv_multiplication
-    ERFA_HELPERS[erfa_ufunc.pvxpv] = helper_pv_multiplication
-    ERFA_HELPERS[erfa_ufunc.pvm] = helper_pvm
-    ERFA_HELPERS[erfa_ufunc.pvmpv] = helper_twoarg_invariant
-    ERFA_HELPERS[erfa_ufunc.pvppv] = helper_twoarg_invariant
-    ERFA_HELPERS[erfa_ufunc.pvstar] = helper_pvstar
-    ERFA_HELPERS[erfa_ufunc.pvtob] = helper_pvtob
-    ERFA_HELPERS[erfa_ufunc.pvu] = helper_pvu
-    ERFA_HELPERS[erfa_ufunc.pvup] = helper_pvup
-    ERFA_HELPERS[erfa_ufunc.pdp] = helper_multiplication
-    ERFA_HELPERS[erfa_ufunc.pxp] = helper_multiplication
-    ERFA_HELPERS[erfa_ufunc.rxp] = helper_multiplication
-    ERFA_HELPERS[erfa_ufunc.rxpv] = helper_multiplication
-    ERFA_HELPERS[erfa_ufunc.s2pv] = helper_s2pv
-    ERFA_HELPERS[erfa_ufunc.s2xpv] = helper_s2xpv
-    ERFA_HELPERS[erfa_ufunc.starpv] = helper_starpv
-    ERFA_HELPERS[erfa_ufunc.sxpv] = helper_multiplication
-    ERFA_HELPERS[erfa_ufunc.trxpv] = helper_multiplication
-    ERFA_HELPERS[erfa_ufunc.gc2gd] = helper_gc2gd
-    ERFA_HELPERS[erfa_ufunc.gc2gde] = helper_gc2gde
-    ERFA_HELPERS[erfa_ufunc.gd2gc] = helper_gd2gc
-    ERFA_HELPERS[erfa_ufunc.gd2gce] = helper_gd2gce
-    ERFA_HELPERS[erfa_ufunc.ldn] = helper_ldn
-    ERFA_HELPERS[erfa_ufunc.aper] = helper_aper
-    ERFA_HELPERS[erfa_ufunc.apio] = helper_apio
-    ERFA_HELPERS[erfa_ufunc.atciq] = helper_atciq
-    ERFA_HELPERS[erfa_ufunc.atciqn] = helper_atciqn
-    ERFA_HELPERS[erfa_ufunc.atciqz] = helper_atciqz_aticq
-    ERFA_HELPERS[erfa_ufunc.aticq] = helper_atciqz_aticq
-    ERFA_HELPERS[erfa_ufunc.aticqn] = helper_aticqn
-    ERFA_HELPERS[erfa_ufunc.atioq] = helper_atioq
-    ERFA_HELPERS[erfa_ufunc.atoiq] = helper_atoiq
-    return ERFA_HELPERS
+    return {
+        erfa_ufunc.apco13: helper_apco13,
+        erfa_ufunc.aper: helper_aper,
+        erfa_ufunc.apio: helper_apio,
+        erfa_ufunc.atciq: helper_atciq,
+        erfa_ufunc.atciqn: helper_atciqn,
+        erfa_ufunc.atciqz: helper_atciqz_aticq,
+        erfa_ufunc.aticq: helper_atciqz_aticq,
+        erfa_ufunc.aticqn: helper_aticqn,
+        erfa_ufunc.atioq: helper_atioq,
+        erfa_ufunc.atoiq: helper_atoiq,
+        erfa_ufunc.c2s: helper_c2s,
+        erfa_ufunc.cpv: helper_invariant,
+        erfa_ufunc.gc2gd: helper_gc2gd,
+        erfa_ufunc.gc2gde: helper_gc2gde,
+        erfa_ufunc.gd2gc: helper_gd2gc,
+        erfa_ufunc.gd2gce: helper_gd2gce,
+        erfa_ufunc.ldn: helper_ldn,
+        erfa_ufunc.p2pv: helper_p2pv,
+        erfa_ufunc.p2s: helper_p2s,
+        erfa_ufunc.pdp: helper_multiplication,
+        erfa_ufunc.pm: helper_invariant,
+        erfa_ufunc.pv2p: helper_pv2p,
+        erfa_ufunc.pv2s: helper_pv2s,
+        erfa_ufunc.pvdpv: helper_pv_multiplication,
+        erfa_ufunc.pvm: helper_pvm,
+        erfa_ufunc.pvmpv: helper_twoarg_invariant,
+        erfa_ufunc.pvppv: helper_twoarg_invariant,
+        erfa_ufunc.pvstar: helper_pvstar,
+        erfa_ufunc.pvtob: helper_pvtob,
+        erfa_ufunc.pvu: helper_pvu,
+        erfa_ufunc.pvup: helper_pvup,
+        erfa_ufunc.pvxpv: helper_pv_multiplication,
+        erfa_ufunc.pxp: helper_multiplication,
+        erfa_ufunc.rxp: helper_multiplication,
+        erfa_ufunc.rxpv: helper_multiplication,
+        erfa_ufunc.s2c: helper_s2c,
+        erfa_ufunc.s2p: helper_s2p,
+        erfa_ufunc.s2pv: helper_s2pv,
+        erfa_ufunc.s2xpv: helper_s2xpv,
+        erfa_ufunc.starpv: helper_starpv,
+        erfa_ufunc.sxpv: helper_multiplication,
+        erfa_ufunc.trxpv: helper_multiplication,
+    }
 
 
 UFUNC_HELPERS.register_module("erfa.ufunc", erfa_ufuncs, get_erfa_helpers)
diff --git a/docs/table/construct_table.rst b/docs/table/construct_table.rst
index 445add1d3f..b77746877f 100644
--- a/docs/table/construct_table.rst
+++ b/docs/table/construct_table.rst
@@ -201,10 +201,8 @@ column where each row element is itself a two-element array.
     ...
   KeyError: 'a_new'
 
-.. _Row data:
-
-List of Rows
-------------
+Row Data
+--------
 
 Row-oriented data can be used to create a table using the ``rows``
 keyword argument.


Test patch:
diff --git a/astropy/coordinates/tests/test_iau_fullstack.py b/astropy/coordinates/tests/test_iau_fullstack.py
index 78280ceab1..db9fef6173 100644
--- a/astropy/coordinates/tests/test_iau_fullstack.py
+++ b/astropy/coordinates/tests/test_iau_fullstack.py
@@ -6,7 +6,7 @@ import pytest
 from numpy import testing as npt
 
 from astropy import units as u
-from astropy.coordinates import Angle, EarthLocation, SkyCoord, golden_spiral_grid
+from astropy.coordinates import EarthLocation, SkyCoord, golden_spiral_grid
 from astropy.coordinates.builtin_frames import ICRS, AltAz
 from astropy.coordinates.builtin_frames.utils import get_jd12
 from astropy.tests.helper import assert_quantity_allclose
@@ -58,21 +58,15 @@ def test_iau_fullstack(
     """
     obsconditions = DEFAULT_OBSCONDITIONS | fullstack_obsconditions
 
-    # create the altaz frame
-    altazframe = AltAz(
-        obstime=fullstack_times, location=fullstack_locations, **obsconditions
+    # Transform to the altaz frame
+    aacoo = FULLSTACK_ICRS.transform_to(
+        AltAz(obstime=fullstack_times, location=fullstack_locations, **obsconditions)
     )
 
-    aacoo = FULLSTACK_ICRS.transform_to(altazframe)
-
     # compare aacoo to the fiducial AltAz - should always be different
     fullstack_fiducial_altaz = FULLSTACK_ICRS.transform_to(ALTAZFRAME)
-    assert np.all(
-        np.abs(aacoo.alt - fullstack_fiducial_altaz.alt) > 50 * u.milliarcsecond
-    )
-    assert np.all(
-        np.abs(aacoo.az - fullstack_fiducial_altaz.az) > 50 * u.milliarcsecond
-    )
+    assert np.all(np.abs(aacoo.alt - fullstack_fiducial_altaz.alt) > 50 * u.mas)
+    assert np.all(np.abs(aacoo.az - fullstack_fiducial_altaz.az) > 50 * u.mas)
 
     # if the refraction correction is included, we *only* do the comparisons
     # where altitude is high enough.  The SOFA guides imply that below 5 deg is
@@ -94,22 +88,21 @@ def test_iau_fullstack(
     # check that we're consistent with the ERFA alt/az result
     astrom, _ = erfa.apco13(
         *get_jd12(fullstack_times, "utc"),
-        fullstack_times.delta_ut1_utc,
-        (geodetic := fullstack_locations.geodetic).lon.rad,
-        geodetic.lat.rad,
-        geodetic.height.to_value(u.m),
-        *Angle(iers.earth_orientation_table.get().pm_xy(fullstack_times)).rad,
-        altazframe.pressure.to_value(u.hPa),
-        altazframe.temperature.to_value(u.deg_C),
-        altazframe.relative_humidity.to_value(u.one),
-        altazframe.obswl.to_value(u.μm),
+        (iers_table := iers.earth_orientation_table.get()).ut1_utc(fullstack_times),
+        (geodetic := fullstack_locations.geodetic).lon,
+        geodetic.lat,
+        geodetic.height,
+        *iers_table.pm_xy(fullstack_times),
+        obsconditions["pressure"],
+        obsconditions["temperature"],
+        obsconditions["relative_humidity"],
+        obsconditions["obswl"],
     )
     erfa_az, erfa_zen, _, _, _ = erfa.atioq(
-        *erfa.atciq(FULLSTACK_ICRS.ra.rad, FULLSTACK_ICRS.dec.rad, 0, 0, 0, 0, astrom),
-        astrom,
+        *erfa.atciqz(FULLSTACK_ICRS.ra, FULLSTACK_ICRS.dec, astrom), astrom
     )
-    npt.assert_allclose(aacoo.alt.rad, np.pi / 2 - erfa_zen, atol=1e-7)
-    npt.assert_allclose(aacoo.az.rad, erfa_az, atol=1e-7)
+    assert_quantity_allclose(aacoo.alt, 90 * u.deg - erfa_zen, atol=1 * u.mas)
+    assert_quantity_allclose(aacoo.az, erfa_az, atol=1 * u.mas)
 
 
 def test_fiducial_roudtrip():
diff --git a/astropy/units/tests/test_quantity_erfa_ufuncs.py b/astropy/units/tests/test_quantity_erfa_ufuncs.py
index a55dc7c970..1dbb936263 100644
--- a/astropy/units/tests/test_quantity_erfa_ufuncs.py
+++ b/astropy/units/tests/test_quantity_erfa_ufuncs.py
@@ -10,6 +10,9 @@ from numpy.testing import assert_array_equal
 
 from astropy import units as u
 from astropy.tests.helper import assert_quantity_allclose
+from astropy.units.quantity_helper.erfa import astrom_unit
+
+ASTROM_UNIT = astrom_unit()
 
 
 def vvd(val, valok, dval, func, test, status):
@@ -18,16 +21,17 @@ def vvd(val, valok, dval, func, test, status):
 
 
 class TestPVUfuncs:
-    def setup_class(self):
-        self.pv_unit = u.Unit("AU,AU/day")
-        self.pv_value = np.array(
+    @classmethod
+    def setup_class(cls):
+        cls.pv_unit = u.Unit("AU,AU/day")
+        cls.pv_value = np.array(
             [
                 ([1.0, 0.0, 0.0], [0.0, 0.0125, 0.0]),
                 ([0.0, 1.0, 0.0], [-0.0125, 0.0, 0.0]),
             ],
             dtype=erfa_ufunc.dt_pv,
         )
-        self.pv = self.pv_value << self.pv_unit
+        cls.pv = cls.pv_value << cls.pv_unit
 
     def test_cpv(self):
         pv_copy = erfa_ufunc.cpv(self.pv)
@@ -318,7 +322,8 @@ class TestPVUfuncs:
 
 
 class TestEraStructUfuncs:
-    def setup_class(self):
+    @classmethod
+    def setup_class(cls):
         ldbody = np.array(
             [
                 (0.00028574, 3e-10, ([-7.81014427, -5.60956681, -1.98079819],
@@ -331,22 +336,19 @@ class TestEraStructUfuncs:
             dtype=erfa_ufunc.dt_eraLDBODY
         )  # fmt: skip
         ldbody_unit = u.StructuredUnit("Msun,radian,(AU,AU/day)", ldbody.dtype)
-        self.ldbody = ldbody << ldbody_unit
-        self.ob = [-0.974170437, -0.2115201, -0.0917583114] << u.AU
-        self.sc = np.array([-0.763276255, -0.608633767, -0.216735543])
+        cls.ldbody = ldbody << ldbody_unit
+        cls.ob = [-0.974170437, -0.2115201, -0.0917583114] << u.AU
+        cls.sc = np.array([-0.763276255, -0.608633767, -0.216735543])
 
         # From t_atciq in t_erfa_c.c
         astrom, eo = erfa_ufunc.apci13(2456165.5, 0.401182685)
-        self.astrom_unit = u.StructuredUnit(
-            "yr,AU,1,AU,1,1,1,rad,rad,rad,rad,1,1,1,rad,rad,rad", astrom.dtype
-        )
-        self.astrom = astrom << self.astrom_unit
-        self.rc = 2.71 * u.rad
-        self.dc = 0.174 * u.rad
-        self.pr = 1e-5 * u.rad / u.year
-        self.pd = 5e-6 * u.rad / u.year
-        self.px = 0.1 * u.arcsec
-        self.rv = 55.0 * u.km / u.s
+        cls.astrom = astrom << ASTROM_UNIT
+        cls.rc = 2.71 * u.rad
+        cls.dc = 0.174 * u.rad
+        cls.pr = 1e-5 * u.rad / u.year
+        cls.pd = 5e-6 * u.rad / u.year
+        cls.px = 0.1 * u.arcsec
+        cls.rv = 55.0 * u.km / u.s
 
     def test_ldn_basic(self):
         sn = erfa_ufunc.ldn(self.ldbody, self.ob, self.sc)
@@ -463,7 +465,7 @@ class TestEraStructUfuncs:
             0.59,
             0.55,
         )
-        astrom = astrom << self.astrom_unit
+        astrom = astrom << ASTROM_UNIT
 
         ri = 2.710121572969038991 * u.rad
         di = 0.1729371367218230438 * u.rad
@@ -489,62 +491,134 @@ class TestEraStructUfuncs:
             di2, 0.17293718391166087785 * u.rad, atol=1e-12 * u.rad
         )
 
-    def test_apio(self):
-        sp = -3.01974337e-11 * u.rad
-        theta = 3.14540971 * u.rad
-        elong = -0.527800806 * u.rad
-        phi = -1.2345856 * u.rad
-        hm = 2738.0 * u.m
-        xp = 2.47230737e-7 * u.rad
-        yp = 1.82640464e-6 * u.rad
-        refa = 0.000201418779 * u.rad
-        refb = -2.36140831e-7 * u.rad
-        astrom = erfa_ufunc.apio(
-            sp.to(u.deg), theta, elong, phi, hm.to(u.km),
-            xp.to(u.arcsec), yp, refa, refb.to(u.deg),
-        )  # fmt: skip
-        assert astrom.unit == self.astrom_unit
-        for name, value in [
-            ("along", -0.5278008060295995734),
-            ("xpl", 0.1133427418130752958e-5),
-            ("ypl", 0.1453347595780646207e-5),
-            ("sphi", -0.9440115679003211329),
-            ("cphi", 0.3299123514971474711),
-            ("diurab", 0.5135843661699913529e-6),
-            ("eral", 2.617608903970400427),
-            ("refa", 0.2014187790000000000e-3),
-            ("refb", -0.2361408310000000000e-6),
-        ]:
-            assert_quantity_allclose(
-                astrom[name],
-                value * self.astrom_unit[name],
-                rtol=1e-12,
-                atol=0 * self.astrom_unit[name],
+
+class TestAp:
+    @classmethod
+    def setup_class(cls):
+        # Values from apco13 and apio test cases in t_erfa_c.c
+        # with units changed so we can check conversion.
+        cls.utc1 = 2456384.5
+        cls.utc2 = 0.969254051
+        cls.dut1 = (0.1550675 * u.s).to(u.ms)
+        cls.sp = (-3.01974337e-11 * u.rad).to(u.deg)
+        cls.theta = (3.14540971 * u.rad).to(u.mdeg)
+        cls.elong = (-0.527800806 * u.rad).to(u.Marcsec)
+        cls.phi = (-1.2345856 * u.rad).to(u.arcmin)
+        cls.hm = (2738.0 * u.m).to(u.km)
+        cls.phpa = (731.0 * u.hPa).to(u.Pa)
+        cls.tc = (12.8 * u.deg_C).to(u.K, equivalencies=u.temperature())
+        cls.rh = (0.59 * u.one).to(u.percent)
+        cls.wl = (0.55 * u.micron).to(u.AA)
+        # For apio.
+        cls.xp = (2.47230737e-7 * u.rad).to(u.arcsec)
+        cls.yp = (1.82640464e-6 * u.rad).to(u.arcmin)
+        cls.refa = (0.000201418779 * u.rad).to(u.mas)
+        cls.refb = (-2.36140831e-7 * u.rad).to(u.uas)
+        cls.apco13_args = [
+            getattr(cls, name)
+            for name in "utc1 utc2 dut1 elong phi hm xp yp phpa tc rh wl".split()
+        ]
+        cls.apio_args = [
+            getattr(cls, name)
+            for name in "sp theta elong phi hm xp yp refa refb".split()
+        ]
+
+    def test_apco13(self):
+        astrom, eo, status = erfa_ufunc.apco13(*self.apco13_args)
+        assert status == 0
+        vvd(eo, -0.003020548354802412839, 1e-14, "eraApco13", "eo", status)
+        assert astrom.unit == ASTROM_UNIT
+        for name, expected in {
+            "pmt": 13.25248468622475727,
+            "eb": [-0.9741827107320875162,
+                   -0.2115130190489716682,
+                   -0.09179840189496755339],
+            "eh": [-0.9736425572586935247,
+                   -0.2092452121603336166,
+                   -0.09075578153885665295],
+            "em": 0.9998233240913898141,
+            "v": [0.2078704994520489246e-4,
+                  -0.8955360133238868938e-4,
+                  -0.3863338993055887398e-4],
+            "bm1": 0.9999999950277561004,
+            "bpn": np.array(
+                [[ 0.9999991390295147999,
+                   0.4978650075315529277e-7,
+                   0.001312227200850293372],
+                 [-0.1136336652812486604e-7,
+                   0.9999999995713154865,
+                  -0.2928086230975367296e-4],
+                 [-0.001312227201745553566,
+                   0.2928082218847679162e-4,
+                   0.9999991386008312212],
+                ]).T,
+            "along": -0.5278008060295995733,
+            "xpl": 0.1133427418130752958e-5,
+            "ypl": 0.1453347595780646207e-5,
+            "sphi": -0.9440115679003211329,
+            "cphi": 0.3299123514971474711,
+            "diurab": 0,
+            "eral": 2.617608909189664000,
+            "refa": 0.2014187785940396921e-3,
+            "refb": -0.2361408314943696227e-6,
+        }.items():  # fmt: skip
+            assert_quantity_allclose(astrom[name], expected * ASTROM_UNIT[name])
+
+    def test_apco13_no_time_units(self):
+        msg = "cannot pass in units for 2-part time in apco13"
+        with pytest.raises(TypeError, match=msg):
+            erfa_ufunc.apco13(self.utc1 * u.day, *self.apco13_args[1:])
+
+        with pytest.raises(TypeError, match=msg):
+            erfa_ufunc.apco13(self.utc1, self.utc2 * u.day, *self.apco13_args[2:])
+
+        with pytest.raises(TypeError, match=msg):
+            erfa_ufunc.apco13(
+                self.utc1 * u.day, self.utc2 * u.day, *self.apco13_args[2:]
             )
 
+    def test_apio(self):
+        astrom = erfa_ufunc.apio(*self.apio_args)
+        assert astrom.unit == ASTROM_UNIT
+        for name, expected in {
+            "along": -0.5278008060295995734,
+            "xpl": 0.1133427418130752958e-5,
+            "ypl": 0.1453347595780646207e-5,
+            "sphi": -0.9440115679003211329,
+            "cphi": 0.3299123514971474711,
+            "diurab": 0.5135843661699913529e-6,
+            "eral": 2.617608903970400427,
+            "refa": 0.2014187790000000000e-3,
+            "refb": -0.2361408310000000000e-6,
+        }.items():
+            assert_quantity_allclose(astrom[name], expected * ASTROM_UNIT[name])
+
 
 class TestGeodetic:
-    def setup_class(self):
-        self.ellipsoid = 1
-        self.length_unit = u.Unit("m")
-        self.equatorial_radius_value = 6378136.0
-        self.equatorial_radius = self.equatorial_radius_value << self.length_unit
-        self.flattening = 0.0033528 * u.dimensionless_unscaled
-        self.lon_value = 0.9827937232473290680
-        self.lon_unit = u.Unit("rad")
-        self.lon = self.lon_value << self.lon_unit
-        self.lat_value = 0.9716018377570411532
-        self.lat_unit = u.Unit("rad")
-        self.lat = self.lat_value << self.lat_unit
-        self.height_value = 332.36862495764397
-        self.height = self.height_value << self.length_unit
-        self.x_value = 2e6
-        self.x = self.x_value << self.length_unit
-        self.y_value = 3e6
-        self.y = self.y_value << self.length_unit
-        self.z_value = 5.244e6
-        self.z = self.z_value << self.length_unit
-        self.xyz = np.stack([self.x, self.y, self.z])
+    @classmethod
+    def setup_class(cls):
+        cls.ellipsoid = 1
+        cls.length_unit = u.Unit("m")
+        cls.equatorial_radius_value = 6378136.0
+        cls.equatorial_radius = (cls.equatorial_radius_value << cls.length_unit).to(
+            u.km
+        )
+        cls.flattening = 0.0033528 * u.dimensionless_unscaled
+        cls.lon_value = 0.9827937232473290680
+        cls.lon_unit = u.Unit("rad")
+        cls.lon = (cls.lon_value << cls.lon_unit).to(u.deg)
+        cls.lat_value = 0.9716018377570411532
+        cls.lat_unit = u.Unit("rad")
+        cls.lat = (cls.lat_value << cls.lat_unit).to(u.deg)
+        cls.height_value = 332.36862495764397
+        cls.height = cls.height_value << cls.length_unit
+        cls.x_value = 2e6
+        cls.x = (cls.x_value << cls.length_unit).to(u.cm)
+        cls.y_value = 3e6
+        cls.y = (cls.y_value << cls.length_unit).to(u.km)
+        cls.z_value = 5.244e6
+        cls.z = (cls.z_value << cls.length_unit).to(u.Mm)
+        cls.xyz = np.stack([cls.x, cls.y, cls.z])
 
     def test_unit_errors(self):
         """Test unit errors when dimensionless parameters are used"""
@@ -585,11 +659,10 @@ class TestGeodetic:
     def test_gc2gde(self):
         """Test that we reproduce erfa/src/t_erfa_c.c t_gc2gd"""
 
-        status = 0
         e, p, h, status = erfa_ufunc.gc2gde(
             self.equatorial_radius, self.flattening, self.xyz
         )
-
+        assert status == 0
         vvd(e, self.lon_value, 1e-14, "eraGc2gde", "e", status)
         vvd(p, self.lat_value, 1e-14, "eraGc2gde", "p", status)
         vvd(h, self.height_value, 1e-8, "eraGc2gde", "h", status)
@@ -597,11 +670,10 @@ class TestGeodetic:
     def test_gd2gce(self):
         """Test that we reproduce erfa/src/t_erfa_c.c t_gc2gd"""
 
-        status = 0
         xyz, status = erfa_ufunc.gd2gce(
             self.equatorial_radius, self.flattening, self.lon, self.lat, self.height
         )
-
+        assert status == 0
         vvd(xyz[0], self.x_value, 1e-7, "eraGd2gce", "e", status)
         vvd(xyz[1], self.y_value, 1e-7, "eraGd2gce", "p", status)
         vvd(xyz[2], self.z_value, 1e-7, "eraGd2gce", "h", status)