diff --git a/glass/algorithm.py b/glass/algorithm.py
index 9f39a594..2f02f4d3 100644
--- a/glass/algorithm.py
+++ b/glass/algorithm.py
@@ -67,28 +67,28 @@ def nnls(
         maxiter = 3 * n
 
     index = np.arange(n)
-    p = np.full(n, fill_value=False)
+    q = np.full(n, fill_value=False)
     x = np.zeros(n)
     for _ in range(maxiter):
-        if np.all(p):
+        if np.all(q):
             break
         w = np.dot(b - a @ x, a)
-        m = index[~p][np.argmax(w[~p])]
+        m = index[~q][np.argmax(w[~q])]
         if w[m] <= tol:
             break
-        p[m] = True
+        q[m] = True
         while True:
-            ap = a[:, p]
-            xp = x[p]
-            sp = np.linalg.solve(ap.T @ ap, b @ ap)
-            t = sp <= 0
+            aq = a[:, q]
+            xq = x[q]
+            sq = np.linalg.solve(aq.T @ aq, b @ aq)
+            t = sq <= 0
             if not np.any(t):
                 break
-            alpha = -np.min(xp[t] / (xp[t] - sp[t]))
-            x[p] += alpha * (sp - xp)
-            p[x <= 0] = False
-        x[p] = sp
-        x[~p] = 0
+            alpha = -np.min(xq[t] / (xq[t] - sq[t]))
+            x[q] += alpha * (sq - xq)
+            q[x <= 0] = False
+        x[q] = sq
+        x[~q] = 0
     return x
 
 
diff --git a/glass/arraytools.py b/glass/arraytools.py
index 28d37b4b..f52925af 100644
--- a/glass/arraytools.py
+++ b/glass/arraytools.py
@@ -91,8 +91,8 @@ def broadcast_leading_axes(
 
 def ndinterp(  # noqa: PLR0913
     x: float | NDArray[np.float64],
-    xp: Sequence[float] | NDArray[np.float64],
-    fp: Sequence[float] | NDArray[np.float64],
+    xq: Sequence[float] | NDArray[np.float64],
+    fq: Sequence[float] | NDArray[np.float64],
     axis: int = -1,
     left: float | None = None,
     right: float | None = None,
@@ -105,16 +105,16 @@ def ndinterp(  # noqa: PLR0913
     ----------
     x
         The x-coordinates.
-    xp
+    xq
         The x-coordinates of the data points.
-    fp
-        The function values corresponding to the x-coordinates in *xp*.
+    fq
+        The function values corresponding to the x-coordinates in *xq*.
     axis
         The axis to interpolate over.
     left
-        The value to return for x < xp[0].
+        The value to return for x < xq[0].
     right
-        The value to return for x > xp[-1].
+        The value to return for x > xq[-1].
     period
         The period of the function, used for interpolating periodic data.
 
@@ -124,9 +124,9 @@ def ndinterp(  # noqa: PLR0913
 
     """
     return np.apply_along_axis(
-        partial(np.interp, x, xp),
+        partial(np.interp, x, xq),
         axis,
-        fp,
+        fq,
         left=left,
         right=right,
         period=period,
diff --git a/tests/test_arraytools.py b/tests/test_arraytools.py
index a26faaff..695680c2 100644
--- a/tests/test_arraytools.py
+++ b/tests/test_arraytools.py
@@ -57,40 +57,40 @@ def test_broadcast_leading_axes() -> None:
 def test_ndinterp() -> None:
     # test 1d interpolation
 
-    xp = np.array([0, 1, 2, 3, 4])
-    yp = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
+    xq = np.array([0, 1, 2, 3, 4])
+    yq = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
 
     x: float | NDArray[np.float64] = 0.5
-    y = glass.arraytools.ndinterp(x, xp, yp)
+    y = glass.arraytools.ndinterp(x, xq, yq)
     assert np.shape(y) == ()
     np.testing.assert_allclose(y, 1.15, atol=1e-15)
 
     x = np.array([0.5, 1.5, 2.5])
-    y = glass.arraytools.ndinterp(x, xp, yp)
+    y = glass.arraytools.ndinterp(x, xq, yq)
     assert np.shape(y) == (3,)
     np.testing.assert_allclose(y, [1.15, 1.25, 1.35], atol=1e-15)
 
     x = np.array([[0.5, 1.5], [2.5, 3.5]])
-    y = glass.arraytools.ndinterp(x, xp, yp)
+    y = glass.arraytools.ndinterp(x, xq, yq)
     assert np.shape(y) == (2, 2)
     np.testing.assert_allclose(y, [[1.15, 1.25], [1.35, 1.45]], atol=1e-15)
 
     # test nd interpolation in final axis
 
-    yp = np.array([[1.1, 1.2, 1.3, 1.4, 1.5], [2.1, 2.2, 2.3, 2.4, 2.5]])
+    yq = np.array([[1.1, 1.2, 1.3, 1.4, 1.5], [2.1, 2.2, 2.3, 2.4, 2.5]])
 
     x = 0.5
-    y = glass.arraytools.ndinterp(x, xp, yp)
+    y = glass.arraytools.ndinterp(x, xq, yq)
     assert np.shape(y) == (2,)
     np.testing.assert_allclose(y, [1.15, 2.15], atol=1e-15)
 
     x = np.array([0.5, 1.5, 2.5])
-    y = glass.arraytools.ndinterp(x, xp, yp)
+    y = glass.arraytools.ndinterp(x, xq, yq)
     assert np.shape(y) == (2, 3)
     np.testing.assert_allclose(y, [[1.15, 1.25, 1.35], [2.15, 2.25, 2.35]], atol=1e-15)
 
     x = np.array([[0.5, 1.5], [2.5, 3.5]])
-    y = glass.arraytools.ndinterp(x, xp, yp)
+    y = glass.arraytools.ndinterp(x, xq, yq)
     assert np.shape(y) == (2, 2, 2)
     np.testing.assert_allclose(
         y,
@@ -100,17 +100,17 @@ def test_ndinterp() -> None:
 
     # test nd interpolation in middle axis
 
-    yp = np.array(
+    yq = np.array(
         [[[1.1], [1.2], [1.3], [1.4], [1.5]], [[2.1], [2.2], [2.3], [2.4], [2.5]]],
     )
 
     x = 0.5
-    y = glass.arraytools.ndinterp(x, xp, yp, axis=1)
+    y = glass.arraytools.ndinterp(x, xq, yq, axis=1)
     assert np.shape(y) == (2, 1)
     np.testing.assert_allclose(y, [[1.15], [2.15]], atol=1e-15)
 
     x = np.array([0.5, 1.5, 2.5])
-    y = glass.arraytools.ndinterp(x, xp, yp, axis=1)
+    y = glass.arraytools.ndinterp(x, xq, yq, axis=1)
     assert np.shape(y) == (2, 3, 1)
     np.testing.assert_allclose(
         y,
@@ -119,7 +119,7 @@ def test_ndinterp() -> None:
     )
 
     x = np.array([[0.5, 1.5, 2.5, 3.5], [3.5, 2.5, 1.5, 0.5], [0.5, 3.5, 1.5, 2.5]])
-    y = glass.arraytools.ndinterp(x, xp, yp, axis=1)
+    y = glass.arraytools.ndinterp(x, xq, yq, axis=1)
     assert np.shape(y) == (2, 3, 4, 1)
     np.testing.assert_allclose(
         y,