diff --git a/AGENTS.md b/AGENTS.md
index cf84d70..fd5d9b0 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -27,7 +27,9 @@ Non-obvious details:
   is a plain class. *Gen 2* is `AbstractSingleInterpolator` → `SingleGridInterpolator` /
   `SingleSplineInterpolator` and `AbstractMultipleInterpolator` → `MultipleGridInterpolator` /
   `MultipleSplineInterpolator`: functional, no mutation, and they accept `chunks=` for lazy dask
-  output. **New work belongs in Gen 2.**
+  output. **New work belongs in Gen 2.** `AbstractMultipleInterpolator.interpolate` — and so
+  `interpolate_to_shape`, both `Multiple*Interpolator` classes, and the `Geo*Interpolator` classes
+  built on them — returns a **tuple**, one array per input array, in input order.
 - **`GeoGridInterpolator` and `GeoSplineInterpolator` are built dynamically** by the
   `_work_with_lonlats(klass)` class factory (`geointerpolator.py:78`). Grepping for the class body only
   finds `GeoKlass`; the names are bound at `geointerpolator.py:111-112`.
@@ -177,10 +179,6 @@ Verified as of this writing; fix them only when the task calls for it.
   `_modis_utils.pyx`) and `6370.997` km (`_modis_interpolator.pyx`) are the same value in different
   units; `6371008.7714` (`viiinterpolator.py`) is the IUGG mean radius, a genuinely different number,
   so unifying it would change `viiinterpolator` results. Deliberate, not an oversight.
-- `AbstractMultipleInterpolator.interpolate` (inherited by both `Multiple*Interpolator` classes)
-  returns a **generator**, not a tuple.
-- `interpolator.py`'s `Interpolator` docstring describes `kx_`/`ky_` as orders "in x and y", but
-  `_interp` passes `kx=self.kx_` to the *row* (first) axis of `RectBivariateSpline`.
 
 ## Roadmap
 
diff --git a/geotiepoints/geointerpolator.py b/geotiepoints/geointerpolator.py
index 356e780..6e48dc9 100644
--- a/geotiepoints/geointerpolator.py
+++ b/geotiepoints/geointerpolator.py
@@ -19,11 +19,14 @@ class GeoInterpolator(Interpolator):
 
     The constructor takes in the tiepointed data as *data*, the
     *tiepoint_grid* and the desired *final_grid*. As optional arguments, one
-    can provide *kx_* and *ky_* as interpolation orders (in x and y directions
-    respectively), and the *chunksize* if the data has to be handled by pieces
-    along the y axis (this affects how the extrapolator behaves). If
-    *chunksize* is set, don't forget to adjust the interpolation orders
-    accordingly: the interpolation is indeed done globaly (not chunkwise).
+    can provide *kx_* and *ky_* as interpolation orders. Note that *kx_*
+    applies to the first (row) axis of the tiepoint grid, the along-track or y
+    direction, and *ky_* to the second (column) axis, the across-track or x
+    direction. One can also provide the *chunksize* if the data has to be
+    handled by pieces along the y axis (this affects how the extrapolator
+    behaves). If *chunksize* is set, don't forget to adjust the interpolation
+    orders accordingly: the interpolation is indeed done globaly (not
+    chunkwise).
 
     """
 
diff --git a/geotiepoints/interpolator.py b/geotiepoints/interpolator.py
index d4b2df4..c81cf64 100644
--- a/geotiepoints/interpolator.py
+++ b/geotiepoints/interpolator.py
@@ -61,10 +61,12 @@ class Interpolator:
     Uses numpy and scipy.
 
     The constructor takes in the tiepointed data as *data*, the *tiepoint_grid* and the desired *final_grid*. As
-    optional arguments, one can provide *kx_* and *ky_* as interpolation orders (in x and y directions respectively),
-    and the *chunksize* if the data has to be handled by pieces along the y axis (this affects how the extrapolator
-    behaves). If *chunksize* is set, don't forget to adjust the interpolation orders accordingly: the interpolation
-    is indeed done globaly (not chunkwise).
+    optional arguments, one can provide *kx_* and *ky_* as interpolation orders. Note that *kx_* applies to the
+    first (row) axis of the tiepoint grid, the along-track or y direction, and *ky_* to the second (column) axis,
+    the across-track or x direction. One can also provide the *chunksize* if the data
+    has to be handled by pieces along the y axis (this affects how the extrapolator behaves). If *chunksize* is
+    set, don't forget to adjust the interpolation orders accordingly: the interpolation is indeed done globaly
+    (not chunkwise).
     """
 
     def __init__(self, data, tiepoint_grid, final_grid,
@@ -352,8 +354,13 @@ def interpolate(self, fine_points, **kwargs):
         """Interpolate the data.
 
         The keyword arguments will be passed on to SingleGridInterpolator's interpolate function.
+
+        Returns:
+            A tuple with one interpolated array per input data array, in input order.
         """
-        return (interpolator.interpolate(fine_points, **kwargs) for interpolator in self.interpolators)
+        return tuple(
+            interpolator.interpolate(fine_points, **kwargs) for interpolator in self.interpolators
+        )
 
     def interpolate_to_shape(self, shape, **interpolator_call_kwargs):
         """Interpolate to a given *shape*."""