diff --git a/docs/user/how-glass-works.rst b/docs/user/how-glass-works.rst
index d873236e..59eec4b7 100644
--- a/docs/user/how-glass-works.rst
+++ b/docs/user/how-glass-works.rst
@@ -30,13 +30,13 @@ which are flat and non-overlapping.
 
 .. plot::
 
-    from glass.shells import redshift_grid, tophat_windows
+    import glass
 
     # create a redshift grid for shell edges
-    zs = redshift_grid(0., 0.5, dz=0.1)
+    zs = glass.redshift_grid(0., 0.5, dz=0.1)
 
     # create the top hat windows
-    ws = tophat_windows(zs)
+    ws = glass.tophat_windows(zs)
 
     # plot each window
     for i, (za, wa, zeff) in enumerate(ws):
@@ -74,17 +74,16 @@ included:
 
 .. plot::
 
-    from glass.shells import (redshift_grid, tophat_windows, linear_windows,
-                              cubic_windows)
+    import glass
 
-    plot_windows = [tophat_windows, linear_windows,
-                    cubic_windows]
+    plot_windows = [glass.tophat_windows, glass.linear_windows,
+                    glass.cubic_windows]
     nr = (len(plot_windows)+1)//2
 
     fig, axes = plt.subplots(nr, 2, figsize=(8, nr*3), layout="constrained",
                              squeeze=False, sharex=False, sharey=True)
 
-    zs = redshift_grid(0., 0.5, dz=0.1)
+    zs = glass.redshift_grid(0., 0.5, dz=0.1)
     zt = np.linspace(0., 0.5, 200)
 
     for ax in axes.flat:
diff --git a/glass/fields.py b/glass/fields.py
index f8a0b7fb..492bc620 100644
--- a/glass/fields.py
+++ b/glass/fields.py
@@ -13,7 +13,8 @@
 from gaussiancl import gaussiancl
 from transformcl import cltovar
 
-from glass import grf
+import glass
+import glass.grf
 
 if TYPE_CHECKING:
     from collections.abc import Callable, Generator, Iterable, Iterator, Sequence
@@ -21,9 +22,7 @@
 
     from numpy.typing import NDArray
 
-    from glass.shells import RadialWindow
-
-    Fields = Sequence[grf.Transformation]
+    Fields = Sequence[glass.grf.Transformation]
     Cls = Sequence[NDArray[Any]]
 
 
@@ -640,7 +639,7 @@ def effective_cls(
     return out
 
 
-def gaussian_fields(shells: Sequence[RadialWindow]) -> Sequence[grf.Normal]:
+def gaussian_fields(shells: Sequence[glass.RadialWindow]) -> Sequence[glass.grf.Normal]:
     """
     Create Gaussian random fields for radial windows *shells*.
 
@@ -654,13 +653,13 @@ def gaussian_fields(shells: Sequence[RadialWindow]) -> Sequence[grf.Normal]:
         A sequence describing the Gaussian random fields.
 
     """
-    return [grf.Normal() for _shell in shells]
+    return [glass.grf.Normal() for _shell in shells]
 
 
 def lognormal_fields(
-    shells: Sequence[RadialWindow],
+    shells: Sequence[glass.RadialWindow],
     shift: Callable[[float], float] | None = None,
-) -> Sequence[grf.Lognormal]:
+) -> Sequence[glass.grf.Lognormal]:
     """
     Create lognormal fields for radial windows *shells*.  If *shifts* is
     given, it must be a callable that returns a lognormal shift (i.e.
@@ -681,7 +680,7 @@ def lognormal_fields(
     if shift is None:
         shift = lambda _z: 1.0  # noqa: E731
 
-    return [grf.Lognormal(shift(shell.zeff)) for shell in shells]
+    return [glass.grf.Lognormal(shift(shell.zeff)) for shell in shells]
 
 
 def compute_gaussian_spectra(fields: Fields, spectra: Cls) -> Cls:
@@ -711,7 +710,7 @@ def compute_gaussian_spectra(fields: Fields, spectra: Cls) -> Cls:
 
     gls = []
     for i, j, cl in enumerate_spectra(spectra):
-        gl = grf.compute(cl, fields[i], fields[j])
+        gl = glass.grf.compute(cl, fields[i], fields[j])
         gls.append(gl)
     return gls
 
@@ -757,7 +756,7 @@ def solve_gaussian_spectra(fields: Fields, spectra: Cls) -> Cls:
             monopole = 0.0 if cl[0] == 0 else None
 
             # call solver
-            gl, _cl_out, info = grf.solve(cl, t1, t2, pad=pad, monopole=monopole)
+            gl, _cl_out, info = glass.grf.solve(cl, t1, t2, pad=pad, monopole=monopole)
 
             # warn if solver didn't converge
             if info == 0:
diff --git a/glass/galaxies.py b/glass/galaxies.py
index fcb2cd78..f63f2e3c 100644
--- a/glass/galaxies.py
+++ b/glass/galaxies.py
@@ -25,19 +25,18 @@
 import healpix
 import numpy as np
 
-from glass.core.array import broadcast_leading_axes, cumulative_trapezoid
+import glass
+import glass.core.array
 
 if TYPE_CHECKING:
     from numpy.typing import NDArray
 
     from cosmology import Cosmology
 
-    from glass.shells import RadialWindow
-
 
 def redshifts(
     n: int | NDArray[np.float64],
-    w: RadialWindow,
+    w: glass.RadialWindow,
     *,
     rng: np.random.Generator | None = None,
 ) -> NDArray[np.float64]:
@@ -118,7 +117,7 @@ def redshifts_from_nz(
         rng = np.random.default_rng()
 
     # bring inputs' leading axes into common shape
-    dims, *rest = broadcast_leading_axes((count, 0), (z, 1), (nz, 1))
+    dims, *rest = glass.core.array.broadcast_leading_axes((count, 0), (z, 1), (nz, 1))
     count_out, z_out, nz_out = rest
 
     # list of results for all dimensions
@@ -130,7 +129,7 @@ def redshifts_from_nz(
     # go through extra dimensions; also works if dims is empty
     for k in np.ndindex(dims):
         # compute the CDF of each galaxy population
-        cdf = cumulative_trapezoid(nz_out[k], z_out[k], dtype=float)
+        cdf = glass.core.array.cumulative_trapezoid(nz_out[k], z_out[k], dtype=float)
         cdf /= cdf[-1]
 
         # sample redshifts and store result
diff --git a/glass/grf/_solver.py b/glass/grf/_solver.py
index 86d57047..1fd72804 100644
--- a/glass/grf/_solver.py
+++ b/glass/grf/_solver.py
@@ -5,13 +5,13 @@
 import numpy as np
 from transformcl import cltocorr, corrtocl
 
+import glass.grf
+
 if TYPE_CHECKING:
     from typing import Any
 
     from numpy.typing import NDArray
 
-from glass.grf._core import Transformation, corr, dcorr, icorr
-
 
 def _relerr(dx: NDArray[Any], x: NDArray[Any]) -> float:
     """Compute the relative error max(|dx/x|)."""
@@ -21,8 +21,8 @@ def _relerr(dx: NDArray[Any], x: NDArray[Any]) -> float:
 
 def solve(  # noqa: PLR0912, PLR0913
     cl: NDArray[Any],
-    t1: Transformation,
-    t2: Transformation | None = None,
+    t1: glass.grf.Transformation,
+    t2: glass.grf.Transformation | None = None,
     *,
     pad: int = 0,
     initial: NDArray[Any] | None = None,
@@ -91,7 +91,7 @@ def solve(  # noqa: PLR0912, PLR0913
         raise ValueError(msg)
 
     if initial is None:
-        gl = corrtocl(icorr(t1, t2, cltocorr(cl)))
+        gl = corrtocl(glass.grf.icorr(t1, t2, cltocorr(cl)))
     else:
         gl = np.zeros(n)
         gl[: len(initial)] = initial[:n]
@@ -100,7 +100,7 @@ def solve(  # noqa: PLR0912, PLR0913
         gl[0] = monopole
 
     gt = cltocorr(np.pad(gl, (0, pad)))
-    rl = corrtocl(corr(t1, t2, gt))
+    rl = corrtocl(glass.grf.corr(t1, t2, gt))
     fl = rl[:n] - cl
     if monopole is not None:
         fl[0] = 0
@@ -116,7 +116,7 @@ def solve(  # noqa: PLR0912, PLR0913
             break
 
         ft = cltocorr(np.pad(fl, (0, pad)))
-        dt = dcorr(t1, t2, gt)
+        dt = glass.grf.dcorr(t1, t2, gt)
         xl = -corrtocl(ft / dt)[:n]
         if monopole is not None:
             xl[0] = 0
@@ -126,7 +126,7 @@ def solve(  # noqa: PLR0912, PLR0913
         while True:
             gl_ = gl + xl
             gt_ = cltocorr(np.pad(gl_, (0, pad)))
-            rl_ = corrtocl(corr(t1, t2, gt_))
+            rl_ = corrtocl(glass.grf.corr(t1, t2, gt_))
             fl_ = rl_[:n] - cl
             if monopole is not None:
                 fl_[0] = 0
diff --git a/glass/grf/_transformations.py b/glass/grf/_transformations.py
index 367fe818..8954e315 100644
--- a/glass/grf/_transformations.py
+++ b/glass/grf/_transformations.py
@@ -8,7 +8,7 @@
 
 from numpy.typing import NDArray  # noqa: TC002
 
-from glass.grf._core import corr, dcorr, icorr
+from glass.grf import corr, dcorr, icorr
 
 
 @dataclass
diff --git a/glass/lensing.py b/glass/lensing.py
index 65eda43e..5ccee900 100644
--- a/glass/lensing.py
+++ b/glass/lensing.py
@@ -43,7 +43,7 @@
 
     from cosmology import Cosmology
 
-    from glass.shells import RadialWindow
+    import glass
 
 
 @overload
@@ -351,7 +351,7 @@ def shear_from_convergence(
     transform.
 
     .. deprecated:: 2023.6
-       Use the more general :func:`from_convergence` function instead.
+       Use the more general :func:`glass.from_convergence` function instead.
 
     Parameters
     ----------
@@ -420,7 +420,7 @@ def __init__(self, cosmo: Cosmology) -> None:
         self.kappa2: NDArray[np.float64] | None = None
         self.kappa3: NDArray[np.float64] | None = None
 
-    def add_window(self, delta: NDArray[np.float64], w: RadialWindow) -> None:
+    def add_window(self, delta: NDArray[np.float64], w: glass.RadialWindow) -> None:
         """
         Add a mass plane from a window function to the convergence.
 
@@ -527,7 +527,7 @@ def wlens(self) -> float:
 
 
 def multi_plane_matrix(
-    shells: Sequence[RadialWindow],
+    shells: Sequence[glass.RadialWindow],
     cosmo: Cosmology,
 ) -> NDArray[np.float64]:
     """
@@ -555,7 +555,7 @@ def multi_plane_matrix(
 
 def multi_plane_weights(
     weights: NDArray[np.float64],
-    shells: Sequence[RadialWindow],
+    shells: Sequence[glass.RadialWindow],
     cosmo: Cosmology,
 ) -> NDArray[np.float64]:
     """
diff --git a/glass/observations.py b/glass/observations.py
index 835a51fc..702467eb 100644
--- a/glass/observations.py
+++ b/glass/observations.py
@@ -34,7 +34,7 @@
 import healpy as hp
 import numpy as np
 
-from glass.core.array import cumulative_trapezoid
+import glass.core.array
 
 if TYPE_CHECKING:
     from numpy.typing import NDArray
@@ -266,7 +266,7 @@ def equal_dens_zbins(
     # first compute the cumulative integral (by trapezoidal rule)
     # then normalise: the first z is at CDF = 0, the last z at CDF = 1
     # interpolate to find the z values at CDF = i/nbins for i = 0, ..., nbins
-    cuml_nz = cumulative_trapezoid(nz, z)
+    cuml_nz = glass.core.array.cumulative_trapezoid(nz, z)
     cuml_nz /= cuml_nz[[-1]]
     zbinedges = np.interp(np.linspace(0, 1, nbins + 1), cuml_nz, z)
 
diff --git a/glass/points.py b/glass/points.py
index 6f87d5ab..a6d81a30 100644
--- a/glass/points.py
+++ b/glass/points.py
@@ -36,15 +36,14 @@
 import healpix
 import numpy as np
 
-from glass.core.array import broadcast_first, broadcast_leading_axes, trapezoid_product
+import glass
+import glass.core.array
 
 if TYPE_CHECKING:
     from collections.abc import Callable, Generator
 
     from numpy.typing import NDArray
 
-    from glass.shells import RadialWindow
-
 
 ARCMIN2_SPHERE = 60**6 // 100 / np.pi
 
@@ -52,7 +51,7 @@
 def effective_bias(
     z: NDArray[np.float64],
     bz: NDArray[np.float64],
-    w: RadialWindow,
+    w: glass.RadialWindow,
 ) -> float | NDArray[np.double]:
     r"""
     Effective bias parameter from a redshift-dependent bias function.
@@ -86,7 +85,7 @@ def effective_bias(
 
     """
     norm = np.trapezoid(w.wa, w.za)
-    return trapezoid_product((z, bz), (w.za, w.wa)) / norm
+    return glass.core.array.trapezoid_product((z, bz), (w.za, w.wa)) / norm
 
 
 def linear_bias(
@@ -189,7 +188,8 @@ def positions_from_delta(  # noqa: PLR0912, PLR0913, PLR0915
     bias_model
         The bias model to apply. If a string, refers to a function in
         the :mod:`~glass.points` module, e.g. ``'linear'`` for
-        :func:`linear_bias()` or ``'loglinear'`` for :func:`loglinear_bias`.
+        :func:`glass.linear_bias()` or ``'glass.loglinear'`` for
+        :func:`glass.loglinear_bias`.
     remove_monopole
         If true, the monopole of the density contrast
         after biasing is fixed to zero.
@@ -232,7 +232,7 @@ def positions_from_delta(  # noqa: PLR0912, PLR0913, PLR0915
         inputs.append((bias, 0))
     if vis is not None:
         inputs.append((vis, 1))
-    dims, *rest = broadcast_leading_axes(*inputs)
+    dims, *rest = glass.core.array.broadcast_leading_axes(*inputs)
     ngal, delta, *rest = rest
     if bias is not None:
         bias, *rest = rest
@@ -409,7 +409,7 @@ def position_weights(
     """
     # bring densities and bias into the same shape
     if bias is not None:
-        densities, bias = broadcast_first(densities, bias)
+        densities, bias = glass.core.array.broadcast_first(densities, bias)
     # normalise densities after shape has been fixed
     densities = densities / np.sum(densities, axis=0)
     # apply bias after normalisation
diff --git a/glass/shells.py b/glass/shells.py
index 2a240759..49867191 100644
--- a/glass/shells.py
+++ b/glass/shells.py
@@ -53,8 +53,8 @@
 import numpy as np
 from numpy.typing import NDArray
 
-from glass.core.algorithm import nnls
-from glass.core.array import ndinterp
+import glass.core.algorithm
+import glass.core.array
 
 if TYPE_CHECKING:
     from cosmology import Cosmology
@@ -144,8 +144,8 @@ class RadialWindow(NamedTuple):
     immutable (however, the array entries may **not** be immutable; do
     not change them in place)::
 
-        >>> from glass import RadialWindow
-        >>> w1 = RadialWindow(..., ..., zeff=0.1)
+        >>> import glass
+        >>> w1 = glass.RadialWindow(..., ..., zeff=0.1)
         >>> w1.zeff = 0.15
         Traceback (most recent call last):
           File "<stdin>", line 1, in <module>
@@ -409,7 +409,9 @@ def restrict(
     """
     z_ = np.compress(np.greater(z, w.za[0]) & np.less(z, w.za[-1]), z)
     zr = np.union1d(w.za, z_)
-    fr = ndinterp(zr, z, f, left=0.0, right=0.0) * ndinterp(zr, w.za, w.wa)
+    fr = glass.core.array.ndinterp(
+        zr, z, f, left=0.0, right=0.0
+    ) * glass.core.array.ndinterp(zr, w.za, w.wa)
     return zr, fr
 
 
@@ -505,7 +507,7 @@ def partition(
     where :math:`\lambda` is a multiplier to enforce the integral
     constraints.
 
-    The :func:`partition()` function implements a number of methods to
+    The :func:`glass.partition()` function implements a number of methods to
     obtain a solution:
 
     If ``method="nnls"`` (the default), obtain a partition from a
@@ -520,7 +522,7 @@ def partition(
     contributions.
 
     If ``method="restrict"``, obtain a partition by integrating the
-    restriction (using :func:`restrict`) of the function :math:`f` to
+    restriction (using :func:`glass.restrict`) of the function :math:`f` to
     each window. For overlapping shells, this method might produce
     results which are far from the input function.
 
@@ -579,7 +581,7 @@ def partition_lstsq(
     a = a * dz
 
     # create the target vector of distribution values
-    b = ndinterp(zp, z, fz, left=0.0, right=0.0)
+    b = glass.core.array.ndinterp(zp, z, fz, left=0.0, right=0.0)
     b = b * dz
 
     # append a constraint for the integral
@@ -654,7 +656,7 @@ def partition_nnls(
     a = a * dz
 
     # create the target vector of distribution values
-    b = ndinterp(zp, z, fz, left=0.0, right=0.0)
+    b = glass.core.array.ndinterp(zp, z, fz, left=0.0, right=0.0)
     b = b * dz
 
     # append a constraint for the integral
@@ -673,7 +675,7 @@ def partition_nnls(
     # for each dim, find non-negative weights x such that y == r @ x
     x = np.empty([len(shells), *dims])
     for i in np.ndindex(*dims):
-        x[(..., *i)] = nnls(r, y[i])  # type: ignore[index]
+        x[(..., *i)] = glass.core.algorithm.nnls(r, y[i])  # type: ignore[index]
 
     # all done
     return x
diff --git a/pyproject.toml b/pyproject.toml
index cf93baf6..74f5b008 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -173,6 +173,8 @@ lint.per-file-ignores = {"__init__.py" = [
     "EM102", # f-string-in-exception
     "PLR2004", # TODO: magic-value-comparison
     "TRY003", # raise-vanilla-args
+], "glass/grf/_*" = [
+    "SLF001", # private-member-access
 ], "noxfile.py" = [
     "T201", # print
 ], "tests*" = [
diff --git a/tests/conftest.py b/tests/conftest.py
index bbb8010e..a9fd0851 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -10,7 +10,7 @@
 
 from cosmology import Cosmology
 
-from glass import RadialWindow
+import glass
 
 # Handling of array backends, inspired by-
 # https://github.com/scipy/scipy/blob/36e349b6afbea057cb713fc314296f10d55194cc/scipy/conftest.py#L139
@@ -156,11 +156,11 @@ def rng() -> np.random.Generator:
 
 
 @pytest.fixture(scope="session")
-def shells() -> list[RadialWindow]:
+def shells() -> list[glass.RadialWindow]:
     return [
-        RadialWindow(np.array([0.0, 1.0, 2.0]), np.array([0.0, 1.0, 0.0]), 1.0),
-        RadialWindow(np.array([1.0, 2.0, 3.0]), np.array([0.0, 1.0, 0.0]), 2.0),
-        RadialWindow(np.array([2.0, 3.0, 4.0]), np.array([0.0, 1.0, 0.0]), 3.0),
-        RadialWindow(np.array([3.0, 4.0, 5.0]), np.array([0.0, 1.0, 0.0]), 4.0),
-        RadialWindow(np.array([4.0, 5.0, 6.0]), np.array([0.0, 1.0, 0.0]), 5.0),
+        glass.RadialWindow(np.array([0.0, 1.0, 2.0]), np.array([0.0, 1.0, 0.0]), 1.0),
+        glass.RadialWindow(np.array([1.0, 2.0, 3.0]), np.array([0.0, 1.0, 0.0]), 2.0),
+        glass.RadialWindow(np.array([2.0, 3.0, 4.0]), np.array([0.0, 1.0, 0.0]), 3.0),
+        glass.RadialWindow(np.array([3.0, 4.0, 5.0]), np.array([0.0, 1.0, 0.0]), 4.0),
+        glass.RadialWindow(np.array([4.0, 5.0, 6.0]), np.array([0.0, 1.0, 0.0]), 5.0),
     ]
diff --git a/tests/core/test_algorithm.py b/tests/core/test_algorithm.py
index 8e6dcc68..8f0541e6 100644
--- a/tests/core/test_algorithm.py
+++ b/tests/core/test_algorithm.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from glass.core.algorithm import nnls
+import glass.core.algorithm
 
 
 def test_nnls(rng: np.random.Generator) -> None:
@@ -10,14 +10,16 @@ def test_nnls(rng: np.random.Generator) -> None:
     a = np.arange(25.0).reshape(-1, 5)
     b = np.arange(5.0)
     y = a @ b
-    res = nnls(a, y)
+    res = glass.core.algorithm.nnls(a, y)
     assert np.linalg.norm((a @ res) - y) < 1e-7
 
     a = rng.uniform(low=-10, high=10, size=[50, 10])
     b = np.abs(rng.uniform(low=-2, high=2, size=[10]))
     b[::2] = 0
     x = a @ b
-    res = nnls(a, x, tol=500 * np.linalg.norm(a, 1) * np.spacing(1.0))
+    res = glass.core.algorithm.nnls(
+        a, x, tol=500 * np.linalg.norm(a, 1) * np.spacing(1.0)
+    )
     np.testing.assert_allclose(res, b, rtol=0.0, atol=1e-10)
 
     # check matrix and vector's shape
@@ -26,8 +28,8 @@ def test_nnls(rng: np.random.Generator) -> None:
     b = rng.standard_normal((100,))
 
     with pytest.raises(ValueError, match="input `a` is not a matrix"):
-        nnls(b, a)
+        glass.core.algorithm.nnls(b, a)
     with pytest.raises(ValueError, match="input `b` is not a vector"):
-        nnls(a, a)
+        glass.core.algorithm.nnls(a, a)
     with pytest.raises(ValueError, match="the shapes of `a` and `b` do not match"):
-        nnls(a.T, b)
+        glass.core.algorithm.nnls(a.T, b)
diff --git a/tests/core/test_array.py b/tests/core/test_array.py
index 10fbcd95..55a0c916 100644
--- a/tests/core/test_array.py
+++ b/tests/core/test_array.py
@@ -3,13 +3,7 @@
 import numpy as np
 import pytest
 
-from glass.core.array import (
-    broadcast_first,
-    broadcast_leading_axes,
-    cumulative_trapezoid,
-    ndinterp,
-    trapezoid_product,
-)
+import glass.core.array
 
 if TYPE_CHECKING:
     from numpy.typing import NDArray
@@ -21,7 +15,7 @@ def test_broadcast_first() -> None:
 
     # arrays with shape ((3, 4, 2)) and ((1, 2)) are passed
     # to np.broadcast_arrays; hence it works
-    a_a, b_a = broadcast_first(a, b)
+    a_a, b_a = glass.core.array.broadcast_first(a, b)
     assert a_a.shape == (2, 3, 4)
     assert b_a.shape == (2, 3, 4)
 
@@ -35,7 +29,7 @@ def test_broadcast_first() -> None:
     b = np.ones((5, 6))
 
     with pytest.raises(ValueError, match="shape mismatch"):
-        broadcast_first(a, b)
+        glass.core.array.broadcast_first(a, b)
 
     # plain np.broadcast_arrays will work
     a_a, b_a = np.broadcast_arrays(a, b)
@@ -49,7 +43,9 @@ def test_broadcast_leading_axes() -> None:
     b_in = np.zeros((4, 10))
     c_in = np.zeros((3, 1, 5, 6))
 
-    dims, *rest = broadcast_leading_axes((a_in, 0), (b_in, 1), (c_in, 2))
+    dims, *rest = glass.core.array.broadcast_leading_axes(
+        (a_in, 0), (b_in, 1), (c_in, 2)
+    )
     a_out, b_out, c_out = rest
 
     assert dims == (3, 4)
@@ -65,17 +61,17 @@ def test_ndinterp() -> None:
     yp = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
 
     x: float | NDArray[np.float64] = 0.5
-    y = ndinterp(x, xp, yp)
+    y = glass.core.array.ndinterp(x, xp, yp)
     assert np.shape(y) == ()
     np.testing.assert_allclose(y, 1.15, atol=1e-15)
 
     x = np.array([0.5, 1.5, 2.5])
-    y = ndinterp(x, xp, yp)
+    y = glass.core.array.ndinterp(x, xp, yp)
     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 = ndinterp(x, xp, yp)
+    y = glass.core.array.ndinterp(x, xp, yp)
     assert np.shape(y) == (2, 2)
     np.testing.assert_allclose(y, [[1.15, 1.25], [1.35, 1.45]], atol=1e-15)
 
@@ -84,17 +80,17 @@ def test_ndinterp() -> None:
     yp = 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 = ndinterp(x, xp, yp)
+    y = glass.core.array.ndinterp(x, xp, yp)
     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 = ndinterp(x, xp, yp)
+    y = glass.core.array.ndinterp(x, xp, yp)
     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 = ndinterp(x, xp, yp)
+    y = glass.core.array.ndinterp(x, xp, yp)
     assert np.shape(y) == (2, 2, 2)
     np.testing.assert_allclose(
         y,
@@ -109,12 +105,12 @@ def test_ndinterp() -> None:
     )
 
     x = 0.5
-    y = ndinterp(x, xp, yp, axis=1)
+    y = glass.core.array.ndinterp(x, xp, yp, 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 = ndinterp(x, xp, yp, axis=1)
+    y = glass.core.array.ndinterp(x, xp, yp, axis=1)
     assert np.shape(y) == (2, 3, 1)
     np.testing.assert_allclose(
         y,
@@ -123,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 = ndinterp(x, xp, yp, axis=1)
+    y = glass.core.array.ndinterp(x, xp, yp, axis=1)
     assert np.shape(y) == (2, 3, 4, 1)
     np.testing.assert_allclose(
         y,
@@ -150,7 +146,7 @@ def test_trapezoid_product() -> None:
     x2 = np.linspace(1, 2, 10)
     f2 = np.full_like(x2, 0.5)
 
-    s = trapezoid_product((x1, f1), (x2, f2))
+    s = glass.core.array.trapezoid_product((x1, f1), (x2, f2))
 
     np.testing.assert_allclose(s, 1.0)
 
@@ -163,18 +159,18 @@ def test_cumulative_trapezoid() -> None:
 
     # default dtype (int)
 
-    ct = cumulative_trapezoid(f, x)
+    ct = glass.core.array.cumulative_trapezoid(f, x)
     np.testing.assert_allclose(ct, np.array([0, 1, 4, 7]))
 
     # explicit dtype (float)
 
-    ct = cumulative_trapezoid(f, x, dtype=float)
+    ct = glass.core.array.cumulative_trapezoid(f, x, dtype=float)
     np.testing.assert_allclose(ct, np.array([0.0, 1.5, 4.0, 7.5]))
 
     # explicit return array
 
     out = np.zeros((4,))
-    ct = cumulative_trapezoid(f, x, dtype=float, out=out)
+    ct = glass.core.array.cumulative_trapezoid(f, x, dtype=float, out=out)
     np.testing.assert_equal(ct, out)
 
     # 2D f and 1D x
@@ -184,12 +180,12 @@ def test_cumulative_trapezoid() -> None:
 
     # default dtype (int)
 
-    ct = cumulative_trapezoid(f, x)
+    ct = glass.core.array.cumulative_trapezoid(f, x)
     np.testing.assert_allclose(ct, np.array([[0, 2, 12, 31], [0, 2, 8, 17]]))
 
     # explicit dtype (float)
 
-    ct = cumulative_trapezoid(f, x, dtype=float)
+    ct = glass.core.array.cumulative_trapezoid(f, x, dtype=float)
     np.testing.assert_allclose(
         ct, np.array([[0.0, 2.5, 12.25, 31.0], [0.0, 2.5, 8.5, 17.5]])
     )
@@ -197,5 +193,5 @@ def test_cumulative_trapezoid() -> None:
     # explicit return array
 
     out = np.zeros((2, 4))
-    ct = cumulative_trapezoid(f, x, dtype=float, out=out)
+    ct = glass.core.array.cumulative_trapezoid(f, x, dtype=float, out=out)
     np.testing.assert_equal(ct, out)
diff --git a/tests/test_fields.py b/tests/test_fields.py
index ed05f337..051168f7 100644
--- a/tests/test_fields.py
+++ b/tests/test_fields.py
@@ -3,19 +3,6 @@
 import pytest
 
 import glass
-from glass import (
-    cls2cov,
-    discretized_cls,
-    effective_cls,
-    generate_gaussian,
-    generate_lognormal,
-    getcl,
-    iternorm,
-    lognormal_gls,
-    multalm,
-    transform_cls,
-)
-from glass.fields import inv_triangle_number
 
 
 def test_iternorm() -> None:
@@ -23,7 +10,9 @@ def test_iternorm() -> None:
 
     k = 2
 
-    generator = iternorm(k, np.array([1.0, 0.5, 0.5, 0.5, 0.2, 0.1, 0.5, 0.1, 0.2]))
+    generator = glass.iternorm(
+        k, np.array([1.0, 0.5, 0.5, 0.5, 0.2, 0.1, 0.5, 0.1, 0.2])
+    )
     result = next(generator)
 
     j, a, s = result
@@ -37,7 +26,7 @@ def test_iternorm() -> None:
 
     size = 3
 
-    generator = iternorm(
+    generator = glass.iternorm(
         k,
         np.array(
             [
@@ -59,18 +48,20 @@ def test_iternorm() -> None:
     # test shape mismatch error
 
     with pytest.raises(TypeError, match="covariance row 0: shape"):
-        list(iternorm(k, np.array([[1.0, 0.5], [0.5, 0.2]])))
+        list(glass.iternorm(k, np.array([[1.0, 0.5], [0.5, 0.2]])))
 
     # test positive definite error
 
     with pytest.raises(ValueError, match="covariance matrix is not positive definite"):
-        list(iternorm(k, np.array([1.0, 0.5, 0.9, 0.5, 0.2, 0.4, 0.9, 0.4, -1.0])))
+        list(
+            glass.iternorm(k, np.array([1.0, 0.5, 0.9, 0.5, 0.2, 0.4, 0.9, 0.4, -1.0]))
+        )
 
     # test multiple iterations
 
     size = (3,)
 
-    generator = iternorm(
+    generator = glass.iternorm(
         k,
         np.array(
             [
@@ -92,7 +83,7 @@ def test_iternorm() -> None:
 
     # test k = 0
 
-    generator = iternorm(0, np.array([1.0]))
+    generator = glass.iternorm(0, np.array([1.0]))
 
     j, a, s = result
 
@@ -107,7 +98,7 @@ def test_cls2cov() -> None:
 
     nl, nf, nc = 3, 2, 2
 
-    generator = cls2cov(
+    generator = glass.cls2cov(
         [np.array([1.0, 0.5, 0.3]), None, np.array([0.7, 0.6, 0.1])],  # type: ignore[list-item]
         nl,
         nf,
@@ -123,7 +114,7 @@ def test_cls2cov() -> None:
 
     # test negative value error
 
-    generator = cls2cov(
+    generator = glass.cls2cov(
         np.array(  # type: ignore[arg-type]
             [
                 [-1.0, 0.5, 0.3],
@@ -142,7 +133,7 @@ def test_cls2cov() -> None:
 
     nl, nf, nc = 3, 3, 2
 
-    generator = cls2cov(
+    generator = glass.cls2cov(
         np.array(  # type: ignore[arg-type]
             [
                 [1.0, 0.5, 0.3],
@@ -177,7 +168,7 @@ def test_multalm() -> None:
     bl = np.array([2.0, 0.5, 1.0])
     alm_copy = np.copy(alm)
 
-    result = multalm(alm, bl, inplace=True)
+    result = glass.multalm(alm, bl, inplace=True)
 
     assert np.array_equal(result, alm)  # in-place
     expected_result = np.array([2.0, 1.0, 3.0, 2.0, 5.0, 6.0])
@@ -188,14 +179,14 @@ def test_multalm() -> None:
 
     bl = np.ones(3)
 
-    result = multalm(alm, bl, inplace=False)
+    result = glass.multalm(alm, bl, inplace=False)
     np.testing.assert_array_equal(result, alm)
 
     # multiple with 0s
 
     bl = np.array([0.0, 1.0, 0.0])
 
-    result = multalm(alm, bl, inplace=False)
+    result = glass.multalm(alm, bl, inplace=False)
 
     expected_result = np.array([0.0, 1.0, 0.0, 2.0, 0.0, 0.0])
     np.testing.assert_allclose(result, expected_result)
@@ -205,7 +196,7 @@ def test_multalm() -> None:
     alm = np.array([])
     bl = np.array([])
 
-    result = multalm(alm, bl, inplace=False)
+    result = glass.multalm(alm, bl, inplace=False)
     np.testing.assert_array_equal(result, alm)
 
 
@@ -216,24 +207,24 @@ def test_transform_cls() -> None:
 
     # empty cls
 
-    assert transform_cls([], tfm, pars) == []
+    assert glass.transform_cls([], tfm, pars) == []
 
     # check output shape
 
-    assert len(transform_cls([sub_cls], tfm, pars)) == 1
-    assert len(transform_cls([sub_cls], tfm, pars)[0]) == 5
+    assert len(glass.transform_cls([sub_cls], tfm, pars)) == 1
+    assert len(glass.transform_cls([sub_cls], tfm, pars)[0]) == 5
 
-    assert len(transform_cls([sub_cls, sub_cls], tfm, pars)) == 2
-    assert len(transform_cls([sub_cls, sub_cls], tfm, pars)[0]) == 5
-    assert len(transform_cls([sub_cls, sub_cls], tfm, pars)[1]) == 5
+    assert len(glass.transform_cls([sub_cls, sub_cls], tfm, pars)) == 2
+    assert len(glass.transform_cls([sub_cls, sub_cls], tfm, pars)[0]) == 5
+    assert len(glass.transform_cls([sub_cls, sub_cls], tfm, pars)[1]) == 5
 
     # one sequence of empty cls
 
-    assert transform_cls([[], sub_cls], tfm, pars)[0] == []
+    assert glass.transform_cls([[], sub_cls], tfm, pars)[0] == []
 
     # monopole behavior
 
-    assert transform_cls([sub_cls, np.linspace(0, 5, 6)], tfm, pars)[1][0] == 0
+    assert glass.transform_cls([sub_cls, np.linspace(0, 5, 6)], tfm, pars)[1][0] == 0
 
 
 def test_lognormal_gls() -> None:
@@ -241,31 +232,38 @@ def test_lognormal_gls() -> None:
 
     # empty cls
 
-    assert lognormal_gls([], shift) == []
+    assert glass.lognormal_gls([], shift) == []
 
     # check output shape
 
-    assert len(lognormal_gls([np.linspace(1, 5, 5)], shift)) == 1
-    assert len(lognormal_gls([np.linspace(1, 5, 5)], shift)[0]) == 5
+    assert len(glass.lognormal_gls([np.linspace(1, 5, 5)], shift)) == 1
+    assert len(glass.lognormal_gls([np.linspace(1, 5, 5)], shift)[0]) == 5
 
-    assert len(lognormal_gls([np.linspace(1, 5, 5), np.linspace(1, 5, 5)], shift)) == 2
     assert (
-        len(lognormal_gls([np.linspace(1, 5, 5), np.linspace(1, 5, 5)], shift)[0]) == 5
+        len(glass.lognormal_gls([np.linspace(1, 5, 5), np.linspace(1, 5, 5)], shift))
+        == 2
+    )
+    assert (
+        len(glass.lognormal_gls([np.linspace(1, 5, 5), np.linspace(1, 5, 5)], shift)[0])
+        == 5
     )
     assert (
-        len(lognormal_gls([np.linspace(1, 5, 5), np.linspace(1, 5, 5)], shift)[1]) == 5
+        len(glass.lognormal_gls([np.linspace(1, 5, 5), np.linspace(1, 5, 5)], shift)[1])
+        == 5
     )
 
 
 def test_discretized_cls() -> None:
     # empty cls
 
-    result = discretized_cls([])
+    result = glass.discretized_cls([])
     assert result == []
 
     # power spectra truncated at lmax + 1 if lmax provided
 
-    result = discretized_cls([np.arange(10), np.arange(10), np.arange(10)], lmax=5)
+    result = glass.discretized_cls(
+        [np.arange(10), np.arange(10), np.arange(10)], lmax=5
+    )
 
     for cl in result:
         assert len(cl) == 6
@@ -275,13 +273,13 @@ def test_discretized_cls() -> None:
     with pytest.raises(
         ValueError, match="length of cls array is not a triangle number"
     ):
-        discretized_cls([np.arange(10), np.arange(10)], ncorr=1)
+        glass.discretized_cls([np.arange(10), np.arange(10)], ncorr=1)
 
     # ncorr not None
 
     cls = [np.arange(10), np.arange(10), np.arange(10)]
     ncorr = 0
-    result = discretized_cls(cls, ncorr=ncorr)
+    result = glass.discretized_cls(cls, ncorr=ncorr)
 
     assert len(result[0]) == 10
     assert len(result[1]) == 10
@@ -293,7 +291,7 @@ def test_discretized_cls() -> None:
 
     pw = hp.pixwin(nside, lmax=7)
 
-    result = discretized_cls([[], np.ones(10), np.ones(10)], nside=nside)
+    result = glass.discretized_cls([[], np.ones(10), np.ones(10)], nside=nside)
 
     for cl in result:
         n = min(len(cl), len(pw))
@@ -304,37 +302,37 @@ def test_discretized_cls() -> None:
 def test_effective_cls() -> None:
     # empty cls
 
-    result = effective_cls([], np.array([]))
+    result = glass.effective_cls([], np.array([]))
     assert result.shape == (0,)
 
     # check ValueError for triangle number
 
     with pytest.raises(ValueError, match="length of cls is not a triangle number"):
-        effective_cls([np.arange(10), np.arange(10)], np.ones((2, 1)))
+        glass.effective_cls([np.arange(10), np.arange(10)], np.ones((2, 1)))
 
     # check ValueError for triangle number
 
     with pytest.raises(ValueError, match="shape mismatch between fields and weights1"):
-        effective_cls([], np.ones((3, 1)))
+        glass.effective_cls([], np.ones((3, 1)))
 
     # check with only weights1
 
     cls = [np.arange(15), np.arange(15), np.arange(15)]
     weights1 = np.ones((2, 1))
 
-    result = effective_cls(cls, weights1)
+    result = glass.effective_cls(cls, weights1)
     assert result.shape == (1, 1, 15)
 
     # check truncation if lmax provided
 
-    result = effective_cls(cls, weights1, lmax=5)
+    result = glass.effective_cls(cls, weights1, lmax=5)
 
     assert result.shape == (1, 1, 6)
     np.testing.assert_allclose(result[..., 6:], 0)
 
     # check with weights1 and weights2 and weights1 is weights2
 
-    result = effective_cls(cls, weights1, weights2=weights1)
+    result = glass.effective_cls(cls, weights1, weights2=weights1)
     assert result.shape == (1, 1, 15)
 
 
@@ -343,26 +341,28 @@ def test_generate_gaussian() -> None:
     nside = 4
     ncorr = 1
 
-    gaussian_fields = list(generate_gaussian(gls, nside))
+    gaussian_fields = list(glass.generate_gaussian(gls, nside))
 
     assert gaussian_fields[0].shape == (hp.nside2npix(nside),)
 
     # requires resetting the RNG for reproducibility
     rng = np.random.default_rng(seed=42)
-    gaussian_fields = list(generate_gaussian(gls, nside, rng=rng))
+    gaussian_fields = list(glass.generate_gaussian(gls, nside, rng=rng))
 
     assert gaussian_fields[0].shape == (hp.nside2npix(nside),)
 
     # requires resetting the RNG for reproducibility
     rng = np.random.default_rng(seed=42)
-    new_gaussian_fields = list(generate_gaussian(gls, nside, ncorr=ncorr, rng=rng))
+    new_gaussian_fields = list(
+        glass.generate_gaussian(gls, nside, ncorr=ncorr, rng=rng)
+    )
 
     assert new_gaussian_fields[0].shape == (hp.nside2npix(nside),)
 
     np.testing.assert_allclose(new_gaussian_fields[0], gaussian_fields[0])
 
     with pytest.raises(ValueError, match="all gls are empty"):
-        list(generate_gaussian([], nside))
+        list(glass.generate_gaussian([], nside))
 
 
 def test_generate_lognormal(rng: np.random.Generator) -> None:
@@ -371,7 +371,7 @@ def test_generate_lognormal(rng: np.random.Generator) -> None:
 
     # check shape and values
 
-    lognormal_fields = list(generate_lognormal(gls, nside, shift=1.0, rng=rng))
+    lognormal_fields = list(glass.generate_lognormal(gls, nside, shift=1.0, rng=rng))
 
     assert len(lognormal_fields) == len(gls)
     for m in lognormal_fields:
@@ -382,10 +382,12 @@ def test_generate_lognormal(rng: np.random.Generator) -> None:
 
     # requires resetting the RNG to obtain exact the same result multiplied by 2 (shift)
     rng = np.random.default_rng(seed=42)
-    lognormal_fields = list(generate_lognormal(gls, nside, shift=1.0, rng=rng))
+    lognormal_fields = list(glass.generate_lognormal(gls, nside, shift=1.0, rng=rng))
 
     rng = np.random.default_rng(seed=42)
-    new_lognormal_fields = list(generate_lognormal(gls, nside, shift=2.0, rng=rng))
+    new_lognormal_fields = list(
+        glass.generate_lognormal(gls, nside, shift=2.0, rng=rng)
+    )
 
     for ms, mu in zip(new_lognormal_fields, lognormal_fields, strict=False):
         np.testing.assert_allclose(ms, mu * 2.0)
@@ -399,18 +401,18 @@ def test_getcl() -> None:
     # make sure indices are retrieved correctly
     for i in range(10):
         for j in range(10):
-            result = getcl(cls, i, j)
+            result = glass.getcl(cls, i, j)
             expected = np.array([min(i, j), max(i, j)], dtype=np.float64)
             np.testing.assert_array_equal(np.sort(result), expected)
 
             # check slicing
-            result = getcl(cls, i, j, lmax=0)
+            result = glass.getcl(cls, i, j, lmax=0)
             expected = np.array([max(i, j)], dtype=np.float64)
             assert len(result) == 1
             np.testing.assert_array_equal(result, expected)
 
             # check padding
-            result = getcl(cls, i, j, lmax=50)
+            result = glass.getcl(cls, i, j, lmax=50)
             expected = np.zeros((49,), dtype=np.float64)
             assert len(result) == 51
             np.testing.assert_array_equal(result[2:], expected)
@@ -418,13 +420,13 @@ def test_getcl() -> None:
 
 def test_inv_triangle_number():
     for n in range(10_000):
-        assert inv_triangle_number(n * (n + 1) // 2) == n
+        assert glass.fields.inv_triangle_number(n * (n + 1) // 2) == n
 
     not_triangle_numbers = [2, 4, 5, 7, 8, 9, 11, 12, 13, 14, 16, 17, 18, 19, 20]
 
     for t in not_triangle_numbers:
         with pytest.raises(ValueError, match="not a triangle number"):
-            inv_triangle_number(t)
+            glass.fields.inv_triangle_number(t)
 
 
 def test_enumerate_spectra():
diff --git a/tests/test_galaxies.py b/tests/test_galaxies.py
index 5c5185bd..b746bf9c 100644
--- a/tests/test_galaxies.py
+++ b/tests/test_galaxies.py
@@ -4,12 +4,7 @@
 import pytest
 import pytest_mock
 
-from glass import (
-    galaxy_shear,
-    gaussian_phz,
-    redshifts,
-    redshifts_from_nz,
-)
+import glass
 
 if TYPE_CHECKING:
     from numpy.typing import NDArray
@@ -22,20 +17,20 @@ def test_redshifts(mocker: pytest_mock.MockerFixture) -> None:
     w.wa = np.exp(-0.5 * (w.za - 0.5) ** 2 / 0.1**2)
 
     # sample redshifts (scalar)
-    z = redshifts(13, w)
+    z = glass.redshifts(13, w)
     assert z.shape == (13,)
     assert z.min() >= 0.0
     assert z.max() <= 1.0
 
     # sample redshifts (array)
-    z = redshifts(np.array([[1, 2], [3, 4]]), w)
+    z = glass.redshifts(np.array([[1, 2], [3, 4]]), w)
     assert z.shape == (10,)
 
 
 def test_redshifts_from_nz(rng: np.random.Generator) -> None:
     # test sampling
 
-    redshifts = redshifts_from_nz(
+    redshifts = glass.redshifts_from_nz(
         10,
         np.array([0, 1, 2, 3, 4]),
         np.array([1, 0, 0, 0, 0]),
@@ -43,7 +38,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
     )
     assert np.all((0 <= redshifts) & (redshifts <= 1))  # noqa: SIM300
 
-    redshifts = redshifts_from_nz(
+    redshifts = glass.redshifts_from_nz(
         10,
         np.array([0, 1, 2, 3, 4]),
         np.array([0, 0, 1, 0, 0]),
@@ -51,7 +46,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
     )
     assert np.all((1 <= redshifts) & (redshifts <= 3))  # noqa: SIM300
 
-    redshifts = redshifts_from_nz(
+    redshifts = glass.redshifts_from_nz(
         10,
         np.array([0, 1, 2, 3, 4]),
         np.array([0, 0, 0, 0, 1]),
@@ -59,7 +54,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
     )
     assert np.all((3 <= redshifts) & (redshifts <= 4))  # noqa: SIM300
 
-    redshifts = redshifts_from_nz(
+    redshifts = glass.redshifts_from_nz(
         10,
         np.array([0, 1, 2, 3, 4]),
         np.array([0, 0, 1, 1, 1]),
@@ -69,7 +64,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
 
     # test with rng
 
-    redshifts = redshifts_from_nz(
+    redshifts = glass.redshifts_from_nz(
         10,
         np.array([0, 1, 2, 3, 4]),
         np.array([0, 0, 1, 1, 1]),
@@ -86,7 +81,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
     z = np.linspace(0, 1, 100)
     nz = z * (1 - z)
 
-    redshifts = redshifts_from_nz(count, z, nz, warn=False)
+    redshifts = glass.redshifts_from_nz(count, z, nz, warn=False)
 
     assert redshifts.shape == (count,)
     assert np.all((0 <= redshifts) & (redshifts <= 1))  # noqa: SIM300
@@ -97,7 +92,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
     z = np.linspace(0, 1, 100)
     nz = z * (1 - z)
 
-    redshifts = redshifts_from_nz(count, z, nz, warn=False)
+    redshifts = glass.redshifts_from_nz(count, z, nz, warn=False)
 
     assert np.shape(redshifts) == (60,)
 
@@ -107,7 +102,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
     z = np.linspace(0, 1, 100)
     nz = np.array([z * (1 - z), (z - 0.5) ** 2])
 
-    redshifts = redshifts_from_nz(count, z, nz, warn=False)
+    redshifts = glass.redshifts_from_nz(count, z, nz, warn=False)
 
     assert redshifts.shape == (20,)
 
@@ -117,7 +112,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
     z = np.linspace(0, 1, 100)
     nz = np.array([z * (1 - z), (z - 0.5) ** 2])
 
-    redshifts = redshifts_from_nz(count, z, nz, warn=False)
+    redshifts = glass.redshifts_from_nz(count, z, nz, warn=False)
 
     assert redshifts.shape == (120,)
 
@@ -128,10 +123,10 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
     nz = np.array([z * (1 - z), (z - 0.5) ** 2])
 
     with pytest.raises(ValueError, match="shape mismatch"):
-        redshifts_from_nz(count, z, nz, warn=False)
+        glass.redshifts_from_nz(count, z, nz, warn=False)
 
     with pytest.warns(UserWarning, match="when sampling galaxies"):
-        redshifts = redshifts_from_nz(
+        redshifts = glass.redshifts_from_nz(
             10,
             np.array([0, 1, 2, 3, 4]),
             np.array([1, 0, 0, 0, 0]),
@@ -147,7 +142,7 @@ def test_galaxy_shear(rng: np.random.Generator) -> None:
         rng.normal(size=(12,)),
     )
 
-    shear = galaxy_shear(
+    shear = glass.galaxy_shear(
         np.array([]),
         np.array([]),
         np.array([]),
@@ -162,12 +157,12 @@ def test_galaxy_shear(rng: np.random.Generator) -> None:
         rng.normal(size=(512,)),
         rng.normal(size=(512,)),
     )
-    shear = galaxy_shear(gal_lon, gal_lat, gal_eps, kappa, gamma1, gamma2)
+    shear = glass.galaxy_shear(gal_lon, gal_lat, gal_eps, kappa, gamma1, gamma2)
     assert np.shape(shear) == (512,)
 
     # shape with no reduced shear
 
-    shear = galaxy_shear(
+    shear = glass.galaxy_shear(
         np.array([]),
         np.array([]),
         np.array([]),
@@ -183,7 +178,7 @@ def test_galaxy_shear(rng: np.random.Generator) -> None:
         rng.normal(size=(512,)),
         rng.normal(size=(512,)),
     )
-    shear = galaxy_shear(
+    shear = glass.galaxy_shear(
         gal_lon,
         gal_lat,
         gal_eps,
@@ -203,13 +198,13 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
     z: float | NDArray[np.float64] = np.linspace(0, 1, 100)
     sigma_0: float | NDArray[np.float64] = 0.0
 
-    phz = gaussian_phz(z, sigma_0)
+    phz = glass.gaussian_phz(z, sigma_0)
 
     np.testing.assert_array_equal(z, phz)
 
     # test with rng
 
-    phz = gaussian_phz(z, sigma_0, rng=rng)
+    phz = glass.gaussian_phz(z, sigma_0, rng=rng)
 
     np.testing.assert_array_equal(z, phz)
 
@@ -218,7 +213,7 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
     z = 0.0
     sigma_0 = np.ones(100)
 
-    phz = gaussian_phz(z, sigma_0)
+    phz = glass.gaussian_phz(z, sigma_0)
 
     assert isinstance(phz, np.ndarray)
     assert phz.shape == (100,)
@@ -229,7 +224,7 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
     z = 1.0
     sigma_0 = np.ones(100)
 
-    phz = gaussian_phz(z, sigma_0, lower=0.5, upper=1.5)
+    phz = glass.gaussian_phz(z, sigma_0, lower=0.5, upper=1.5)
 
     assert isinstance(phz, np.ndarray)
     assert phz.shape == (100,)
@@ -243,7 +238,7 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
     z = 1.0
     sigma_0 = 0.0
 
-    phz = gaussian_phz(z, sigma_0)
+    phz = glass.gaussian_phz(z, sigma_0)
 
     assert np.ndim(phz) == 0
     assert phz == z
@@ -253,7 +248,7 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
     z = np.linspace(0, 1, 10)
     sigma_0 = 0.0
 
-    phz = gaussian_phz(z, sigma_0)
+    phz = glass.gaussian_phz(z, sigma_0)
 
     assert isinstance(phz, np.ndarray)
     assert phz.shape == (10,)
@@ -264,7 +259,7 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
     z = 1.0
     sigma_0 = np.zeros(10)
 
-    phz = gaussian_phz(z, sigma_0)
+    phz = glass.gaussian_phz(z, sigma_0)
 
     assert isinstance(phz, np.ndarray)
     assert phz.shape == (10,)
@@ -275,7 +270,7 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
     z = np.linspace(0, 1, 10)
     sigma_0 = np.zeros((11, 1))
 
-    phz = gaussian_phz(z, sigma_0)
+    phz = glass.gaussian_phz(z, sigma_0)
 
     assert isinstance(phz, np.ndarray)
     assert phz.shape == (11, 10)
@@ -283,13 +278,13 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
 
     # test resampling
 
-    phz = gaussian_phz(np.array(0.0), np.array(1.0), lower=np.array([0]))
+    phz = glass.gaussian_phz(np.array(0.0), np.array(1.0), lower=np.array([0]))
     assert isinstance(phz, float)
 
-    phz = gaussian_phz(np.array(0.0), np.array(1.0), upper=np.array([1]))
+    phz = glass.gaussian_phz(np.array(0.0), np.array(1.0), upper=np.array([1]))
     assert isinstance(phz, float)
 
     # test error
 
     with pytest.raises(ValueError, match="requires lower < upper"):
-        phz = gaussian_phz(z, sigma_0, lower=np.array([1]), upper=np.array([0]))
+        phz = glass.gaussian_phz(z, sigma_0, lower=np.array([1]), upper=np.array([0]))
diff --git a/tests/test_lensing.py b/tests/test_lensing.py
index 93c9644f..10d2a711 100644
--- a/tests/test_lensing.py
+++ b/tests/test_lensing.py
@@ -6,22 +6,14 @@
 import numpy as np
 import pytest
 
-from glass import (
-    MultiPlaneConvergence,
-    RadialWindow,
-    deflect,
-    from_convergence,
-    multi_plane_matrix,
-    multi_plane_weights,
-    shear_from_convergence,  # noqa: F401
-)
+import glass
 
 if TYPE_CHECKING:
     from cosmology import Cosmology
 
 
 def test_from_convergence(rng: np.random.Generator) -> None:
-    """Add unit tests for :func:`from_convergence`."""
+    """Add unit tests for :func:`glass.from_convergence`."""
     # l_max = 100  # noqa: ERA001
     n_side = 32
 
@@ -30,48 +22,48 @@ def test_from_convergence(rng: np.random.Generator) -> None:
 
     # check with all False
 
-    results = from_convergence(kappa)
+    results = glass.from_convergence(kappa)
     np.testing.assert_array_equal(results, ())
 
     # check all combinations of potential, deflection, shear being True
 
-    results = from_convergence(kappa, potential=True)
+    results = glass.from_convergence(kappa, potential=True)
     np.testing.assert_array_equal(len(results), 1)
 
-    results = from_convergence(kappa, deflection=True)
+    results = glass.from_convergence(kappa, deflection=True)
     np.testing.assert_array_equal(len(results), 1)
 
-    results = from_convergence(kappa, shear=True)
+    results = glass.from_convergence(kappa, shear=True)
     np.testing.assert_array_equal(len(results), 1)
 
-    results = from_convergence(kappa, potential=True, deflection=True)
+    results = glass.from_convergence(kappa, potential=True, deflection=True)
     np.testing.assert_array_equal(len(results), 2)
 
-    results = from_convergence(kappa, potential=True, shear=True)
+    results = glass.from_convergence(kappa, potential=True, shear=True)
     np.testing.assert_array_equal(len(results), 2)
 
-    results = from_convergence(kappa, deflection=True, shear=True)
+    results = glass.from_convergence(kappa, deflection=True, shear=True)
     np.testing.assert_array_equal(len(results), 2)
 
-    results = from_convergence(kappa, potential=True, deflection=True, shear=True)
+    results = glass.from_convergence(kappa, potential=True, deflection=True, shear=True)
     np.testing.assert_array_equal(len(results), 3)
 
 
 def test_shear_from_convergence() -> None:
-    """Add unit tests for :func:`shear_from_convergence`."""
+    """Add unit tests for :func:`glass.shear_from_convergence`."""
 
 
 def test_multi_plane_matrix(
-    shells: list[RadialWindow],
+    shells: list[glass.RadialWindow],
     cosmo: Cosmology,
     rng: np.random.Generator,
 ) -> None:
-    mat = multi_plane_matrix(shells, cosmo)
+    mat = glass.multi_plane_matrix(shells, cosmo)
 
     np.testing.assert_array_equal(mat, np.tril(mat))
     np.testing.assert_array_equal(np.triu(mat, 1), 0)
 
-    convergence = MultiPlaneConvergence(cosmo)
+    convergence = glass.MultiPlaneConvergence(cosmo)
 
     deltas = rng.random((len(shells), 10))
     kappas = []
@@ -84,17 +76,17 @@ def test_multi_plane_matrix(
 
 
 def test_multi_plane_weights(
-    shells: list[RadialWindow],
+    shells: list[glass.RadialWindow],
     cosmo: Cosmology,
     rng: np.random.Generator,
 ) -> None:
     w_in = np.eye(len(shells))
-    w_out = multi_plane_weights(w_in, shells, cosmo)
+    w_out = glass.multi_plane_weights(w_in, shells, cosmo)
 
     np.testing.assert_array_equal(w_out, np.triu(w_out, 1))
     np.testing.assert_array_equal(np.tril(w_out), 0)
 
-    convergence = MultiPlaneConvergence(cosmo)
+    convergence = glass.MultiPlaneConvergence(cosmo)
 
     deltas = rng.random((len(shells), 10))
     weights = rng.random((len(shells), 3))
@@ -104,7 +96,7 @@ def test_multi_plane_weights(
         kappa = kappa + weight[..., None] * convergence.kappa
     kappa /= weights.sum(axis=0)[..., None]
 
-    wmat = multi_plane_weights(weights, shells, cosmo)
+    wmat = glass.multi_plane_weights(weights, shells, cosmo)
 
     np.testing.assert_allclose(np.einsum("ij,ik", wmat, deltas), kappa)
 
@@ -118,19 +110,19 @@ def alpha(re: float, im: float, *, usecomplex: bool) -> complex | list[float]:
         return re + 1j * im if usecomplex else [re, im]
 
     # north
-    lon, lat = deflect(0.0, 0.0, alpha(r, 0, usecomplex=usecomplex))
+    lon, lat = glass.deflect(0.0, 0.0, alpha(r, 0, usecomplex=usecomplex))
     np.testing.assert_allclose([lon, lat], [0.0, d], atol=1e-15)
 
     # south
-    lon, lat = deflect(0.0, 0.0, alpha(-r, 0, usecomplex=usecomplex))
+    lon, lat = glass.deflect(0.0, 0.0, alpha(-r, 0, usecomplex=usecomplex))
     np.testing.assert_allclose([lon, lat], [0.0, -d], atol=1e-15)
 
     # east
-    lon, lat = deflect(0.0, 0.0, alpha(0, r, usecomplex=usecomplex))
+    lon, lat = glass.deflect(0.0, 0.0, alpha(0, r, usecomplex=usecomplex))
     np.testing.assert_allclose([lon, lat], [-d, 0.0], atol=1e-15)
 
     # west
-    lon, lat = deflect(0.0, 0.0, alpha(0, -r, usecomplex=usecomplex))
+    lon, lat = glass.deflect(0.0, 0.0, alpha(0, -r, usecomplex=usecomplex))
     np.testing.assert_allclose([lon, lat], [d, 0.0], atol=1e-15)
 
 
@@ -142,7 +134,7 @@ def test_deflect_many(rng: np.random.Generator) -> None:
     lon_ = np.degrees(rng.uniform(-np.pi, np.pi, size=n))
     lat_ = np.degrees(np.arcsin(rng.uniform(-1, 1, size=n)))
 
-    lon, lat = deflect(lon_, lat_, abs_alpha * np.exp(1j * arg_alpha))
+    lon, lat = glass.deflect(lon_, lat_, abs_alpha * np.exp(1j * arg_alpha))
 
     x_, y_, z_ = healpix.ang2vec(lon_, lat_, lonlat=True)
     x, y, z = healpix.ang2vec(lon, lat, lonlat=True)
diff --git a/tests/test_observations.py b/tests/test_observations.py
index 797dc871..a1066413 100644
--- a/tests/test_observations.py
+++ b/tests/test_observations.py
@@ -2,65 +2,58 @@
 import numpy as np
 import pytest
 
-from glass import (
-    equal_dens_zbins,
-    fixed_zbins,
-    gaussian_nz,
-    smail_nz,
-    tomo_nz_gausserr,
-    vmap_galactic_ecliptic,
-)
+import glass
 
 
 def test_vmap_galactic_ecliptic() -> None:
-    """Add unit tests for :func:`vmap_galactic_ecliptic`."""
+    """Add unit tests for :func:`glass.vmap_galactic_ecliptic`."""
     n_side = 4
 
     # check shape
 
-    vmap = vmap_galactic_ecliptic(n_side)
+    vmap = glass.vmap_galactic_ecliptic(n_side)
     np.testing.assert_array_equal(len(vmap), healpix.nside2npix(n_side))
 
     # no rotation
 
-    vmap = vmap_galactic_ecliptic(n_side, galactic=(0, 0), ecliptic=(0, 0))
+    vmap = glass.vmap_galactic_ecliptic(n_side, galactic=(0, 0), ecliptic=(0, 0))
     np.testing.assert_array_equal(vmap, np.zeros_like(vmap))
 
     # check errors raised
 
     with pytest.raises(TypeError, match="galactic stripe must be a pair of numbers"):
-        vmap_galactic_ecliptic(n_side, galactic=(1,))  # type: ignore[arg-type]
+        glass.vmap_galactic_ecliptic(n_side, galactic=(1,))  # type: ignore[arg-type]
 
     with pytest.raises(TypeError, match="ecliptic stripe must be a pair of numbers"):
-        vmap_galactic_ecliptic(n_side, ecliptic=(1,))  # type: ignore[arg-type]
+        glass.vmap_galactic_ecliptic(n_side, ecliptic=(1,))  # type: ignore[arg-type]
 
     with pytest.raises(TypeError, match="galactic stripe must be a pair of numbers"):
-        vmap_galactic_ecliptic(n_side, galactic=(1, 2, 3))  # type: ignore[arg-type]
+        glass.vmap_galactic_ecliptic(n_side, galactic=(1, 2, 3))  # type: ignore[arg-type]
 
     with pytest.raises(TypeError, match="ecliptic stripe must be a pair of numbers"):
-        vmap_galactic_ecliptic(n_side, ecliptic=(1, 2, 3))  # type: ignore[arg-type]
+        glass.vmap_galactic_ecliptic(n_side, ecliptic=(1, 2, 3))  # type: ignore[arg-type]
 
 
 def test_gaussian_nz(rng: np.random.Generator) -> None:
-    """Add unit tests for :func:`gaussian_nz`."""
+    """Add unit tests for :func:`glass.gaussian_nz`."""
     mean = 0
     sigma = 1
     z = np.linspace(0, 1, 11)
 
     # check passing in the norm
 
-    nz = gaussian_nz(z, mean, sigma, norm=0)
+    nz = glass.gaussian_nz(z, mean, sigma, norm=0)
     np.testing.assert_array_equal(nz, np.zeros_like(nz))
 
     # check the value of each entry is close to the norm
 
     norm = 1
-    nz = gaussian_nz(z, mean, sigma, norm=norm)
+    nz = glass.gaussian_nz(z, mean, sigma, norm=norm)
     np.testing.assert_allclose(nz.sum() / nz.shape, norm, rtol=1e-2)
 
     # check multidimensionality size
 
-    nz = gaussian_nz(
+    nz = glass.gaussian_nz(
         z,
         np.tile(mean, z.shape),
         np.tile(sigma, z.shape),
@@ -70,7 +63,7 @@ def test_gaussian_nz(rng: np.random.Generator) -> None:
 
 
 def test_smail_nz() -> None:
-    """Add unit tests for :func:`smail_nz`."""
+    """Add unit tests for :func:`glass.smail_nz`."""
     alpha = 1
     beta = 1
     mode = 1
@@ -78,12 +71,12 @@ def test_smail_nz() -> None:
 
     # check passing in the norm
 
-    pz = smail_nz(z, mode, alpha, beta, norm=0)
+    pz = glass.smail_nz(z, mode, alpha, beta, norm=0)
     np.testing.assert_array_equal(pz, np.zeros_like(pz))
 
 
 def test_fixed_zbins() -> None:
-    """Add unit tests for :func:`fixed_zbins`."""
+    """Add unit tests for :func:`glass.fixed_zbins`."""
     zmin = 0
     zmax = 1
 
@@ -91,37 +84,37 @@ def test_fixed_zbins() -> None:
 
     nbins = 5
     expected_zbins = [(0.0, 0.2), (0.2, 0.4), (0.4, 0.6), (0.6, 0.8), (0.8, 1.0)]
-    zbins = fixed_zbins(zmin, zmax, nbins=nbins)
+    zbins = glass.fixed_zbins(zmin, zmax, nbins=nbins)
     np.testing.assert_array_equal(len(zbins), nbins)
     np.testing.assert_allclose(zbins, expected_zbins, rtol=1e-15)
 
     # check dz input
 
     dz = 0.2
-    zbins = fixed_zbins(zmin, zmax, dz=dz)
+    zbins = glass.fixed_zbins(zmin, zmax, dz=dz)
     np.testing.assert_array_equal(len(zbins), np.ceil((zmax - zmin) / dz))
     np.testing.assert_allclose(zbins, expected_zbins, rtol=1e-15)
 
     # check dz for spacing which results in a max value above zmax
 
-    zbins = fixed_zbins(zmin, zmax, dz=0.3)
+    zbins = glass.fixed_zbins(zmin, zmax, dz=0.3)
     np.testing.assert_array_less(zmax, zbins[-1][1])
 
     # check error raised
 
     with pytest.raises(ValueError, match="exactly one of nbins and dz must be given"):
-        fixed_zbins(zmin, zmax, nbins=nbins, dz=dz)
+        glass.fixed_zbins(zmin, zmax, nbins=nbins, dz=dz)
 
 
 def test_equal_dens_zbins() -> None:
-    """Add unit tests for :func:`equal_dens_zbins`."""
+    """Add unit tests for :func:`glass.equal_dens_zbins`."""
     z = np.linspace(0, 1, 11)
     nbins = 5
 
     # check expected zbins returned
 
     expected_zbins = [(0.0, 0.2), (0.2, 0.4), (0.4, 0.6), (0.6, 0.8), (0.8, 1.0)]
-    zbins = equal_dens_zbins(z, np.ones_like(z), nbins)
+    zbins = glass.equal_dens_zbins(z, np.ones_like(z), nbins)
     np.testing.assert_allclose(zbins, expected_zbins, rtol=1e-15)
 
     # check output shape
@@ -130,14 +123,14 @@ def test_equal_dens_zbins() -> None:
 
 
 def test_tomo_nz_gausserr() -> None:
-    """Add unit tests for :func:`tomo_nz_gausserr`."""
+    """Add unit tests for :func:`glass.tomo_nz_gausserr`."""
     sigma_0 = 0.1
     z = np.linspace(0, 1, 11)
     zbins = [(0, 0.2), (0.2, 0.4), (0.4, 0.6), (0.6, 0.8), (0.8, 1.0)]
 
     # check zeros returned
 
-    binned_nz = tomo_nz_gausserr(z, np.zeros_like(z), sigma_0, zbins)
+    binned_nz = glass.tomo_nz_gausserr(z, np.zeros_like(z), sigma_0, zbins)
     np.testing.assert_array_equal(binned_nz, np.zeros_like(binned_nz))
 
     # check the shape of the output
diff --git a/tests/test_points.py b/tests/test_points.py
index 4a05d316..acc56846 100644
--- a/tests/test_points.py
+++ b/tests/test_points.py
@@ -6,14 +6,7 @@
 import numpy as np
 import pytest
 
-from glass import (
-    effective_bias,
-    linear_bias,
-    loglinear_bias,
-    position_weights,
-    positions_from_delta,
-    uniform_positions,
-)
+import glass
 
 if TYPE_CHECKING:
     from collections.abc import Generator
@@ -54,17 +47,17 @@ def test_effective_bias(mocker: pytest_mock.MockerFixture) -> None:
     z = np.linspace(0, 1, 10)
     bz = np.zeros((10,))
 
-    np.testing.assert_allclose(effective_bias(z, bz, w), np.zeros((10,)))
+    np.testing.assert_allclose(glass.effective_bias(z, bz, w), np.zeros((10,)))
 
     z = np.zeros((10,))
     bz = np.full_like(z, 0.5)
 
-    np.testing.assert_allclose(effective_bias(z, bz, w), np.zeros((10,)))
+    np.testing.assert_allclose(glass.effective_bias(z, bz, w), np.zeros((10,)))
 
     z = np.linspace(0, 1, 10)
     bz = np.full_like(z, 0.5)
 
-    np.testing.assert_allclose(effective_bias(z, bz, w), 0.25)
+    np.testing.assert_allclose(glass.effective_bias(z, bz, w), 0.25)
 
 
 def test_linear_bias(rng: np.random.Generator) -> None:
@@ -73,21 +66,21 @@ def test_linear_bias(rng: np.random.Generator) -> None:
     delta = np.zeros((2, 2))
     b = 2.0
 
-    np.testing.assert_allclose(linear_bias(delta, b), np.zeros((2, 2)))
+    np.testing.assert_allclose(glass.linear_bias(delta, b), np.zeros((2, 2)))
 
     # test with 0 b
 
     delta = rng.normal(5, 1, size=(2, 2))
     b = 0.0
 
-    np.testing.assert_allclose(linear_bias(delta, b), np.zeros((2, 2)))
+    np.testing.assert_allclose(glass.linear_bias(delta, b), np.zeros((2, 2)))
 
     # compare with original implementation
 
     delta = rng.normal(5, 1, size=(2, 2))
     b = 2.0
 
-    np.testing.assert_allclose(linear_bias(delta, b), b * delta)
+    np.testing.assert_allclose(glass.linear_bias(delta, b), b * delta)
 
 
 def test_loglinear_bias(rng: np.random.Generator) -> None:
@@ -96,21 +89,23 @@ def test_loglinear_bias(rng: np.random.Generator) -> None:
     delta = np.zeros((2, 2))
     b = 2.0
 
-    np.testing.assert_allclose(loglinear_bias(delta, b), np.zeros((2, 2)))
+    np.testing.assert_allclose(glass.loglinear_bias(delta, b), np.zeros((2, 2)))
 
     # test with 0 b
 
     delta = rng.normal(5, 1, size=(2, 2))
     b = 0.0
 
-    np.testing.assert_allclose(loglinear_bias(delta, b), np.zeros((2, 2)))
+    np.testing.assert_allclose(glass.loglinear_bias(delta, b), np.zeros((2, 2)))
 
     # compare with numpy implementation
 
     delta = rng.normal(5, 1, size=(2, 2))
     b = 2.0
 
-    np.testing.assert_allclose(loglinear_bias(delta, b), np.expm1(b * np.log1p(delta)))
+    np.testing.assert_allclose(
+        glass.loglinear_bias(delta, b), np.expm1(b * np.log1p(delta))
+    )
 
 
 def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR0915
@@ -125,14 +120,14 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
     bias = 0.8
     vis = np.ones(npix)
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))
+    lon, lat, cnt = catpos(glass.positions_from_delta(ngal, delta, bias, vis))
 
     assert isinstance(cnt, int)
     assert lon.shape == lat.shape == (cnt,)
 
     # test with rng
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis, rng=rng))
+    lon, lat, cnt = catpos(glass.positions_from_delta(ngal, delta, bias, vis, rng=rng))
 
     assert isinstance(cnt, int)
     assert lon.shape == lat.shape == (cnt,)
@@ -140,7 +135,7 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
     # case: Nons bias and callable bias model
 
     lon, lat, cnt = catpos(
-        positions_from_delta(ngal, delta, None, vis, bias_model=lambda x: x)
+        glass.positions_from_delta(ngal, delta, None, vis, bias_model=lambda x: x)
     )
 
     assert isinstance(cnt, int)
@@ -148,7 +143,7 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
 
     # case: None vis
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, None))
+    lon, lat, cnt = catpos(glass.positions_from_delta(ngal, delta, bias, None))
 
     assert isinstance(cnt, int)
     assert lon.shape == lat.shape == (cnt,)
@@ -156,7 +151,7 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
     # case: remove monopole
 
     lon, lat, cnt = catpos(
-        positions_from_delta(ngal, delta, bias, vis, remove_monopole=True)
+        glass.positions_from_delta(ngal, delta, bias, vis, remove_monopole=True)
     )
 
     assert isinstance(cnt, int)
@@ -165,7 +160,7 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
     # case: negative delta
 
     lon, lat, cnt = catpos(
-        positions_from_delta(ngal, np.linspace(-1, -1, npix), None, vis)
+        glass.positions_from_delta(ngal, np.linspace(-1, -1, npix), None, vis)
     )
 
     assert isinstance(cnt, int)
@@ -175,7 +170,7 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
     # case: large delta
 
     lon, lat, cnt = catpos(
-        positions_from_delta(ngal, rng.normal(100, 1, size=(npix,)), bias, vis)
+        glass.positions_from_delta(ngal, rng.normal(100, 1, size=(npix,)), bias, vis)
     )
 
     assert isinstance(cnt, int)
@@ -188,7 +183,7 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
     bias = 0.8
     vis = np.ones(12)
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))
+    lon, lat, cnt = catpos(glass.positions_from_delta(ngal, delta, bias, vis))
 
     assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (2,)
@@ -202,7 +197,7 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
     bias = 0.8
     vis = np.ones(12)
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))
+    lon, lat, cnt = catpos(glass.positions_from_delta(ngal, delta, bias, vis))
 
     assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (3, 2)
@@ -216,7 +211,7 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
     bias = 0.8
     vis = np.ones(12)
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))
+    lon, lat, cnt = catpos(glass.positions_from_delta(ngal, delta, bias, vis))
 
     assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (3, 2)
@@ -227,7 +222,7 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
 
     vis[: vis.size // 2] = 0.0
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))
+    lon, lat, cnt = catpos(glass.positions_from_delta(ngal, delta, bias, vis))
 
     assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (3, 2)
@@ -237,7 +232,7 @@ def test_positions_from_delta(rng: np.random.Generator) -> None:  # noqa: PLR091
     # test TypeError
 
     with pytest.raises(TypeError, match="bias_model must be string or callable"):
-        next(positions_from_delta(ngal, delta, bias, vis, bias_model=0))  # type: ignore[arg-type]
+        next(glass.positions_from_delta(ngal, delta, bias, vis, bias_model=0))  # type: ignore[arg-type]
 
 
 def test_uniform_positions(rng: np.random.Generator) -> None:
@@ -245,11 +240,11 @@ def test_uniform_positions(rng: np.random.Generator) -> None:
 
     ngal: float | NDArray[np.float64] = 1e-3
 
-    lon, lat, cnt = catpos(uniform_positions(ngal))
+    lon, lat, cnt = catpos(glass.uniform_positions(ngal))
 
     # test with rng
 
-    lon, lat, cnt = catpos(uniform_positions(ngal, rng=rng))
+    lon, lat, cnt = catpos(glass.uniform_positions(ngal, rng=rng))
 
     assert isinstance(cnt, int)
     assert lon.shape == lat.shape == (cnt,)
@@ -258,7 +253,7 @@ def test_uniform_positions(rng: np.random.Generator) -> None:
 
     ngal = np.array([1e-3, 2e-3, 3e-3])
 
-    lon, lat, cnt = catpos(uniform_positions(ngal))
+    lon, lat, cnt = catpos(glass.uniform_positions(ngal))
 
     assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (3,)
@@ -268,7 +263,7 @@ def test_uniform_positions(rng: np.random.Generator) -> None:
 
     ngal = np.array([[1e-3, 2e-3], [3e-3, 4e-3], [5e-3, 6e-3]])
 
-    lon, lat, cnt = catpos(uniform_positions(ngal))
+    lon, lat, cnt = catpos(glass.uniform_positions(ngal))
 
     assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (3, 2)
@@ -281,7 +276,7 @@ def test_position_weights(rng: np.random.Generator) -> None:
             counts = rng.random(cshape)
             bias = None if bshape is None else rng.random(bshape)
 
-            weights = position_weights(counts, bias)
+            weights = glass.position_weights(counts, bias)
 
             expected = counts / counts.sum(axis=0, keepdims=True)
             if bias is not None:
diff --git a/tests/test_shapes.py b/tests/test_shapes.py
index 71892455..d01f2eb5 100644
--- a/tests/test_shapes.py
+++ b/tests/test_shapes.py
@@ -1,39 +1,36 @@
 import numpy as np
 import pytest
 
-from glass import (
-    ellipticity_gaussian,
-    ellipticity_intnorm,
-    ellipticity_ryden04,
-    triaxial_axis_ratio,
-)
+import glass
 
 
 def test_triaxial_axis_ratio(rng: np.random.Generator) -> None:
     # single axis ratio
 
-    q = triaxial_axis_ratio(0.8, 0.4)
+    q = glass.triaxial_axis_ratio(0.8, 0.4)
     assert np.isscalar(q)
 
     # many axis ratios
 
-    q = triaxial_axis_ratio(0.8, 0.4, size=1000)
+    q = glass.triaxial_axis_ratio(0.8, 0.4, size=1000)
     assert np.shape(q) == (1000,)
 
     # explicit shape
 
-    q = triaxial_axis_ratio(0.8, 0.4, size=(10, 10))
+    q = glass.triaxial_axis_ratio(0.8, 0.4, size=(10, 10))
     assert np.shape(q) == (10, 10)
 
     # implicit size
 
-    q1 = triaxial_axis_ratio(np.array([0.8, 0.9]), 0.4)
-    q2 = triaxial_axis_ratio(0.8, np.array([0.4, 0.5]))
+    q1 = glass.triaxial_axis_ratio(np.array([0.8, 0.9]), 0.4)
+    q2 = glass.triaxial_axis_ratio(0.8, np.array([0.4, 0.5]))
     assert np.shape(q1) == np.shape(q2) == (2,)
 
     # broadcasting rule
 
-    q = triaxial_axis_ratio(np.array([[0.6, 0.7], [0.8, 0.9]]), np.array([0.4, 0.5]))
+    q = glass.triaxial_axis_ratio(
+        np.array([[0.6, 0.7], [0.8, 0.9]]), np.array([0.4, 0.5])
+    )
     assert np.shape(q) == (2, 2)
 
     # random parameters and check that projection is
@@ -42,42 +39,42 @@ def test_triaxial_axis_ratio(rng: np.random.Generator) -> None:
     zeta, xi = np.sort(rng.uniform(0, 1, size=(2, 1000)), axis=0)
     qmin = np.min([zeta, xi, xi / zeta], axis=0)
     qmax = np.max([zeta, xi, xi / zeta], axis=0)
-    q = triaxial_axis_ratio(zeta, xi)
+    q = glass.triaxial_axis_ratio(zeta, xi)
     assert np.all((qmax >= q) & (q >= qmin))
 
 
 def test_ellipticity_ryden04(rng: np.random.Generator) -> None:
     # single ellipticity
 
-    e = ellipticity_ryden04(-1.85, 0.89, 0.222, 0.056)
+    e = glass.ellipticity_ryden04(-1.85, 0.89, 0.222, 0.056)
     assert np.isscalar(e)
 
     # test with rng
 
-    e = ellipticity_ryden04(-1.85, 0.89, 0.222, 0.056, rng=rng)
+    e = glass.ellipticity_ryden04(-1.85, 0.89, 0.222, 0.056, rng=rng)
     assert np.isscalar(e)
 
     # many ellipticities
 
-    e = ellipticity_ryden04(-1.85, 0.89, 0.222, 0.056, size=1000)
+    e = glass.ellipticity_ryden04(-1.85, 0.89, 0.222, 0.056, size=1000)
     assert np.shape(e) == (1000,)
 
     # explicit shape
 
-    e = ellipticity_ryden04(-1.85, 0.89, 0.222, 0.056, size=(10, 10))
+    e = glass.ellipticity_ryden04(-1.85, 0.89, 0.222, 0.056, size=(10, 10))
     assert np.shape(e) == (10, 10)
 
     # implicit size
 
-    e1 = ellipticity_ryden04(-1.85, 0.89, np.array([0.222, 0.333]), 0.056)
-    e2 = ellipticity_ryden04(-1.85, 0.89, 0.222, np.array([0.056, 0.067]))
-    e3 = ellipticity_ryden04(np.array([-1.85, -2.85]), 0.89, 0.222, 0.056)
-    e4 = ellipticity_ryden04(-1.85, np.array([0.89, 1.001]), 0.222, 0.056)
+    e1 = glass.ellipticity_ryden04(-1.85, 0.89, np.array([0.222, 0.333]), 0.056)
+    e2 = glass.ellipticity_ryden04(-1.85, 0.89, 0.222, np.array([0.056, 0.067]))
+    e3 = glass.ellipticity_ryden04(np.array([-1.85, -2.85]), 0.89, 0.222, 0.056)
+    e4 = glass.ellipticity_ryden04(-1.85, np.array([0.89, 1.001]), 0.222, 0.056)
     assert np.shape(e1) == np.shape(e2) == np.shape(e3) == np.shape(e4) == (2,)
 
     # broadcasting rule
 
-    e = ellipticity_ryden04(
+    e = glass.ellipticity_ryden04(
         np.array([-1.9, -2.9]),
         0.9,
         np.array([[0.2, 0.3], [0.4, 0.5]]),
@@ -87,10 +84,10 @@ def test_ellipticity_ryden04(rng: np.random.Generator) -> None:
 
     # check that result is in the specified range
 
-    e = ellipticity_ryden04(0.0, 1.0, 0.222, 0.056, size=10)
+    e = glass.ellipticity_ryden04(0.0, 1.0, 0.222, 0.056, size=10)
     assert np.all((e.real >= -1.0) & (e.real <= 1.0))
 
-    e = ellipticity_ryden04(0.0, 1.0, 0.0, 1.0, size=10)
+    e = glass.ellipticity_ryden04(0.0, 1.0, 0.0, 1.0, size=10)
     assert np.all((e.real >= -1.0) & (e.real <= 1.0))
 
 
@@ -98,13 +95,13 @@ def test_ellipticity_ryden04(rng: np.random.Generator) -> None:
 def test_ellipticity_gaussian(rng: np.random.Generator) -> None:
     n = 1_000_000
 
-    eps = ellipticity_gaussian(n, 0.256)
+    eps = glass.ellipticity_gaussian(n, 0.256)
 
     assert eps.shape == (n,)
 
     # test with rng
 
-    eps = ellipticity_gaussian(n, 0.256, rng=rng)
+    eps = glass.ellipticity_gaussian(n, 0.256, rng=rng)
 
     assert eps.shape == (n,)
 
@@ -113,7 +110,7 @@ def test_ellipticity_gaussian(rng: np.random.Generator) -> None:
     np.testing.assert_allclose(np.std(eps.real), 0.256, atol=1e-3, rtol=0)
     np.testing.assert_allclose(np.std(eps.imag), 0.256, atol=1e-3, rtol=0)
 
-    eps = ellipticity_gaussian(np.array([n, n]), np.array([0.128, 0.256]))
+    eps = glass.ellipticity_gaussian(np.array([n, n]), np.array([0.128, 0.256]))
 
     assert eps.shape == (2 * n,)
 
@@ -128,13 +125,13 @@ def test_ellipticity_gaussian(rng: np.random.Generator) -> None:
 def test_ellipticity_intnorm(rng: np.random.Generator) -> None:
     n = 1_000_000
 
-    eps = ellipticity_intnorm(n, 0.256)
+    eps = glass.ellipticity_intnorm(n, 0.256)
 
     assert eps.shape == (n,)
 
     # test with rng
 
-    eps = ellipticity_intnorm(n, 0.256, rng=rng)
+    eps = glass.ellipticity_intnorm(n, 0.256, rng=rng)
 
     assert eps.shape == (n,)
 
@@ -143,7 +140,7 @@ def test_ellipticity_intnorm(rng: np.random.Generator) -> None:
     np.testing.assert_allclose(np.std(eps.real), 0.256, atol=1e-3, rtol=0)
     np.testing.assert_allclose(np.std(eps.imag), 0.256, atol=1e-3, rtol=0)
 
-    eps = ellipticity_intnorm(np.array([n, n]), np.array([0.128, 0.256]))
+    eps = glass.ellipticity_intnorm(np.array([n, n]), np.array([0.128, 0.256]))
 
     assert eps.shape == (2 * n,)
 
@@ -155,4 +152,4 @@ def test_ellipticity_intnorm(rng: np.random.Generator) -> None:
     np.testing.assert_allclose(np.std(eps.imag[n:]), 0.256, atol=1e-3, rtol=0)
 
     with pytest.raises(ValueError, match="sigma must be between"):
-        ellipticity_intnorm(1, 0.71)
+        glass.ellipticity_intnorm(1, 0.71)
diff --git a/tests/test_shells.py b/tests/test_shells.py
index 8f87bdf8..ab97f83e 100644
--- a/tests/test_shells.py
+++ b/tests/test_shells.py
@@ -3,29 +3,16 @@
 
 from cosmology import Cosmology
 
-from glass import (
-    RadialWindow,
-    combine,  # noqa: F401
-    cubic_windows,
-    density_weight,
-    distance_grid,
-    distance_weight,
-    linear_windows,
-    partition,
-    redshift_grid,
-    restrict,
-    tophat_windows,
-    volume_weight,
-)
+import glass
 
 
 def test_distance_weight(cosmo: Cosmology) -> None:
-    """Add unit tests for :func:`distance_weight`."""
+    """Add unit tests for :func:`glass.distance_weight`."""
     z = np.linspace(0, 1, 6)
 
     # check shape
 
-    w = distance_weight(z, cosmo)
+    w = glass.distance_weight(z, cosmo)
     np.testing.assert_array_equal(w.shape, z.shape)
 
     # check first value is 1
@@ -38,12 +25,12 @@ def test_distance_weight(cosmo: Cosmology) -> None:
 
 
 def test_volume_weight(cosmo: Cosmology) -> None:
-    """Add unit tests for :func:`volume_weight`."""
+    """Add unit tests for :func:`glass.volume_weight`."""
     z = np.linspace(0, 1, 6)
 
     # check shape
 
-    w = volume_weight(z, cosmo)
+    w = glass.volume_weight(z, cosmo)
     np.testing.assert_array_equal(w.shape, z.shape)
 
     # check first value is 0
@@ -56,12 +43,12 @@ def test_volume_weight(cosmo: Cosmology) -> None:
 
 
 def test_density_weight(cosmo: Cosmology) -> None:
-    """Add unit tests for :func:`density_weight`."""
+    """Add unit tests for :func:`glass.density_weight`."""
     z = np.linspace(0, 1, 6)
 
     # check shape
 
-    w = density_weight(z, cosmo)
+    w = glass.density_weight(z, cosmo)
     np.testing.assert_array_equal(w.shape, z.shape)
 
     # check first value is 0
@@ -74,11 +61,11 @@ def test_density_weight(cosmo: Cosmology) -> None:
 
 
 def test_tophat_windows() -> None:
-    """Add unit tests for :func:`tophat_windows`."""
+    """Add unit tests for :func:`glass.tophat_windows`."""
     zb = np.array([0.0, 0.1, 0.2, 0.5, 1.0, 2.0])
     dz = 0.005
 
-    ws = tophat_windows(zb, dz)
+    ws = glass.tophat_windows(zb, dz)
 
     assert len(ws) == len(zb) - 1
 
@@ -96,7 +83,7 @@ def test_tophat_windows() -> None:
 
 
 def test_linear_windows() -> None:
-    """Add unit tests for :func:`linear_windows`."""
+    """Add unit tests for :func:`glass.linear_windows`."""
     dz = 1e-2
     zgrid = [
         0.0,
@@ -108,7 +95,7 @@ def test_linear_windows() -> None:
 
     # check spacing of redshift grid
 
-    ws = linear_windows(zgrid)
+    ws = glass.linear_windows(zgrid)
     np.testing.assert_allclose(dz, np.diff(ws[0].za).mean(), atol=1e-2)
 
     # check number of windows
@@ -121,7 +108,7 @@ def test_linear_windows() -> None:
 
     # check weight function input
 
-    ws = linear_windows(
+    ws = glass.linear_windows(
         zgrid,
         weight=lambda _: 0,  # type: ignore[arg-type, return-value]
     )
@@ -131,18 +118,18 @@ def test_linear_windows() -> None:
     # check error raised
 
     with pytest.raises(ValueError, match="nodes must have at least 3 entries"):
-        linear_windows([])
+        glass.linear_windows([])
 
     # check warning issued
 
     with pytest.warns(
         UserWarning, match="first triangular window does not start at z=0"
     ):
-        linear_windows([0.1, 0.2, 0.3])
+        glass.linear_windows([0.1, 0.2, 0.3])
 
 
 def test_cubic_windows() -> None:
-    """Add unit tests for :func:`cubic_windows`."""
+    """Add unit tests for :func:`glass.cubic_windows`."""
     dz = 1e-2
     zgrid = [
         0.0,
@@ -154,7 +141,7 @@ def test_cubic_windows() -> None:
 
     # check spacing of redshift grid
 
-    ws = cubic_windows(zgrid)
+    ws = glass.cubic_windows(zgrid)
     np.testing.assert_allclose(dz, np.diff(ws[0].za).mean(), atol=1e-2)
 
     # check number of windows
@@ -167,7 +154,7 @@ def test_cubic_windows() -> None:
 
     # check weight function input
 
-    ws = cubic_windows(
+    ws = glass.cubic_windows(
         zgrid,
         weight=lambda _: 0,  # type: ignore[arg-type, return-value]
     )
@@ -177,29 +164,29 @@ def test_cubic_windows() -> None:
     # check error raised
 
     with pytest.raises(ValueError, match="nodes must have at least 3 entries"):
-        cubic_windows([])
+        glass.cubic_windows([])
 
     # check warning issued
 
     with pytest.warns(
         UserWarning, match="first cubic spline window does not start at z=0"
     ):
-        cubic_windows([0.1, 0.2, 0.3])
+        glass.cubic_windows([0.1, 0.2, 0.3])
 
 
 def test_restrict() -> None:
-    """Add unit tests for :func:`restrict`."""
+    """Add unit tests for :func:`glass.restrict`."""
     # Gaussian test function
     z = np.linspace(0.0, 5.0, 1000)
     f = np.exp(-(((z - 2.0) / 0.5) ** 2) / 2)
 
     # window for restriction
-    w = RadialWindow(
+    w = glass.RadialWindow(
         za=np.array([1.0, 2.0, 3.0, 4.0]),
         wa=np.array([0.0, 0.5, 0.5, 0.0]),
     )
 
-    zr, fr = restrict(z, f, w)
+    zr, fr = glass.restrict(z, f, w)
 
     assert zr[0] == w.za[0]
     assert zr[-1] == w.za[-1]
@@ -220,14 +207,14 @@ def test_restrict() -> None:
 
 @pytest.mark.parametrize("method", ["lstsq", "nnls", "restrict"])
 def test_partition(method: str) -> None:
-    """Add unit tests for :func:`partition`."""
+    """Add unit tests for :func:`glass.partition`."""
     shells = [
-        RadialWindow(np.array([0.0, 1.0]), np.array([1.0, 0.0]), 0.0),
-        RadialWindow(np.array([0.0, 1.0, 2.0]), np.array([0.0, 1.0, 0.0]), 0.5),
-        RadialWindow(np.array([1.0, 2.0, 3.0]), np.array([0.0, 1.0, 0.0]), 1.5),
-        RadialWindow(np.array([2.0, 3.0, 4.0]), np.array([0.0, 1.0, 0.0]), 2.5),
-        RadialWindow(np.array([3.0, 4.0, 5.0]), np.array([0.0, 1.0, 0.0]), 3.5),
-        RadialWindow(np.array([4.0, 5.0]), np.array([0.0, 1.0]), 5.0),
+        glass.RadialWindow(np.array([0.0, 1.0]), np.array([1.0, 0.0]), 0.0),
+        glass.RadialWindow(np.array([0.0, 1.0, 2.0]), np.array([0.0, 1.0, 0.0]), 0.5),
+        glass.RadialWindow(np.array([1.0, 2.0, 3.0]), np.array([0.0, 1.0, 0.0]), 1.5),
+        glass.RadialWindow(np.array([2.0, 3.0, 4.0]), np.array([0.0, 1.0, 0.0]), 2.5),
+        glass.RadialWindow(np.array([3.0, 4.0, 5.0]), np.array([0.0, 1.0, 0.0]), 3.5),
+        glass.RadialWindow(np.array([4.0, 5.0]), np.array([0.0, 1.0]), 5.0),
     ]
 
     z = np.linspace(0.0, 5.0, 1000)
@@ -236,7 +223,7 @@ def test_partition(method: str) -> None:
 
     assert fz.shape == (3, 2, 1000)
 
-    part = partition(z, fz, shells, method=method)
+    part = glass.partition(z, fz, shells, method=method)
 
     assert part.shape == (len(shells), 3, 2)
 
@@ -244,25 +231,25 @@ def test_partition(method: str) -> None:
 
 
 def test_redshift_grid() -> None:
-    """Add unit tests for :func:`redshift_grid`."""
+    """Add unit tests for :func:`glass.redshift_grid`."""
     zmin = 0
     zmax = 1
 
     # check num input
 
     num = 5
-    z = redshift_grid(zmin, zmax, num=5)
+    z = glass.redshift_grid(zmin, zmax, num=5)
     assert len(z) == num + 1
 
     # check dz input
 
     dz = 0.2
-    z = redshift_grid(zmin, zmax, dz=dz)
+    z = glass.redshift_grid(zmin, zmax, dz=dz)
     assert len(z) == np.ceil((zmax - zmin) / dz) + 1
 
     # check dz for spacing which results in a max value above zmax
 
-    z = redshift_grid(zmin, zmax, dz=0.3)
+    z = glass.redshift_grid(zmin, zmax, dz=0.3)
     assert zmax < z[-1]
 
     # check error raised
@@ -271,35 +258,35 @@ def test_redshift_grid() -> None:
         ValueError,
         match="exactly one of grid step size or number of steps must be given",
     ):
-        redshift_grid(zmin, zmax)
+        glass.redshift_grid(zmin, zmax)
 
     with pytest.raises(
         ValueError,
         match="exactly one of grid step size or number of steps must be given",
     ):
-        redshift_grid(zmin, zmax, dz=dz, num=num)
+        glass.redshift_grid(zmin, zmax, dz=dz, num=num)
 
 
 def test_distance_grid(cosmo: Cosmology) -> None:
-    """Add unit tests for :func:`distance_grid`."""
+    """Add unit tests for :func:`glass.distance_grid`."""
     zmin = 0
     zmax = 1
 
     # check num input
 
     num = 5
-    x = distance_grid(cosmo, zmin, zmax, num=5)
+    x = glass.distance_grid(cosmo, zmin, zmax, num=5)
     assert len(x) == num + 1
 
     # check dz input
 
     dx = 0.2
-    x = distance_grid(cosmo, zmin, zmax, dx=dx)
+    x = glass.distance_grid(cosmo, zmin, zmax, dx=dx)
     assert len(x) == np.ceil((zmax - zmin) / dx) + 1
 
     # check decrease in distance
 
-    x = distance_grid(cosmo, zmin, zmax, dx=0.3)
+    x = glass.distance_grid(cosmo, zmin, zmax, dx=0.3)
     np.testing.assert_array_less(x[1:], x[:-1])
 
     # check error raised
@@ -308,14 +295,14 @@ def test_distance_grid(cosmo: Cosmology) -> None:
         ValueError,
         match="exactly one of grid step size or number of steps must be given",
     ):
-        distance_grid(cosmo, zmin, zmax)
+        glass.distance_grid(cosmo, zmin, zmax)
 
     with pytest.raises(
         ValueError,
         match="exactly one of grid step size or number of steps must be given",
     ):
-        distance_grid(cosmo, zmin, zmax, dx=dx, num=num)
+        glass.distance_grid(cosmo, zmin, zmax, dx=dx, num=num)
 
 
 def test_combine() -> None:
-    """Add unit tests for :func:`combine`."""
+    """Add unit tests for :func:`glass.combine`."""
diff --git a/tests/test_user.py b/tests/test_user.py
index 754c7c1e..02f5a69a 100644
--- a/tests/test_user.py
+++ b/tests/test_user.py
@@ -5,7 +5,7 @@
 import pytest
 from numpy.typing import NDArray
 
-from glass import load_cls, save_cls, write_catalog
+import glass
 
 # check if fitsio is available for testing
 HAVE_FITSIO = importlib.util.find_spec("fitsio") is not None
@@ -19,7 +19,7 @@
 
 def test_read_write_cls(rng: np.random.Generator, tmp_path: pathlib.Path) -> None:
     cls = rng.normal(size=(10, 10))
-    save_cls(tmp_path / cls_file, cls)
+    glass.save_cls(tmp_path / cls_file, cls)
 
     assert pathlib.Path.exists(tmp_path / cls_file)
 
@@ -31,7 +31,7 @@ def test_read_write_cls(rng: np.random.Generator, tmp_path: pathlib.Path) -> Non
     np.testing.assert_array_equal(split, np.cumsum([len(cl) for cl in cls[:-1]]))
     np.testing.assert_array_equal(cls, np.split(values, split))
 
-    npz = load_cls(tmp_path / cls_file)
+    npz = glass.load_cls(tmp_path / cls_file)
     np.testing.assert_array_equal(npz, cls)
 
 
@@ -56,7 +56,7 @@ def _test_append(
             hdu.write(data, names=names, firstrow=hdu.get_nrows())
 
     with (
-        write_catalog(tmp_path / filename_gfits, ext="CATALOG") as out,
+        glass.write_catalog(tmp_path / filename_gfits, ext="CATALOG") as out,
         fitsio.FITS(tmp_path / filename_tfits, "rw", clobber=True) as my_fits,
     ):
         for i in range(my_max):
@@ -99,7 +99,7 @@ def raise_error(msg: str) -> None:
         raise TestWriteError(msg)
 
     try:
-        with write_catalog(tmp_path / filename, ext="CATALOG") as out:
+        with glass.write_catalog(tmp_path / filename, ext="CATALOG") as out:
             for i in range(my_max):
                 if i == except_int:
                     msg = "Unhandled exception"