diff --git a/.gitignore b/.gitignore
index 93b4e32e..856bbb0e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,3 +9,4 @@ dist
 .env
 .coverage*
 coverage*
+.ipynb_checkpoints
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 85a19922..ae311683 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -69,5 +69,6 @@ repos:
     rev: v1.13.0
     hooks:
       - id: mypy
+        files: ^glass/
         additional_dependencies:
           - numpy
diff --git a/docs/conf.py b/docs/conf.py
index 05b15944..874ec590 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -54,7 +54,7 @@
 
 html_logo = "_static/logo.png"
 html_favicon = "_static/favicon.ico"
-html_css_files = []  # type: ignore[var-annotated]
+html_css_files: list[str] = []
 
 
 # -- Intersphinx -------------------------------------------------------------
diff --git a/glass/core/algorithm.py b/glass/core/algorithm.py
index f3d2d0a7..4cd6b953 100644
--- a/glass/core/algorithm.py
+++ b/glass/core/algorithm.py
@@ -7,12 +7,12 @@
 
 
 def nnls(
-    a: npt.ArrayLike,
-    b: npt.ArrayLike,
+    a: npt.NDArray[np.float64],
+    b: npt.NDArray[np.float64],
     *,
     tol: float = 0.0,
     maxiter: int | None = None,
-) -> npt.ArrayLike:
+) -> npt.NDArray[np.float64]:
     """
     Compute a non-negative least squares solution.
 
diff --git a/glass/core/array.py b/glass/core/array.py
index a30e3fb7..c963c35b 100644
--- a/glass/core/array.py
+++ b/glass/core/array.py
@@ -1,18 +1,35 @@
 """Module for array utilities."""
 
+from __future__ import annotations
+
+import typing
 from functools import partial
 
 import numpy as np
+import numpy.typing as npt
+
+if typing.TYPE_CHECKING:
+    import collections.abc
 
 
-def broadcast_first(*arrays):  # type: ignore[no-untyped-def]
+def broadcast_first(
+    *arrays: npt.NDArray[np.float64],
+) -> tuple[npt.NDArray[np.float64], ...]:
     """Broadcast arrays, treating the first axis as common."""
     arrays = tuple(np.moveaxis(a, 0, -1) if np.ndim(a) else a for a in arrays)
     arrays = np.broadcast_arrays(*arrays)
     return tuple(np.moveaxis(a, -1, 0) if np.ndim(a) else a for a in arrays)
 
 
-def broadcast_leading_axes(*args):  # type: ignore[no-untyped-def]
+def broadcast_leading_axes(
+    *args: tuple[
+        float | npt.NDArray[np.float64],
+        int,
+    ],
+) -> tuple[
+    tuple[int, ...],
+    typing.Unpack[tuple[npt.NDArray[np.float64], ...]],
+]:
     """
     Broadcast all but the last N axes.
 
@@ -49,7 +66,15 @@ def broadcast_leading_axes(*args):  # type: ignore[no-untyped-def]
     return (dims, *arrs)
 
 
-def ndinterp(x, xp, fp, axis=-1, left=None, right=None, period=None):  # type: ignore[no-untyped-def] # noqa: PLR0913
+def ndinterp(  # noqa: PLR0913
+    x: float | npt.NDArray[np.float64],
+    xp: collections.abc.Sequence[float] | npt.NDArray[np.float64],
+    fp: collections.abc.Sequence[float] | npt.NDArray[np.float64],
+    axis: int = -1,
+    left: float | None = None,
+    right: float | None = None,
+    period: float | None = None,
+) -> npt.NDArray[np.float64]:
     """Interpolate multi-dimensional array over axis."""
     return np.apply_along_axis(
         partial(np.interp, x, xp),
@@ -61,8 +86,16 @@ def ndinterp(x, xp, fp, axis=-1, left=None, right=None, period=None):  # type: i
     )
 
 
-def trapz_product(f, *ff, axis=-1):  # type: ignore[no-untyped-def]
+def trapz_product(
+    f: tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]],
+    *ff: tuple[
+        npt.NDArray[np.float64],
+        npt.NDArray[np.float64],
+    ],
+    axis: int = -1,
+) -> npt.NDArray[np.float64]:
     """Trapezoidal rule for a product of functions."""
+    x: npt.NDArray[np.float64]
     x, _ = f
     for x_, _ in ff:
         x = np.union1d(
@@ -72,10 +105,19 @@ def trapz_product(f, *ff, axis=-1):  # type: ignore[no-untyped-def]
     y = np.interp(x, *f)
     for f_ in ff:
         y *= np.interp(x, *f_)
-    return np.trapz(y, x, axis=axis)  # type: ignore[attr-defined]
+    return np.trapz(  # type: ignore[attr-defined, no-any-return]
+        y,
+        x,
+        axis=axis,
+    )
 
 
-def cumtrapz(f, x, dtype=None, out=None):  # type: ignore[no-untyped-def]
+def cumtrapz(
+    f: npt.NDArray[np.int_] | npt.NDArray[np.float64],
+    x: npt.NDArray[np.int_] | npt.NDArray[np.float64],
+    dtype: npt.DTypeLike | None = None,
+    out: npt.NDArray[np.float64] | None = None,
+) -> npt.NDArray[np.float64]:
     """Cumulative trapezoidal rule along last axis."""
     if out is None:
         out = np.empty_like(f, dtype=dtype)
diff --git a/glass/ext/__init__.py b/glass/ext/__init__.py
index eba224bb..2b55ab0f 100644
--- a/glass/ext/__init__.py
+++ b/glass/ext/__init__.py
@@ -7,7 +7,7 @@
 """
 
 
-def _extend_path(path, name) -> list:  # type: ignore[no-untyped-def, type-arg]
+def _extend_path(path: list[str], name: str) -> list[str]:
     import os.path
     from pkgutil import extend_path
 
diff --git a/glass/fields.py b/glass/fields.py
index 0994b1eb..613a939e 100644
--- a/glass/fields.py
+++ b/glass/fields.py
@@ -35,32 +35,20 @@
 import numpy.typing as npt
 from gaussiancl import gaussiancl
 
-# types
-Size = typing.Optional[typing.Union[int, tuple[int, ...]]]
-Iternorm = tuple[typing.Optional[int], npt.NDArray[typing.Any], npt.NDArray[typing.Any]]
-ClTransform = typing.Union[
-    str,
-    typing.Callable[[npt.NDArray[typing.Any]], npt.NDArray[typing.Any]],
-]
 Cls = collections.abc.Sequence[
-    typing.Union[npt.NDArray[typing.Any], collections.abc.Sequence[float]]
+    typing.Union[npt.NDArray[np.float64], collections.abc.Sequence[float]]
 ]
-Alms = npt.NDArray[typing.Any]
 
 
 def iternorm(
     k: int,
-    cov: collections.abc.Iterable[npt.NDArray[typing.Any]],
-    size: Size = None,
-) -> collections.abc.Generator[Iternorm, None, None]:
+    cov: collections.abc.Iterable[npt.NDArray[np.float64]],
+    size: int | tuple[int, ...] = (),
+) -> collections.abc.Generator[
+    tuple[int | None, npt.NDArray[np.float64], npt.NDArray[np.float64]]
+]:
     """Return the vector a and variance sigma^2 for iterative normal sampling."""
-    n: tuple[int, ...]
-    if size is None:
-        n = ()
-    elif isinstance(size, int):
-        n = (size,)
-    else:
-        n = size
+    n = (size,) if isinstance(size, int) else size
 
     m = np.zeros((*n, k, k))
     a = np.zeros((*n, k))
@@ -115,27 +103,29 @@ def cls2cov(
     nl: int,
     nf: int,
     nc: int,
-) -> collections.abc.Generator[npt.NDArray[typing.Any], None, None]:
+) -> collections.abc.Generator[npt.NDArray[np.float64]]:
     """Return array of cls as a covariance matrix for iterative sampling."""
     cov = np.zeros((nl, nc + 1))
     end = 0
     for j in range(nf):
         begin, end = end, end + j + 1
         for i, cl in enumerate(cls[begin:end][: nc + 1]):
-            if cl is None:
-                cov[:, i] = 0  # type: ignore[unreachable]
-            else:
-                if i == 0 and np.any(np.less(cl, 0)):
-                    msg = "negative values in cl"
-                    raise ValueError(msg)
-                n = len(cl)
-                cov[:n, i] = cl
-                cov[n:, i] = 0
+            if i == 0 and np.any(np.less(cl, 0)):
+                msg = "negative values in cl"
+                raise ValueError(msg)
+            n = len(cl)
+            cov[:n, i] = cl
+            cov[n:, i] = 0
         cov /= 2
         yield cov
 
 
-def multalm(alm: Alms, bl: npt.NDArray[typing.Any], *, inplace: bool = False) -> Alms:
+def multalm(
+    alm: npt.NDArray[np.complex128],
+    bl: npt.NDArray[np.float64],
+    *,
+    inplace: bool = False,
+) -> npt.NDArray[np.complex128]:
     """Multiply alm by bl."""
     n = len(bl)
     out = np.asanyarray(alm) if inplace else np.copy(alm)
@@ -144,11 +134,15 @@ def multalm(alm: Alms, bl: npt.NDArray[typing.Any], *, inplace: bool = False) ->
     return out
 
 
-def transform_cls(cls: Cls, tfm: ClTransform, pars: tuple[typing.Any, ...] = ()) -> Cls:
+def transform_cls(
+    cls: Cls,
+    tfm: str | typing.Callable[[npt.NDArray[np.float64]], npt.NDArray[np.float64]],
+    pars: tuple[typing.Any, ...] = (),
+) -> Cls:
     """Transform Cls to Gaussian Cls."""
     gls = []
     for cl in cls:
-        if cl is not None and len(cl) > 0:  # type: ignore[redundant-expr]
+        if len(cl) > 0:
             monopole = 0.0 if cl[0] == 0 else None
             gl, info, _, _ = gaussiancl(cl, tfm, pars, monopole=monopole)
             if info == 0:
@@ -194,7 +188,7 @@ def discretized_cls(
 
     gls = []
     for cl in cls:
-        if cl is not None and len(cl) > 0:  # type: ignore[redundant-expr]
+        if len(cl) > 0:
             if lmax is not None:
                 cl = cl[: lmax + 1]  # noqa: PLW2901
             if nside is not None:
@@ -218,7 +212,7 @@ def generate_gaussian(
     *,
     ncorr: int | None = None,
     rng: np.random.Generator | None = None,
-) -> collections.abc.Generator[npt.NDArray[typing.Any], None, None]:
+) -> collections.abc.Generator[npt.NDArray[np.float64]]:
     """
     Sample Gaussian random fields from Cls iteratively.
 
@@ -255,7 +249,7 @@ def generate_gaussian(
         ncorr = ngrf - 1
 
     # number of modes
-    n = max((len(gl) for gl in gls if gl is not None), default=0)  # type: ignore[redundant-expr]
+    n = max((len(gl) for gl in gls), default=0)
     if n == 0:
         msg = "all gls are empty"
         raise ValueError(msg)
@@ -304,7 +298,7 @@ def generate_lognormal(
     *,
     ncorr: int | None = None,
     rng: np.random.Generator | None = None,
-) -> collections.abc.Generator[npt.NDArray[typing.Any], None, None]:
+) -> collections.abc.Generator[npt.NDArray[np.float64]]:
     """Sample lognormal random fields from Gaussian Cls iteratively."""
     for i, m in enumerate(generate_gaussian(gls, nside, ncorr=ncorr, rng=rng)):
         # compute the variance of the auto-correlation
@@ -326,7 +320,14 @@ def generate_lognormal(
         yield m
 
 
-def getcl(cls, i, j, lmax=None):  # type: ignore[no-untyped-def]
+def getcl(
+    cls: collections.abc.Sequence[
+        npt.NDArray[np.float64] | collections.abc.Sequence[float]
+    ],
+    i: int,
+    j: int,
+    lmax: int | None = None,
+) -> npt.NDArray[np.float64] | collections.abc.Sequence[float]:
     """
     Return a specific angular power spectrum from an array.
 
@@ -356,12 +357,14 @@ def getcl(cls, i, j, lmax=None):  # type: ignore[no-untyped-def]
     return cl
 
 
-def effective_cls(  # type: ignore[no-untyped-def]
-    cls,
-    weights1,
-    weights2=None,
+def effective_cls(
+    cls: collections.abc.Sequence[
+        npt.NDArray[np.float64] | collections.abc.Sequence[float]
+    ],
+    weights1: npt.NDArray[np.float64],
+    weights2: npt.NDArray[np.float64] | None = None,
     *,
-    lmax=None,
+    lmax: int | None = None,
 ) -> npt.NDArray[np.float64]:
     r"""
     Compute effective angular power spectra from weights.
@@ -411,10 +414,11 @@ def effective_cls(  # type: ignore[no-untyped-def]
 
     # get the iterator over leading weight axes
     # auto-spectra do not repeat identical computations
-    if weights2 is weights1:
-        pairs = combinations_with_replacement(np.ndindex(shape1[1:]), 2)
-    else:
-        pairs = product(np.ndindex(shape1[1:]), np.ndindex(shape2[1:]))  # type: ignore[assignment]
+    pairs = (
+        combinations_with_replacement(np.ndindex(shape1[1:]), 2)
+        if weights2 is weights1
+        else product(np.ndindex(shape1[1:]), np.ndindex(shape2[1:]))
+    )
 
     # create the output array: axes for all input axes plus lmax+1
     out = np.empty(shape1[1:] + shape2[1:] + (lmax + 1,))
@@ -427,7 +431,7 @@ def effective_cls(  # type: ignore[no-untyped-def]
     for j1, j2 in pairs:
         w1, w2 = weights1[c + j1], weights2[c + j2]
         cl = sum(
-            w1[i1] * w2[i2] * getcl(cls, i1, i2, lmax=lmax)  # type: ignore[no-untyped-call]
+            w1[i1] * w2[i2] * getcl(cls, i1, i2, lmax=lmax)
             for i1, i2 in np.ndindex(n, n)
         )
         out[j1 + j2] = cl
diff --git a/glass/galaxies.py b/glass/galaxies.py
index 6757435b..16c26c9a 100644
--- a/glass/galaxies.py
+++ b/glass/galaxies.py
@@ -39,11 +39,11 @@
 
 
 def redshifts(
-    n: int | npt.ArrayLike,
+    n: int | npt.NDArray[np.float64],
     w: RadialWindow,
     *,
     rng: np.random.Generator | None = None,
-) -> npt.NDArray[typing.Any]:
+) -> npt.NDArray[np.float64]:
     """
     Sample redshifts from a radial window function.
 
@@ -67,13 +67,13 @@ def redshifts(
 
 
 def redshifts_from_nz(
-    count: int | npt.ArrayLike,
-    z: npt.ArrayLike,
-    nz: npt.ArrayLike,
+    count: int | npt.NDArray[np.float64],
+    z: npt.NDArray[np.float64],
+    nz: npt.NDArray[np.float64],
     *,
     rng: np.random.Generator | None = None,
     warn: bool = True,
-) -> npt.NDArray[typing.Any]:
+) -> npt.NDArray[np.float64]:
     """
     Generate galaxy redshifts from a source distribution.
 
@@ -117,10 +117,11 @@ def redshifts_from_nz(
         rng = np.random.default_rng()
 
     # bring inputs' leading axes into common shape
-    dims, count, z, nz = broadcast_leading_axes((count, 0), (z, 1), (nz, 1))  # type: ignore[no-untyped-call]
+    dims, *rest = broadcast_leading_axes((count, 0), (z, 1), (nz, 1))
+    count_out, z_out, nz_out = rest
 
     # list of results for all dimensions
-    redshifts = np.empty(count.sum())  # type: ignore[union-attr]
+    redshifts = np.empty(count_out.sum())
 
     # keep track of the number of sampled redshifts
     total = 0
@@ -128,16 +129,16 @@ 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 = cumtrapz(nz[k], z[k], dtype=float)  # type: ignore[call-overload, index, no-untyped-call]
+        cdf = cumtrapz(nz_out[k], z_out[k], dtype=float)
         cdf /= cdf[-1]
 
         # sample redshifts and store result
-        redshifts[total : total + count[k]] = np.interp(  # type: ignore[call-overload, index, misc, operator]
-            rng.uniform(0, 1, size=count[k]),  # type: ignore[arg-type, call-overload, index]
+        redshifts[total : total + count_out[k]] = np.interp(
+            rng.uniform(0, 1, size=count_out[k]),
             cdf,
-            z[k],  # type: ignore[arg-type, call-overload, index]
+            z_out[k],
         )
-        total += count[k]  # type: ignore[assignment, call-overload, index, operator]
+        total += count_out[k]
 
     assert total == redshifts.size  # noqa: S101
 
@@ -145,15 +146,15 @@ def redshifts_from_nz(
 
 
 def galaxy_shear(  # noqa: PLR0913
-    lon: npt.NDArray[typing.Any],
-    lat: npt.NDArray[typing.Any],
-    eps: npt.NDArray[typing.Any],
-    kappa: npt.NDArray[typing.Any],
-    gamma1: npt.NDArray[typing.Any],
-    gamma2: npt.NDArray[typing.Any],
+    lon: npt.NDArray[np.float64],
+    lat: npt.NDArray[np.float64],
+    eps: npt.NDArray[np.float64],
+    kappa: npt.NDArray[np.float64],
+    gamma1: npt.NDArray[np.float64],
+    gamma2: npt.NDArray[np.float64],
     *,
     reduced_shear: bool = True,
-) -> npt.NDArray[typing.Any]:
+) -> npt.NDArray[np.float64]:
     """
     Observed galaxy shears from weak lensing.
 
@@ -212,13 +213,13 @@ def galaxy_shear(  # noqa: PLR0913
 
 
 def gaussian_phz(
-    z: npt.ArrayLike,
-    sigma_0: float | npt.ArrayLike,
+    z: float | npt.NDArray[np.float64],
+    sigma_0: float | npt.NDArray[np.float64],
     *,
-    lower: npt.ArrayLike | None = None,
-    upper: npt.ArrayLike | None = None,
+    lower: float | npt.NDArray[np.float64] | None = None,
+    upper: float | npt.NDArray[np.float64] | None = None,
     rng: np.random.Generator | None = None,
-) -> npt.NDArray[typing.Any]:
+) -> float | npt.NDArray[np.float64]:
     r"""
     Photometric redshifts assuming a Gaussian error.
 
@@ -270,33 +271,33 @@ def gaussian_phz(
     sigma = np.add(1, z) * sigma_0
     dims = np.shape(sigma)
 
-    zphot = rng.normal(z, sigma)  # type: ignore[arg-type]
+    zphot = rng.normal(z, sigma)
 
     if lower is None:
         lower = 0.0
     if upper is None:
         upper = np.inf
 
-    if not np.all(lower < upper):  # type: ignore[operator]
+    if not np.all(lower < upper):
         msg = "requires lower < upper"
         raise ValueError(msg)
 
     if not dims:
-        while zphot < lower or zphot > upper:  # type: ignore[operator]
-            zphot = rng.normal(z, sigma)  # type: ignore[arg-type]
+        while zphot < lower or zphot > upper:
+            zphot = rng.normal(z, sigma)
     else:
         z = np.broadcast_to(z, dims)
-        trunc = np.where((zphot < lower) | (zphot > upper))[0]  # type: ignore[operator]
+        trunc = np.where((zphot < lower) | (zphot > upper))[0]
         while trunc.size:
-            znew = rng.normal(z[trunc], sigma[trunc])  # type: ignore[arg-type]
-            zphot[trunc] = znew  # type: ignore[index]
-            trunc = trunc[(znew < lower) | (znew > upper)]  # type: ignore[operator]
+            znew = rng.normal(z[trunc], sigma[trunc])
+            zphot[trunc] = znew
+            trunc = trunc[(znew < lower) | (znew > upper)]
 
-    return zphot  # type: ignore[return-value]
+    return zphot
 
 
 def kappa_ia_nla(  # noqa: PLR0913
-    delta: npt.NDArray[typing.Any],
+    delta: npt.NDArray[np.float64],
     zeff: float,
     a_ia: float,
     cosmo: Cosmology,
@@ -306,7 +307,7 @@ def kappa_ia_nla(  # noqa: PLR0913
     lbar: float = 0.0,
     l0: float = 1e-9,
     beta: float = 0.0,
-) -> npt.NDArray[typing.Any]:
+) -> npt.NDArray[np.float64]:
     r"""
     Effective convergence from intrinsic alignments using the NLA model.
 
diff --git a/glass/lensing.py b/glass/lensing.py
index b644f752..4d8b19e6 100644
--- a/glass/lensing.py
+++ b/glass/lensing.py
@@ -46,14 +46,14 @@
 
 
 def from_convergence(  # noqa: PLR0913
-    kappa: npt.NDArray[typing.Any],
+    kappa: npt.NDArray[np.float64],
     lmax: int | None = None,
     *,
     potential: bool = False,
     deflection: bool = False,
     shear: bool = False,
     discretized: bool = True,
-) -> tuple[npt.NDArray[typing.Any], ...]:
+) -> tuple[npt.NDArray[np.float64], ...]:
     r"""
     Compute other weak lensing maps from the convergence.
 
@@ -171,7 +171,7 @@ def from_convergence(  # noqa: PLR0913
     l = np.arange(lmax + 1)  # noqa: E741
 
     # this tuple will be returned
-    results = ()
+    results: tuple[npt.NDArray[np.float64], ...] = ()
 
     # convert convergence to potential
     fl = np.divide(-2, l * (l + 1), where=(l > 0), out=np.zeros(lmax + 1))
@@ -180,7 +180,7 @@ def from_convergence(  # noqa: PLR0913
     # if potential is requested, compute map and add to output
     if potential:
         psi = hp.alm2map(alm, nside, lmax=lmax)
-        results += (psi,)  # type: ignore[assignment]
+        results += (psi,)
 
     # if no spin-weighted maps are requested, stop here
     if not (deflection or shear):
@@ -199,7 +199,7 @@ def from_convergence(  # noqa: PLR0913
     if deflection:
         alpha = hp.alm2map_spin([alm, blm], nside, 1, lmax)
         alpha = alpha[0] + 1j * alpha[1]
-        results += (alpha,)  # type: ignore[assignment]
+        results += (alpha,)
 
     # if no shear is requested, stop here
     if not shear:
@@ -217,18 +217,18 @@ def from_convergence(  # noqa: PLR0913
     # transform to shear maps
     gamma = hp.alm2map_spin([alm, blm], nside, 2, lmax)
     gamma = gamma[0] + 1j * gamma[1]
-    results += (gamma,)  # type: ignore[assignment]
+    results += (gamma,)
 
     # all done
     return results
 
 
 def shear_from_convergence(
-    kappa: npt.NDArray[typing.Any],
+    kappa: npt.NDArray[np.float64],
     lmax: int | None = None,
     *,
     discretized: bool = True,
-) -> npt.NDArray[typing.Any]:
+) -> npt.NDArray[np.float64]:
     r"""
     Weak lensing shear from convergence.
 
@@ -280,11 +280,11 @@ def __init__(self, cosmo: Cosmology) -> None:
         self.x3: float = 0.0
         self.w3: float = 0.0
         self.r23: float = 1.0
-        self.delta3: npt.NDArray[typing.Any] = np.array(0.0)
-        self.kappa2: npt.NDArray[typing.Any] | None = None
-        self.kappa3: npt.NDArray[typing.Any] | None = None
+        self.delta3: npt.NDArray[np.float64] = np.array(0.0)
+        self.kappa2: npt.NDArray[np.float64] | None = None
+        self.kappa3: npt.NDArray[np.float64] | None = None
 
-    def add_window(self, delta: npt.NDArray[typing.Any], w: RadialWindow) -> None:
+    def add_window(self, delta: npt.NDArray[np.float64], w: RadialWindow) -> None:
         """
         Add a mass plane from a window function to the convergence.
 
@@ -293,13 +293,23 @@ def add_window(self, delta: npt.NDArray[typing.Any], w: RadialWindow) -> None:
 
         """
         zsrc = w.zeff
-        lens_weight = np.trapz(w.wa, w.za) / np.interp(zsrc, w.za, w.wa)  # type: ignore[attr-defined]
+        lens_weight = float(
+            np.trapz(  # type: ignore[attr-defined]
+                w.wa,
+                w.za,
+            )
+            / np.interp(
+                zsrc,
+                w.za,
+                w.wa,
+            )
+        )
 
         self.add_plane(delta, zsrc, lens_weight)
 
     def add_plane(
         self,
-        delta: npt.NDArray[typing.Any],
+        delta: npt.NDArray[np.float64],
         zsrc: float,
         wlens: float = 1.0,
     ) -> None:
@@ -351,12 +361,12 @@ def zsrc(self) -> float:
         return self.z3
 
     @property
-    def kappa(self) -> npt.NDArray[typing.Any] | None:
+    def kappa(self) -> npt.NDArray[np.float64] | None:
         """The current convergence plane."""
         return self.kappa3
 
     @property
-    def delta(self) -> npt.NDArray[typing.Any]:
+    def delta(self) -> npt.NDArray[np.float64]:
         """The current matter plane."""
         return self.delta3
 
@@ -369,7 +379,7 @@ def wlens(self) -> float:
 def multi_plane_matrix(
     shells: collections.abc.Sequence[RadialWindow],
     cosmo: Cosmology,
-) -> npt.ArrayLike:
+) -> npt.NDArray[np.float64]:
     """Compute the matrix of lensing contributions from each shell."""
     mpc = MultiPlaneConvergence(cosmo)
     wmat = np.eye(len(shells))
@@ -380,10 +390,10 @@ def multi_plane_matrix(
 
 
 def multi_plane_weights(
-    weights: npt.ArrayLike,
+    weights: npt.NDArray[np.float64],
     shells: collections.abc.Sequence[RadialWindow],
     cosmo: Cosmology,
-) -> npt.ArrayLike:
+) -> npt.NDArray[np.float64]:
     """
     Compute effective weights for multi-plane convergence.
 
@@ -416,14 +426,17 @@ def multi_plane_weights(
     weights = weights / np.sum(weights, axis=0)
     # combine weights and the matrix of lensing contributions
     mat = multi_plane_matrix(shells, cosmo)
-    return np.matmul(mat.T, weights)  # type: ignore[no-any-return, union-attr]
+    return np.matmul(mat.T, weights)  # type: ignore[no-any-return]
 
 
 def deflect(
-    lon: npt.ArrayLike,
-    lat: npt.ArrayLike,
-    alpha: npt.ArrayLike,
-) -> npt.NDArray[typing.Any]:
+    lon: float | npt.NDArray[np.float64],
+    lat: float | npt.NDArray[np.float64],
+    alpha: complex | list[float] | npt.NDArray[np.complex128] | npt.NDArray[np.float64],
+) -> tuple[
+    npt.NDArray[np.float64],
+    npt.NDArray[np.float64],
+]:
     r"""
     Apply deflections to positions.
 
@@ -477,4 +490,4 @@ def deflect(
 
     d = np.arctan2(sa * sg, st * ca - ct * sa * cg)
 
-    return lon - np.degrees(d), np.degrees(tp)  # type: ignore[return-value]
+    return lon - np.degrees(d), np.degrees(tp)
diff --git a/glass/observations.py b/glass/observations.py
index 5e98b146..bd09a92f 100644
--- a/glass/observations.py
+++ b/glass/observations.py
@@ -29,7 +29,6 @@
 from __future__ import annotations
 
 import math
-import typing
 
 import healpy as hp
 import numpy as np
@@ -42,7 +41,7 @@ def vmap_galactic_ecliptic(
     nside: int,
     galactic: tuple[float, float] = (30, 90),
     ecliptic: tuple[float, float] = (20, 80),
-) -> npt.NDArray[typing.Any]:
+) -> npt.NDArray[np.float64]:
     """
     Visibility map masking galactic and ecliptic plane.
 
@@ -82,12 +81,12 @@ def vmap_galactic_ecliptic(
 
 
 def gaussian_nz(
-    z: npt.NDArray[typing.Any],
-    mean: npt.ArrayLike,
-    sigma: npt.ArrayLike,
+    z: npt.NDArray[np.float64],
+    mean: float | npt.NDArray[np.float64],
+    sigma: float | npt.NDArray[np.float64],
     *,
-    norm: npt.ArrayLike | None = None,
-) -> npt.NDArray[typing.Any]:
+    norm: npt.NDArray[np.float64] | None = None,
+) -> npt.NDArray[np.float64]:
     r"""
     Gaussian redshift distribution.
 
@@ -115,7 +114,11 @@ def gaussian_nz(
     sigma = np.reshape(sigma, np.shape(sigma) + (1,) * np.ndim(z))
 
     nz = np.exp(-(((z - mean) / sigma) ** 2) / 2)
-    nz /= np.trapz(nz, z, axis=-1)[..., np.newaxis]  # type: ignore[attr-defined]
+    nz /= np.trapz(  # type: ignore[attr-defined]
+        nz,
+        z,
+        axis=-1,
+    )[..., np.newaxis]
 
     if norm is not None:
         nz *= norm
@@ -124,13 +127,13 @@ def gaussian_nz(
 
 
 def smail_nz(
-    z: npt.NDArray[typing.Any],
-    z_mode: npt.ArrayLike,
-    alpha: npt.ArrayLike,
-    beta: npt.ArrayLike,
+    z: npt.NDArray[np.float64],
+    z_mode: npt.NDArray[np.float64],
+    alpha: npt.NDArray[np.float64],
+    beta: npt.NDArray[np.float64],
     *,
-    norm: npt.ArrayLike | None = None,
-) -> npt.NDArray[typing.Any]:
+    norm: npt.NDArray[np.float64] | None = None,
+) -> npt.NDArray[np.float64]:
     r"""
     Redshift distribution following Smail et al. (1994).
 
@@ -174,7 +177,11 @@ def smail_nz(
     beta = np.asanyarray(beta)[..., np.newaxis]
 
     pz = z**alpha * np.exp(-alpha / beta * (z / z_mode) ** beta)
-    pz /= np.trapz(pz, z, axis=-1)[..., np.newaxis]  # type: ignore[attr-defined]
+    pz /= np.trapz(  # type: ignore[attr-defined]
+        pz,
+        z,
+        axis=-1,
+    )[..., np.newaxis]
 
     if norm is not None:
         pz *= norm
@@ -221,8 +228,8 @@ def fixed_zbins(
 
 
 def equal_dens_zbins(
-    z: npt.NDArray[typing.Any],
-    nz: npt.NDArray[typing.Any],
+    z: npt.NDArray[np.float64],
+    nz: npt.NDArray[np.float64],
     nbins: int,
 ) -> list[tuple[float, float]]:
     """
@@ -247,7 +254,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 = cumtrapz(nz, z)  # type: ignore[no-untyped-call]
+    cuml_nz = cumtrapz(nz, z)
     cuml_nz /= cuml_nz[[-1]]
     zbinedges = np.interp(np.linspace(0, 1, nbins + 1), cuml_nz, z)
 
@@ -255,11 +262,11 @@ def equal_dens_zbins(
 
 
 def tomo_nz_gausserr(
-    z: npt.NDArray[typing.Any],
-    nz: npt.NDArray[typing.Any],
+    z: npt.NDArray[np.float64],
+    nz: npt.NDArray[np.float64],
     sigma_0: float,
     zbins: list[tuple[float, float]],
-) -> npt.NDArray[typing.Any]:
+) -> npt.NDArray[np.float64]:
     """
     Tomographic redshift bins with a Gaussian redshift error.
 
diff --git a/glass/points.py b/glass/points.py
index 44ee3e05..7c0a5c2b 100644
--- a/glass/points.py
+++ b/glass/points.py
@@ -29,15 +29,30 @@
 
 """  # noqa: D205, D400
 
+from __future__ import annotations
+
+import typing
+
 import healpix
 import numpy as np
+import numpy.typing as npt
 
 from glass.core.array import broadcast_first, broadcast_leading_axes, trapz_product
 
+if typing.TYPE_CHECKING:
+    import collections.abc
+
+    from glass.shells import RadialWindow
+
+
 ARCMIN2_SPHERE = 60**6 // 100 / np.pi
 
 
-def effective_bias(z, bz, w):  # type: ignore[no-untyped-def]
+def effective_bias(
+    z: npt.NDArray[np.float64],
+    bz: npt.NDArray[np.float64],
+    w: RadialWindow,
+) -> npt.NDArray[np.float64]:
     r"""
     Effective bias parameter from a redshift-dependent bias function.
 
@@ -67,16 +82,25 @@ def effective_bias(z, bz, w):  # type: ignore[no-untyped-def]
         \;.
 
     """
-    norm = np.trapz(w.wa, w.za)  # type: ignore[attr-defined]
-    return trapz_product((z, bz), (w.za, w.wa)) / norm  # type: ignore[no-untyped-call]
+    norm = np.trapz(  # type: ignore[attr-defined]
+        w.wa,
+        w.za,
+    )
+    return trapz_product((z, bz), (w.za, w.wa)) / norm  # type: ignore[no-any-return]
 
 
-def linear_bias(delta, b):  # type: ignore[no-untyped-def]
+def linear_bias(
+    delta: npt.NDArray[np.float64],
+    b: float | npt.NDArray[np.float64],
+) -> npt.NDArray[np.float64]:
     r"""Linear bias model :math:`\delta_g = b \, \delta`."""
     return b * delta
 
 
-def loglinear_bias(delta, b):  # type: ignore[no-untyped-def]
+def loglinear_bias(
+    delta: npt.NDArray[np.float64],
+    b: float | npt.NDArray[np.float64],
+) -> npt.NDArray[np.float64]:
     r"""log-linear bias model :math:`\ln(1 + \delta_g) = b \ln(1 + \delta)`."""
     delta_g = np.log1p(delta)
     delta_g *= b
@@ -84,17 +108,23 @@ def loglinear_bias(delta, b):  # type: ignore[no-untyped-def]
     return delta_g
 
 
-def positions_from_delta(  # type: ignore[no-untyped-def] # noqa: PLR0912, PLR0913, PLR0915
-    ngal,
-    delta,
-    bias=None,
-    vis=None,
+def positions_from_delta(  # noqa: PLR0912, PLR0913, PLR0915
+    ngal: float | npt.NDArray[np.float64],
+    delta: npt.NDArray[np.float64],
+    bias: float | npt.NDArray[np.float64] | None = None,
+    vis: npt.NDArray[np.float64] | None = None,
     *,
-    bias_model="linear",
-    remove_monopole=False,
-    batch=1_000_000,
-    rng=None,
-):
+    bias_model: str | typing.Callable[..., typing.Any] = "linear",
+    remove_monopole: bool = False,
+    batch: int = 1_000_000,
+    rng: np.random.Generator | None = None,
+) -> collections.abc.Generator[
+    tuple[
+        npt.NDArray[np.float64],
+        npt.NDArray[np.float64],
+        int | npt.NDArray[np.int_],
+    ]
+]:
     """
     Generate positions tracing a density contrast.
 
@@ -157,18 +187,20 @@ def positions_from_delta(  # type: ignore[no-untyped-def] # noqa: PLR0912, PLR09
 
     # get the bias model
     if isinstance(bias_model, str):
-        bias_model = globals()[f"{bias_model}_bias"]
+        bias_model_callable = globals()[f"{bias_model}_bias"]
     elif not callable(bias_model):
-        msg = "bias_model must be string or callable"
-        raise TypeError(msg)
+        raise TypeError("bias_model must be string or callable")  # noqa: EM101,TRY003
+    else:
+        bias_model_callable = bias_model
 
     # broadcast inputs to common shape of extra dimensions
-    inputs = [(ngal, 0), (delta, 1)]
+    inputs: list[tuple[float | npt.NDArray[np.float64], int]] = [(ngal, 0), (delta, 1)]
     if bias is not None:
-        inputs += [(bias, 0)]
+        inputs.append((bias, 0))
     if vis is not None:
-        inputs += [(vis, 1)]
-    dims, ngal, delta, *rest = broadcast_leading_axes(*inputs)  # type: ignore[no-untyped-call]
+        inputs.append((vis, 1))
+    dims, *rest = broadcast_leading_axes(*inputs)
+    ngal, delta, *rest = rest
     if bias is not None:
         bias, *rest = rest
     if vis is not None:
@@ -177,7 +209,11 @@ def positions_from_delta(  # type: ignore[no-untyped-def] # noqa: PLR0912, PLR09
     # iterate the leading dimensions
     for k in np.ndindex(dims):
         # compute density contrast from bias model, or copy
-        n = np.copy(delta[k]) if bias is None else bias_model(delta[k], bias[k])
+        n = (
+            np.copy(delta[k])
+            if bias is None
+            else bias_model_callable(delta[k], bias[k])
+        )
 
         # remove monopole if asked to
         if remove_monopole:
@@ -208,11 +244,12 @@ def positions_from_delta(  # type: ignore[no-untyped-def] # noqa: PLR0912, PLR09
         nside = healpix.npix2nside(npix)
 
         # create a mask to report the count in the right axis
+        cmask: int | npt.NDArray[np.int_]
         if dims:
             cmask = np.zeros(dims, dtype=int)
             cmask[k] = 1
         else:
-            cmask = 1  # type: ignore[assignment]
+            cmask = 1
 
         # sample the map in batches
         step = 1000
@@ -227,7 +264,7 @@ def positions_from_delta(  # type: ignore[no-untyped-def] # noqa: PLR0912, PLR09
                 size += q[-1]
             else:
                 # how many pixels from this group do we need?
-                stop += np.searchsorted(q, batch - size, side="right")
+                stop += int(np.searchsorted(q, batch - size, side="right"))
                 # if the first pixel alone is too much, use it anyway
                 if stop == start:
                     stop += 1
@@ -245,7 +282,17 @@ def positions_from_delta(  # type: ignore[no-untyped-def] # noqa: PLR0912, PLR09
         assert np.sum(n[stop:]) == 0  # noqa: S101
 
 
-def uniform_positions(ngal, *, rng=None):  # type: ignore[no-untyped-def]
+def uniform_positions(
+    ngal: float | npt.NDArray[np.int_] | npt.NDArray[np.float64],
+    *,
+    rng: np.random.Generator | None = None,
+) -> collections.abc.Generator[
+    tuple[
+        npt.NDArray[np.float64],
+        npt.NDArray[np.float64],
+        int | npt.NDArray[np.int_],
+    ]
+]:
     """
     Generate positions uniformly over the sphere.
 
@@ -289,16 +336,20 @@ def uniform_positions(ngal, *, rng=None):  # type: ignore[no-untyped-def]
         lat = np.rad2deg(np.arcsin(rng.uniform(-1, 1, size=ngal[k])))
 
         # report count
+        count: int | npt.NDArray[np.int_]
         if dims:
             count = np.zeros(dims, dtype=int)
             count[k] = ngal[k]
         else:
-            count = int(ngal[k])  # type: ignore[assignment]
+            count = int(ngal[k])
 
         yield lon, lat, count
 
 
-def position_weights(densities, bias=None):  # type: ignore[no-untyped-def]
+def position_weights(
+    densities: npt.NDArray[np.float64],
+    bias: npt.NDArray[np.float64] | None = None,
+) -> npt.NDArray[np.float64]:
     r"""
     Compute relative weights for angular clustering.
 
@@ -323,7 +374,7 @@ def position_weights(densities, bias=None):  # type: ignore[no-untyped-def]
     """
     # bring densities and bias into the same shape
     if bias is not None:
-        densities, bias = broadcast_first(densities, bias)  # type: ignore[no-untyped-call]
+        densities, bias = 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/py.typed b/glass/py.typed
new file mode 100644
index 00000000..e69de29b
diff --git a/glass/shapes.py b/glass/shapes.py
index bf0869cb..25dc6f45 100644
--- a/glass/shapes.py
+++ b/glass/shapes.py
@@ -24,13 +24,17 @@
 
 from __future__ import annotations
 
-import typing
-
 import numpy as np
 import numpy.typing as npt
 
 
-def triaxial_axis_ratio(zeta, xi, size=None, *, rng=None):  # type: ignore[no-untyped-def]
+def triaxial_axis_ratio(
+    zeta: float | npt.NDArray[np.float64],
+    xi: float | npt.NDArray[np.float64],
+    size: int | tuple[int, ...] | None = None,
+    *,
+    rng: np.random.Generator | None = None,
+) -> npt.NDArray[np.float64]:
     r"""
     Axis ratio of a randomly projected triaxial ellipsoid.
 
@@ -87,12 +91,20 @@ def triaxial_axis_ratio(zeta, xi, size=None, *, rng=None):  # type: ignore[no-un
     C = 1 + z2m1 * cos2_phi  # noqa: N806
 
     # eq. (12)
-    return np.sqrt(
+    return np.sqrt(  # type: ignore[no-any-return]
         (A + C - np.sqrt((A - C) ** 2 + B2)) / (A + C + np.sqrt((A - C) ** 2 + B2)),
     )
 
 
-def ellipticity_ryden04(mu, sigma, gamma, sigma_gamma, size=None, *, rng=None):  # type: ignore[no-untyped-def] # noqa: PLR0913
+def ellipticity_ryden04(  # noqa: PLR0913
+    mu: float | npt.NDArray[np.float64],
+    sigma: float | npt.NDArray[np.float64],
+    gamma: float | npt.NDArray[np.float64],
+    sigma_gamma: float | npt.NDArray[np.float64],
+    size: int | tuple[int, ...] | None = None,
+    *,
+    rng: np.random.Generator | None = None,
+) -> npt.NDArray[np.float64]:
     r"""
     Ellipticity distribution following Ryden (2004).
 
@@ -155,22 +167,22 @@ def ellipticity_ryden04(mu, sigma, gamma, sigma_gamma, size=None, *, rng=None):
     xi = (1 - gam) * zeta
 
     # random projection of random triaxial ellipsoid
-    q = triaxial_axis_ratio(zeta, xi, rng=rng)  # type: ignore[no-untyped-call]
+    q = triaxial_axis_ratio(zeta, xi, rng=rng)
 
     # assemble ellipticity with random complex phase
     e = np.exp(1j * rng.uniform(0, 2 * np.pi, size=np.shape(q)))
     e *= (1 - q) / (1 + q)
 
     # return the ellipticity
-    return e
+    return e  # type: ignore[no-any-return]
 
 
 def ellipticity_gaussian(
-    count: int | npt.ArrayLike,
-    sigma: npt.ArrayLike,
+    count: int | npt.NDArray[np.int_],
+    sigma: float | npt.NDArray[np.float64],
     *,
     rng: np.random.Generator | None = None,
-) -> npt.NDArray[typing.Any]:
+) -> npt.NDArray[np.complex128]:
     r"""
     Sample Gaussian galaxy ellipticities.
 
@@ -220,11 +232,11 @@ def ellipticity_gaussian(
 
 
 def ellipticity_intnorm(
-    count: int | npt.ArrayLike,
-    sigma: npt.ArrayLike,
+    count: int | npt.NDArray[np.int_],
+    sigma: float | npt.NDArray[np.float64],
     *,
     rng: np.random.Generator | None = None,
-) -> npt.NDArray[typing.Any]:
+) -> npt.NDArray[np.complex128]:
     r"""
     Sample galaxy ellipticities with intrinsic normal distribution.
 
diff --git a/glass/shells.py b/glass/shells.py
index a1309b95..479b589f 100644
--- a/glass/shells.py
+++ b/glass/shells.py
@@ -56,22 +56,30 @@
 if typing.TYPE_CHECKING:
     from cosmology import Cosmology
 
-# types
-ArrayLike1D = typing.Union[collections.abc.Sequence[float], npt.NDArray[typing.Any]]
-WeightFunc = typing.Callable[[ArrayLike1D], npt.NDArray[typing.Any]]
+ArrayLike1D = typing.Union[collections.abc.Sequence[float], npt.NDArray[np.float64]]
+WeightFunc = typing.Callable[[ArrayLike1D], npt.NDArray[np.float64]]
 
 
-def distance_weight(z: npt.ArrayLike, cosmo: Cosmology) -> npt.NDArray[typing.Any]:
+def distance_weight(
+    z: npt.NDArray[np.float64],
+    cosmo: Cosmology,
+) -> npt.NDArray[np.float64]:
     """Uniform weight in comoving distance."""
     return 1 / cosmo.ef(z)  # type: ignore[no-any-return]
 
 
-def volume_weight(z: npt.ArrayLike, cosmo: Cosmology) -> npt.NDArray[typing.Any]:
+def volume_weight(
+    z: npt.NDArray[np.float64],
+    cosmo: Cosmology,
+) -> npt.NDArray[np.float64]:
     """Uniform weight in comoving volume."""
     return cosmo.xm(z) ** 2 / cosmo.ef(z)  # type: ignore[no-any-return]
 
 
-def density_weight(z: npt.ArrayLike, cosmo: Cosmology) -> npt.NDArray[typing.Any]:
+def density_weight(
+    z: npt.NDArray[np.float64],
+    cosmo: Cosmology,
+) -> npt.NDArray[np.float64]:
     """Uniform weight in matter density."""
     return cosmo.rho_m_z(z) * cosmo.xm(z) ** 2 / cosmo.ef(z)  # type: ignore[no-any-return]
 
@@ -120,9 +128,9 @@ class RadialWindow(typing.NamedTuple):
 
     """
 
-    za: collections.abc.Sequence[float]
-    wa: collections.abc.Sequence[float]
-    zeff: float
+    za: npt.NDArray[np.float64]
+    wa: npt.NDArray[np.float64]
+    zeff: float = 0
 
 
 def tophat_windows(
@@ -176,8 +184,14 @@ def tophat_windows(
         n = max(round((zmax - zmin) / dz), 2)
         z = np.linspace(zmin, zmax, n)
         w = wht(z)
-        zeff = np.trapz(w * z, z) / np.trapz(w, z)  # type: ignore[attr-defined]
-        ws.append(RadialWindow(z, w, zeff))  # type: ignore[arg-type]
+        zeff = np.trapz(  # type: ignore[attr-defined]
+            w * z,
+            z,
+        ) / np.trapz(  # type: ignore[attr-defined]
+            w,
+            z,
+        )
+        ws.append(RadialWindow(z, w, zeff))
     return ws
 
 
@@ -300,7 +314,7 @@ def restrict(
     z: ArrayLike1D,
     f: ArrayLike1D,
     w: RadialWindow,
-) -> tuple[npt.NDArray[typing.Any], npt.NDArray[typing.Any]]:
+) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
     """
     Restrict a function to a redshift window.
 
@@ -332,17 +346,17 @@ 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)  # type: ignore[no-untyped-call]
+    fr = ndinterp(zr, z, f, left=0.0, right=0.0) * ndinterp(zr, w.za, w.wa)
     return zr, fr
 
 
 def partition(
-    z: npt.ArrayLike,
-    fz: npt.ArrayLike,
+    z: npt.NDArray[np.float64],
+    fz: npt.NDArray[np.float64],
     shells: collections.abc.Sequence[RadialWindow],
     *,
     method: str = "nnls",
-) -> npt.ArrayLike:
+) -> npt.NDArray[np.float64]:
     r"""
     Partition a function by a sequence of windows.
 
@@ -448,12 +462,12 @@ def partition(
 
 
 def partition_lstsq(
-    z: npt.ArrayLike,
-    fz: npt.ArrayLike,
+    z: npt.NDArray[np.float64],
+    fz: npt.NDArray[np.float64],
     shells: collections.abc.Sequence[RadialWindow],
     *,
     sumtol: float = 0.01,
-) -> npt.ArrayLike:
+) -> npt.NDArray[np.float64]:
     """Least-squares partition."""
     # make sure nothing breaks
     sumtol = max(sumtol, 1e-4)
@@ -470,36 +484,53 @@ def partition_lstsq(
     dz = np.gradient(zp)
 
     # create the window function matrix
-    a = [np.interp(zp, za, wa, left=0.0, right=0.0) for za, wa, _ in shells]  # type: ignore[arg-type]
-    a /= np.trapz(a, zp, axis=-1)[..., None]  # type: ignore[attr-defined]
+    a = np.array([np.interp(zp, za, wa, left=0.0, right=0.0) for za, wa, _ in shells])
+    a /= np.trapz(  # type: ignore[attr-defined]
+        a,
+        zp,
+        axis=-1,
+    )[..., None]
     a = a * dz
 
     # create the target vector of distribution values
-    b = ndinterp(zp, z, fz, left=0.0, right=0.0)  # type: ignore[no-untyped-call]
+    b = ndinterp(zp, z, fz, left=0.0, right=0.0)
     b = b * dz
 
     # append a constraint for the integral
     mult = 1 / sumtol
-    a = np.concatenate([a, mult * np.ones((len(shells), 1))], axis=-1)  # type: ignore[assignment]
-    b = np.concatenate([b, mult * np.reshape(np.trapz(fz, z), (*dims, 1))], axis=-1)  # type: ignore[attr-defined]
+    a = np.concatenate([a, mult * np.ones((len(shells), 1))], axis=-1)
+    b = np.concatenate(
+        [
+            b,
+            mult
+            * np.reshape(
+                np.trapz(  # type: ignore[attr-defined]
+                    fz,
+                    z,
+                ),
+                (*dims, 1),
+            ),
+        ],
+        axis=-1,
+    )
 
     # now a is a matrix of shape (len(shells), len(zp) + 1)
     # and b is a matrix of shape (*dims, len(zp) + 1)
     # need to find weights x such that b == x @ a over all axes of b
     # do the least-squares fit over partially flattened b, then reshape
-    x = np.linalg.lstsq(a.T, b.reshape(-1, zp.size + 1).T, rcond=None)[0]  # type: ignore[attr-defined, union-attr]
+    x = np.linalg.lstsq(a.T, b.reshape(-1, zp.size + 1).T, rcond=None)[0]
     x = x.T.reshape(*dims, len(shells))
     # roll the last axis of size len(shells) to the front
     return np.moveaxis(x, -1, 0)
 
 
 def partition_nnls(
-    z: npt.ArrayLike,
-    fz: npt.ArrayLike,
+    z: npt.NDArray[np.float64],
+    fz: npt.NDArray[np.float64],
     shells: collections.abc.Sequence[RadialWindow],
     *,
     sumtol: float = 0.01,
-) -> npt.ArrayLike:
+) -> npt.NDArray[np.float64]:
     """
     Non-negative least-squares partition.
 
@@ -524,25 +555,53 @@ def partition_nnls(
     dz = np.gradient(zp)
 
     # create the window function matrix
-    a = [np.interp(zp, za, wa, left=0.0, right=0.0) for za, wa, _ in shells]  # type: ignore[arg-type]
-    a /= np.trapz(a, zp, axis=-1)[..., None]  # type: ignore[attr-defined]
+    a = np.array(
+        [
+            np.interp(
+                zp,
+                za,
+                wa,
+                left=0.0,
+                right=0.0,
+            )
+            for za, wa, _ in shells
+        ],
+    )
+    a /= np.trapz(  # type: ignore[attr-defined]
+        a,
+        zp,
+        axis=-1,
+    )[..., None]
     a = a * dz
 
     # create the target vector of distribution values
-    b = ndinterp(zp, z, fz, left=0.0, right=0.0)  # type: ignore[no-untyped-call]
+    b = ndinterp(zp, z, fz, left=0.0, right=0.0)
     b = b * dz
 
     # append a constraint for the integral
     mult = 1 / sumtol
-    a = np.concatenate([a, mult * np.ones((len(shells), 1))], axis=-1)  # type: ignore[assignment]
-    b = np.concatenate([b, mult * np.reshape(np.trapz(fz, z), (*dims, 1))], axis=-1)  # type: ignore[attr-defined]
+    a = np.concatenate([a, mult * np.ones((len(shells), 1))], axis=-1)
+    b = np.concatenate(
+        [
+            b,
+            mult
+            * np.reshape(
+                np.trapz(  # type: ignore[attr-defined]
+                    fz,
+                    z,
+                ),
+                (*dims, 1),
+            ),
+        ],
+        axis=-1,
+    )
 
     # now a is a matrix of shape (len(shells), len(zp) + 1)
     # and b is a matrix of shape (*dims, len(zp) + 1)
     # for each dim, find non-negative weights x such that b == a.T @ x
 
     # reduce the dimensionality of the problem using a thin QR decomposition
-    q, r = np.linalg.qr(a.T)  # type: ignore[attr-defined]
+    q, r = np.linalg.qr(a.T)
     y = np.einsum("ji,...j", q, b)
 
     # for each dim, find non-negative weights x such that y == r @ x
@@ -555,19 +614,29 @@ def partition_nnls(
 
 
 def partition_restrict(
-    z: npt.ArrayLike,
-    fz: npt.ArrayLike,
+    z: npt.NDArray[np.float64],
+    fz: npt.NDArray[np.float64],
     shells: collections.abc.Sequence[RadialWindow],
-) -> npt.ArrayLike:
+) -> npt.NDArray[np.float64]:
     """Partition by restriction and integration."""
     part = np.empty((len(shells),) + np.shape(fz)[:-1])
     for i, w in enumerate(shells):
-        zr, fr = restrict(z, fz, w)  # type: ignore[arg-type]
-        part[i] = np.trapz(fr, zr, axis=-1)  # type: ignore[attr-defined]
+        zr, fr = restrict(z, fz, w)
+        part[i] = np.trapz(  # type: ignore[attr-defined]
+            fr,
+            zr,
+            axis=-1,
+        )
     return part
 
 
-def redshift_grid(zmin, zmax, *, dz=None, num=None):  # type: ignore[no-untyped-def]
+def redshift_grid(
+    zmin: float,
+    zmax: float,
+    *,
+    dz: float | None = None,
+    num: int | None = None,
+) -> npt.NDArray[np.float64]:
     """Redshift grid with uniform spacing in redshift."""
     if dz is not None and num is None:
         z = np.arange(zmin, np.nextafter(zmax + dz, zmax), dz)
@@ -579,7 +648,14 @@ def redshift_grid(zmin, zmax, *, dz=None, num=None):  # type: ignore[no-untyped-
     return z
 
 
-def distance_grid(cosmo, zmin, zmax, *, dx=None, num=None):  # type: ignore[no-untyped-def]
+def distance_grid(
+    cosmo: Cosmology,
+    zmin: float,
+    zmax: float,
+    *,
+    dx: float | None = None,
+    num: int | None = None,
+) -> npt.NDArray[np.float64]:
     """Redshift grid with uniform spacing in comoving distance."""
     xmin, xmax = cosmo.dc(zmin), cosmo.dc(zmax)
     if dx is not None and num is None:
@@ -589,14 +665,14 @@ def distance_grid(cosmo, zmin, zmax, *, dx=None, num=None):  # type: ignore[no-u
     else:
         msg = 'exactly one of "dx" or "num" must be given'
         raise ValueError(msg)
-    return cosmo.dc_inv(x)
+    return cosmo.dc_inv(x)  # type: ignore[no-any-return]
 
 
 def combine(
-    z: npt.ArrayLike,
-    weights: npt.ArrayLike,
+    z: npt.NDArray[np.float64],
+    weights: npt.NDArray[np.float64],
     shells: collections.abc.Sequence[RadialWindow],
-) -> npt.ArrayLike:
+) -> npt.NDArray[np.float64]:
     r"""
     Evaluate a linear combination of window functions.
 
@@ -626,14 +702,21 @@ def combine(
         Find weights for a given function.
 
     """
-    return sum(
-        np.expand_dims(weight, -1)
-        * np.interp(
-            z,  # type: ignore[arg-type]
-            shell.za,
-            shell.wa / np.trapz(shell.wa, shell.za),  # type: ignore[attr-defined]
-            left=0.0,
-            right=0.0,
-        )
-        for shell, weight in zip(shells, weights)  # type: ignore[arg-type]
+    return np.sum(  # type: ignore[no-any-return]
+        [
+            np.expand_dims(weight, -1)
+            * np.interp(
+                z,
+                shell.za,
+                shell.wa
+                / np.trapz(  # type: ignore[attr-defined]
+                    shell.wa,
+                    shell.za,
+                ),
+                left=0.0,
+                right=0.0,
+            )
+            for shell, weight in zip(shells, weights)
+        ],
+        axis=0,
     )
diff --git a/glass/user.py b/glass/user.py
index 876e70d4..11edd66c 100644
--- a/glass/user.py
+++ b/glass/user.py
@@ -17,12 +17,29 @@
 
 """  # noqa: D205, D400
 
+from __future__ import annotations
+
+import typing
 from contextlib import contextmanager
 
 import numpy as np
+import numpy.typing as npt
+
+if typing.TYPE_CHECKING:
+    import collections.abc
+    import importlib.util
+    import pathlib
+
+    if importlib.util.find_spec("fitsio") is not None:
+        import fitsio
 
 
-def save_cls(filename, cls) -> None:  # type: ignore[no-untyped-def]
+def save_cls(
+    filename: str,
+    cls: collections.abc.Sequence[
+        npt.NDArray[np.float64] | collections.abc.Sequence[float]
+    ],
+) -> None:
     """
     Save a list of Cls to file.
 
@@ -30,12 +47,14 @@ def save_cls(filename, cls) -> None:  # type: ignore[no-untyped-def]
     ``.npz`` suffix, or it will be given one.
 
     """
-    split = np.cumsum([len(cl) if cl is not None else 0 for cl in cls[:-1]])
-    values = np.concatenate([cl for cl in cls if cl is not None])
+    split = np.cumsum([len(cl) for cl in cls[:-1]])
+    values = np.concatenate(cls)
     np.savez(filename, values=values, split=split)
 
 
-def load_cls(filename):  # type: ignore[no-untyped-def]
+def load_cls(
+    filename: str,
+) -> list[npt.NDArray[np.float64] | collections.abc.Sequence[float]]:
     """
     Load a list of Cls from file.
 
@@ -55,12 +74,16 @@ class _FitsWriter:
     Initialised with the fits object and extension name.
     """
 
-    def __init__(self, fits, ext=None) -> None:  # type: ignore[no-untyped-def]
+    def __init__(self, fits: fitsio.FITS, ext: str | None = None) -> None:
         """Create a new, uninitialised writer."""
         self.fits = fits
         self.ext = ext
 
-    def _append(self, data, names=None) -> None:  # type: ignore[no-untyped-def]
+    def _append(
+        self,
+        data: npt.NDArray[np.float64] | list[npt.NDArray[np.float64]],
+        names: list[str] | None = None,
+    ) -> None:
         """Write the FITS file."""
         if self.ext is None or self.ext not in self.fits:
             self.fits.write_table(data, names=names, extname=self.ext)
@@ -71,7 +94,12 @@ def _append(self, data, names=None) -> None:  # type: ignore[no-untyped-def]
             # not using hdu.append here because of incompatibilities
             hdu.write(data, names=names, firstrow=hdu.get_nrows())
 
-    def write(self, data=None, /, **columns) -> None:  # type: ignore[no-untyped-def]
+    def write(
+        self,
+        data: npt.NDArray[np.float64] | None = None,
+        /,
+        **columns: npt.NDArray[np.float64],
+    ) -> None:
         """
         Write to FITS by calling the internal _append method.
 
@@ -89,7 +117,11 @@ def write(self, data=None, /, **columns) -> None:  # type: ignore[no-untyped-def
 
 
 @contextmanager
-def write_catalog(filename, *, ext=None):  # type: ignore[no-untyped-def]
+def write_catalog(
+    filename: pathlib.Path,
+    *,
+    ext: str | None = None,
+) -> collections.abc.Generator[_FitsWriter]:
     """
     Write a catalogue into a FITS file.
 
diff --git a/pyproject.toml b/pyproject.toml
index 40a3d6c9..87788927 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -17,6 +17,7 @@ classifiers = [
     "Programming Language :: Python :: 3.12",
     "Topic :: Scientific/Engineering :: Astronomy",
     "Topic :: Scientific/Engineering :: Physics",
+    "Typing :: Typed",
 ]
 dependencies = [
     "cosmology>=2022.10.9",
@@ -131,10 +132,6 @@ src = [
     "glass",
 ]
 lint.ignore = [
-    "ANN001", # TODO: missing-type-function-argument
-    "ANN002", # TODO: missing-type-args
-    "ANN003", # TODO: missing-type-kwargs
-    "ANN201", # TODO: missing-return-type-undocumented-public-function
     "COM812", # missing-trailing-comma (ruff-format recommended)
     "D203", # one-blank-line-before-class
     "D212", # blank-line-before-class
@@ -168,9 +165,6 @@ lint.per-file-ignores = {"__init__.py" = [
 ], "noxfile.py" = [
     "T201", # print
 ], "tests*" = [
-    "ANN001", # TODO: missing-type-function-argument
-    "ANN201", # TODO: issing-return-type-undocumented-public-function
-    "ANN202", # TODO: missing-return-type-private-function
     "D100", # undocumented-public-module
     "D103", # TODO: undocumented-public-function
     "D104", # undocumented-public-package
diff --git a/tests/conftest.py b/tests/conftest.py
index 511a14b1..b2b6a9fd 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,8 +1,9 @@
+import numpy as np
 import pytest
 
 
 @pytest.fixture(scope="session")
-def rng():  # type: ignore[no-untyped-def]
+def rng() -> np.random.Generator:
     import numpy as np
 
     return np.random.default_rng(seed=42)
diff --git a/tests/core/test_algorithm.py b/tests/core/test_algorithm.py
index 1438c95c..fa1badf9 100644
--- a/tests/core/test_algorithm.py
+++ b/tests/core/test_algorithm.py
@@ -10,7 +10,7 @@
 
 
 @pytest.mark.skipif(not HAVE_SCIPY, reason="test requires SciPy")
-def test_nnls(rng):  # type: ignore[no-untyped-def]
+def test_nnls(rng: np.random.Generator) -> None:
     from scipy.optimize import nnls as nnls_scipy
 
     # cross-check output with scipy's nnls
@@ -21,7 +21,7 @@ def test_nnls(rng):  # type: ignore[no-untyped-def]
     x_glass = nnls_glass(a, b)
     x_scipy, _ = nnls_scipy(a, b)
 
-    np.testing.assert_allclose(x_glass, x_scipy)  # type: ignore[arg-type]
+    np.testing.assert_allclose(x_glass, x_scipy)
 
     # check matrix and vector's shape
 
diff --git a/tests/core/test_array.py b/tests/core/test_array.py
index 315dd204..8258b75e 100644
--- a/tests/core/test_array.py
+++ b/tests/core/test_array.py
@@ -1,6 +1,7 @@
 import importlib.util
 
 import numpy as np
+import numpy.typing as npt
 import pytest
 
 from glass.core.array import (
@@ -15,13 +16,13 @@
 HAVE_SCIPY = importlib.util.find_spec("scipy") is not None
 
 
-def test_broadcast_first():  # type: ignore[no-untyped-def]
+def test_broadcast_first() -> None:
     a = np.ones((2, 3, 4))
     b = np.ones((2, 1))
 
     # 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)  # type: ignore[no-untyped-call]
+    a_a, b_a = broadcast_first(a, b)
     assert a_a.shape == (2, 3, 4)
     assert b_a.shape == (2, 3, 4)
 
@@ -35,7 +36,7 @@ def test_broadcast_first():  # type: ignore[no-untyped-def]
     b = np.ones((5, 6))
 
     with pytest.raises(ValueError, match="shape mismatch"):
-        broadcast_first(a, b)  # type: ignore[no-untyped-call]
+        broadcast_first(a, b)
 
     # plain np.broadcast_arrays will work
     a_a, b_a = np.broadcast_arrays(a, b)
@@ -44,56 +45,57 @@ def test_broadcast_first():  # type: ignore[no-untyped-def]
     assert b_a.shape == (4, 5, 6)
 
 
-def test_broadcast_leading_axes():  # type: ignore[no-untyped-def]
-    a = 0
-    b = np.zeros((4, 10))
-    c = np.zeros((3, 1, 5, 6))
+def test_broadcast_leading_axes() -> None:
+    a_in = 0
+    b_in = np.zeros((4, 10))
+    c_in = np.zeros((3, 1, 5, 6))
 
-    dims, a, b, c = broadcast_leading_axes((a, 0), (b, 1), (c, 2))  # type: ignore[no-untyped-call]
+    dims, *rest = broadcast_leading_axes((a_in, 0), (b_in, 1), (c_in, 2))
+    a_out, b_out, c_out = rest
 
     assert dims == (3, 4)
-    assert a.shape == (3, 4)  # type: ignore[attr-defined]
-    assert b.shape == (3, 4, 10)
-    assert c.shape == (3, 4, 5, 6)
+    assert a_out.shape == (3, 4)
+    assert b_out.shape == (3, 4, 10)
+    assert c_out.shape == (3, 4, 5, 6)
 
 
-def test_ndinterp():  # type: ignore[no-untyped-def]
+def test_ndinterp() -> None:
     # test 1d interpolation
 
-    xp = [0, 1, 2, 3, 4]
-    yp = [1.1, 1.2, 1.3, 1.4, 1.5]
+    xp = np.array([0, 1, 2, 3, 4])
+    yp = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
 
-    x = 0.5
-    y = ndinterp(x, xp, yp)  # type: ignore[no-untyped-call]
+    x: float | npt.NDArray[np.float64] = 0.5
+    y = ndinterp(x, xp, yp)
     assert np.shape(y) == ()
     np.testing.assert_allclose(y, 1.15, atol=1e-15)
 
-    x = [0.5, 1.5, 2.5]  # type: ignore[assignment]
-    y = ndinterp(x, xp, yp)  # type: ignore[no-untyped-call]
+    x = np.array([0.5, 1.5, 2.5])
+    y = ndinterp(x, xp, yp)
     assert np.shape(y) == (3,)
     np.testing.assert_allclose(y, [1.15, 1.25, 1.35], atol=1e-15)
 
-    x = [[0.5, 1.5], [2.5, 3.5]]  # type: ignore[assignment]
-    y = ndinterp(x, xp, yp)  # type: ignore[no-untyped-call]
+    x = np.array([[0.5, 1.5], [2.5, 3.5]])
+    y = 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)
 
     # test nd interpolation in final axis
 
-    yp = [[1.1, 1.2, 1.3, 1.4, 1.5], [2.1, 2.2, 2.3, 2.4, 2.5]]  # type: ignore[list-item]
+    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)  # type: ignore[no-untyped-call]
+    y = ndinterp(x, xp, yp)
     assert np.shape(y) == (2,)
     np.testing.assert_allclose(y, [1.15, 2.15], atol=1e-15)
 
-    x = [0.5, 1.5, 2.5]  # type: ignore[assignment]
-    y = ndinterp(x, xp, yp)  # type: ignore[no-untyped-call]
+    x = np.array([0.5, 1.5, 2.5])
+    y = 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 = [[0.5, 1.5], [2.5, 3.5]]  # type: ignore[assignment]
-    y = ndinterp(x, xp, yp)  # type: ignore[no-untyped-call]
+    x = np.array([[0.5, 1.5], [2.5, 3.5]])
+    y = ndinterp(x, xp, yp)
     assert np.shape(y) == (2, 2, 2)
     np.testing.assert_allclose(
         y,
@@ -103,15 +105,17 @@ def test_ndinterp():  # type: ignore[no-untyped-def]
 
     # test nd interpolation in middle axis
 
-    yp = [[[1.1], [1.2], [1.3], [1.4], [1.5]], [[2.1], [2.2], [2.3], [2.4], [2.5]]]  # type: ignore[list-item]
+    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, axis=1)  # type: ignore[no-untyped-call]
+    y = 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 = [0.5, 1.5, 2.5]  # type: ignore[assignment]
-    y = ndinterp(x, xp, yp, axis=1)  # type: ignore[no-untyped-call]
+    x = np.array([0.5, 1.5, 2.5])
+    y = ndinterp(x, xp, yp, axis=1)
     assert np.shape(y) == (2, 3, 1)
     np.testing.assert_allclose(
         y,
@@ -119,8 +123,8 @@ def test_ndinterp():  # type: ignore[no-untyped-def]
         atol=1e-15,
     )
 
-    x = [[0.5, 1.5, 2.5, 3.5], [3.5, 2.5, 1.5, 0.5], [0.5, 3.5, 1.5, 2.5]]  # type: ignore[assignment]
-    y = ndinterp(x, xp, yp, axis=1)  # type: ignore[no-untyped-call]
+    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)
     assert np.shape(y) == (2, 3, 4, 1)
     np.testing.assert_allclose(
         y,
@@ -140,20 +144,20 @@ def test_ndinterp():  # type: ignore[no-untyped-def]
     )
 
 
-def test_trapz_product():  # type: ignore[no-untyped-def]
+def test_trapz_product() -> None:
     x1 = np.linspace(0, 2, 100)
     f1 = np.full_like(x1, 2.0)
 
     x2 = np.linspace(1, 2, 10)
     f2 = np.full_like(x2, 0.5)
 
-    s = trapz_product((x1, f1), (x2, f2))  # type: ignore[no-untyped-call]
+    s = trapz_product((x1, f1), (x2, f2))
 
     np.testing.assert_allclose(s, 1.0)
 
 
 @pytest.mark.skipif(not HAVE_SCIPY, reason="test requires SciPy")
-def test_cumtrapz():  # type: ignore[no-untyped-def]
+def test_cumtrapz() -> None:
     from scipy.integrate import cumulative_trapezoid
 
     # 1D f and x
@@ -163,18 +167,18 @@ def test_cumtrapz():  # type: ignore[no-untyped-def]
 
     # default dtype (int - not supported by scipy)
 
-    glass_ct = cumtrapz(f, x)  # type: ignore[no-untyped-call]
+    glass_ct = cumtrapz(f, x)
     np.testing.assert_allclose(glass_ct, np.array([0, 1, 4, 7]))
 
     # explicit dtype (float)
 
-    glass_ct = cumtrapz(f, x, dtype=float)  # type: ignore[no-untyped-call]
+    glass_ct = cumtrapz(f, x, dtype=float)
     scipy_ct = cumulative_trapezoid(f, x, initial=0)
     np.testing.assert_allclose(glass_ct, scipy_ct)
 
     # explicit return array
 
-    result = cumtrapz(f, x, dtype=float, out=np.zeros((4,)))  # type: ignore[no-untyped-call]
+    result = cumtrapz(f, x, dtype=float, out=np.zeros((4,)))
     scipy_ct = cumulative_trapezoid(f, x, initial=0)
     np.testing.assert_allclose(result, scipy_ct)
 
@@ -185,17 +189,17 @@ def test_cumtrapz():  # type: ignore[no-untyped-def]
 
     # default dtype (int - not supported by scipy)
 
-    glass_ct = cumtrapz(f, x)  # type: ignore[no-untyped-call]
+    glass_ct = cumtrapz(f, x)
     np.testing.assert_allclose(glass_ct, np.array([[0, 2, 12, 31], [0, 2, 8, 17]]))
 
     # explicit dtype (float)
 
-    glass_ct = cumtrapz(f, x, dtype=float)  # type: ignore[no-untyped-call]
+    glass_ct = cumtrapz(f, x, dtype=float)
     scipy_ct = cumulative_trapezoid(f, x, initial=0)
     np.testing.assert_allclose(glass_ct, scipy_ct)
 
     # explicit return array
 
-    glass_ct = cumtrapz(f, x, dtype=float, out=np.zeros((2, 4)))  # type: ignore[no-untyped-call]
+    glass_ct = cumtrapz(f, x, dtype=float, out=np.zeros((2, 4)))
     scipy_ct = cumulative_trapezoid(f, x, initial=0)
     np.testing.assert_allclose(glass_ct, scipy_ct)
diff --git a/tests/test_fields.py b/tests/test_fields.py
index 18091f1e..46b32843 100644
--- a/tests/test_fields.py
+++ b/tests/test_fields.py
@@ -1,10 +1,16 @@
+import numpy as np
+
 from glass.fields import getcl
 
 
-def test_getcl():  # type: ignore[no-untyped-def]
+def test_getcl() -> None:
     # make a mock Cls array with the index pairs as entries
-    cls = [{i, j} for i in range(10) for j in range(i, -1, -1)]
+    cls = [
+        np.array([i, j], dtype=np.float64) for i in range(10) for j in range(i, -1, -1)
+    ]
     # make sure indices are retrieved correctly
     for i in range(10):
         for j in range(10):
-            assert getcl(cls, i, j) == {i, j}  # type: ignore[no-untyped-call]
+            result = 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)
diff --git a/tests/test_fits.py b/tests/test_fits.py
index 3a6d6229..396bcbae 100644
--- a/tests/test_fits.py
+++ b/tests/test_fits.py
@@ -1,6 +1,8 @@
 import importlib.util
+import pathlib
 
 import numpy as np
+import numpy.typing as npt
 import pytest
 
 from glass import user
@@ -8,17 +10,6 @@
 # check if fitsio is available for testing
 HAVE_FITSIO = importlib.util.find_spec("fitsio") is not None
 
-
-def _test_append(fits, data, names):  # type: ignore[no-untyped-def]
-    """Write routine for FITS test cases."""
-    cat_name = "CATALOG"
-    if cat_name not in fits:
-        fits.write_table(data, names=names, extname=cat_name)
-    else:
-        hdu = fits[cat_name]
-        hdu.write(data, names=names, firstrow=hdu.get_nrows())
-
-
 delta = 0.001  # Number of points in arrays
 my_max = 1000  # Typically number of galaxies in loop
 except_int = 750  # Where test exception occurs in loop
@@ -26,12 +17,25 @@ def _test_append(fits, data, names):  # type: ignore[no-untyped-def]
 
 
 @pytest.mark.skipif(not HAVE_FITSIO, reason="test requires fitsio")
-def test_basic_write(tmp_path):  # type: ignore[no-untyped-def]
+def test_basic_write(tmp_path: pathlib.Path) -> None:
     import fitsio
 
     filename_gfits = "gfits.fits"  # what GLASS creates
     filename_tfits = "tfits.fits"  # file created on the fly to test against
 
+    def _test_append(
+        fits: fitsio.FITS,
+        data: list[npt.NDArray[np.float64]],
+        names: list[str],
+    ) -> None:
+        """Write routine for FITS test cases."""
+        cat_name = "CATALOG"
+        if cat_name not in fits:
+            fits.write_table(data, names=names, extname=cat_name)
+        else:
+            hdu = fits[cat_name]
+            hdu.write(data, names=names, firstrow=hdu.get_nrows())
+
     with (
         user.write_catalog(tmp_path / filename_gfits, ext="CATALOG") as out,
         fitsio.FITS(tmp_path / filename_tfits, "rw", clobber=True) as my_fits,
@@ -42,7 +46,7 @@ def test_basic_write(tmp_path):  # type: ignore[no-untyped-def]
             out.write(RA=array, RB=array2)
             arrays = [array, array2]
             names = ["RA", "RB"]
-            _test_append(my_fits, arrays, names)  # type: ignore[no-untyped-call]
+            _test_append(my_fits, arrays, names)
 
     with (
         fitsio.FITS(tmp_path / filename_gfits) as g_fits,
@@ -55,7 +59,7 @@ def test_basic_write(tmp_path):  # type: ignore[no-untyped-def]
 
 
 @pytest.mark.skipif(not HAVE_FITSIO, reason="test requires fitsio")
-def test_write_exception(tmp_path):  # type: ignore[no-untyped-def]
+def test_write_exception(tmp_path: pathlib.Path) -> None:
     try:
         with user.write_catalog(tmp_path / filename, ext="CATALOG") as out:
             for i in range(my_max):
diff --git a/tests/test_galaxies.py b/tests/test_galaxies.py
index 2f9b2820..a2a347eb 100644
--- a/tests/test_galaxies.py
+++ b/tests/test_galaxies.py
@@ -1,10 +1,12 @@
 import numpy as np
+import numpy.typing as npt
 import pytest
+import pytest_mock
 
 from glass.galaxies import galaxy_shear, gaussian_phz, redshifts, redshifts_from_nz
 
 
-def test_redshifts(mocker):  # type: ignore[no-untyped-def]
+def test_redshifts(mocker: pytest_mock.MockerFixture) -> None:
     # create a mock radial window function
     w = mocker.Mock()
     w.za = np.linspace(0.0, 1.0, 20)
@@ -17,31 +19,51 @@ def test_redshifts(mocker):  # type: ignore[no-untyped-def]
     assert z.max() <= 1.0
 
     # sample redshifts (array)
-    z = redshifts([[1, 2], [3, 4]], w)
+    z = 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(10, [0, 1, 2, 3, 4], [1, 0, 0, 0, 0], warn=False)
+    redshifts = redshifts_from_nz(
+        10,
+        np.array([0, 1, 2, 3, 4]),
+        np.array([1, 0, 0, 0, 0]),
+        warn=False,
+    )
     assert np.all((0 <= redshifts) & (redshifts <= 1))  # noqa: SIM300
 
-    redshifts = redshifts_from_nz(10, [0, 1, 2, 3, 4], [0, 0, 1, 0, 0], warn=False)
+    redshifts = redshifts_from_nz(
+        10,
+        np.array([0, 1, 2, 3, 4]),
+        np.array([0, 0, 1, 0, 0]),
+        warn=False,
+    )
     assert np.all((1 <= redshifts) & (redshifts <= 3))  # noqa: SIM300
 
-    redshifts = redshifts_from_nz(10, [0, 1, 2, 3, 4], [0, 0, 0, 0, 1], warn=False)
+    redshifts = redshifts_from_nz(
+        10,
+        np.array([0, 1, 2, 3, 4]),
+        np.array([0, 0, 0, 0, 1]),
+        warn=False,
+    )
     assert np.all((3 <= redshifts) & (redshifts <= 4))  # noqa: SIM300
 
-    redshifts = redshifts_from_nz(10, [0, 1, 2, 3, 4], [0, 0, 1, 1, 1], warn=False)
+    redshifts = redshifts_from_nz(
+        10,
+        np.array([0, 1, 2, 3, 4]),
+        np.array([0, 0, 1, 1, 1]),
+        warn=False,
+    )
     assert not np.any(redshifts <= 1)
 
     # test with rng
 
     redshifts = redshifts_from_nz(
         10,
-        [0, 1, 2, 3, 4],
-        [0, 0, 1, 1, 1],
+        np.array([0, 1, 2, 3, 4]),
+        np.array([0, 0, 1, 1, 1]),
         warn=False,
         rng=rng,
     )
@@ -51,7 +73,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
 
     # case: no extra dimensions
 
-    count = 10
+    count: int | npt.NDArray[np.float64] = 10
     z = np.linspace(0, 1, 100)
     nz = z * (1 - z)
 
@@ -62,7 +84,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
 
     # case: extra dimensions from count
 
-    count = [10, 20, 30]  # type: ignore[assignment]
+    count = np.array([10, 20, 30])
     z = np.linspace(0, 1, 100)
     nz = z * (1 - z)
 
@@ -74,7 +96,7 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
 
     count = 10
     z = np.linspace(0, 1, 100)
-    nz = [z * (1 - z), (z - 0.5) ** 2]  # type: ignore[assignment]
+    nz = np.array([z * (1 - z), (z - 0.5) ** 2])
 
     redshifts = redshifts_from_nz(count, z, nz, warn=False)
 
@@ -82,9 +104,9 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
 
     # case: extra dimensions from count and nz
 
-    count = [[10], [20], [30]]  # type: ignore[assignment]
+    count = np.array([[10], [20], [30]])
     z = np.linspace(0, 1, 100)
-    nz = [z * (1 - z), (z - 0.5) ** 2]  # type: ignore[assignment]
+    nz = np.array([z * (1 - z), (z - 0.5) ** 2])
 
     redshifts = redshifts_from_nz(count, z, nz, warn=False)
 
@@ -92,18 +114,22 @@ def test_redshifts_from_nz(rng: np.random.Generator) -> None:
 
     # case: incompatible input shapes
 
-    count = [10, 20, 30]  # type: ignore[assignment]
+    count = np.array([10, 20, 30])
     z = np.linspace(0, 1, 100)
-    nz = [z * (1 - z), (z - 0.5) ** 2]  # type: ignore[assignment]
+    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)
 
     with pytest.warns(UserWarning, match="when sampling galaxies"):
-        redshifts = redshifts_from_nz(10, [0, 1, 2, 3, 4], [1, 0, 0, 0, 0])
+        redshifts = redshifts_from_nz(
+            10,
+            np.array([0, 1, 2, 3, 4]),
+            np.array([1, 0, 0, 0, 0]),
+        )
 
 
-def test_galaxy_shear(rng):  # type: ignore[no-untyped-def]
+def test_galaxy_shear(rng: np.random.Generator) -> None:
     # check shape of the output
 
     kappa, gamma1, gamma2 = (
@@ -112,7 +138,14 @@ def test_galaxy_shear(rng):  # type: ignore[no-untyped-def]
         rng.normal(size=(12,)),
     )
 
-    shear = galaxy_shear([], [], [], kappa, gamma1, gamma2)  # type: ignore[arg-type]
+    shear = galaxy_shear(
+        np.array([]),
+        np.array([]),
+        np.array([]),
+        kappa,
+        gamma1,
+        gamma2,
+    )
     np.testing.assert_equal(shear, [])
 
     gal_lon, gal_lat, gal_eps = (
@@ -125,7 +158,15 @@ def test_galaxy_shear(rng):  # type: ignore[no-untyped-def]
 
     # shape with no reduced shear
 
-    shear = galaxy_shear([], [], [], kappa, gamma1, gamma2, reduced_shear=False)  # type: ignore[arg-type]
+    shear = galaxy_shear(
+        np.array([]),
+        np.array([]),
+        np.array([]),
+        kappa,
+        gamma1,
+        gamma2,
+        reduced_shear=False,
+    )
     np.testing.assert_equal(shear, [])
 
     gal_lon, gal_lat, gal_eps = (
@@ -150,8 +191,8 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
 
     # case: zero variance
 
-    z = np.linspace(0, 1, 100)
-    sigma_0 = 0.0
+    z: float | npt.NDArray[np.float64] = np.linspace(0, 1, 100)
+    sigma_0: float | npt.NDArray[np.float64] = 0.0
 
     phz = gaussian_phz(z, sigma_0)
 
@@ -165,21 +206,23 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
 
     # case: truncated normal
 
-    z = 0.0  # type: ignore[assignment]
-    sigma_0 = np.ones(100)  # type: ignore[assignment]
+    z = 0.0
+    sigma_0 = np.ones(100)
 
     phz = gaussian_phz(z, sigma_0)
 
+    assert isinstance(phz, np.ndarray)
     assert phz.shape == (100,)
     assert np.all(phz >= 0)
 
     # case: upper and lower bound
 
-    z = 1.0  # type: ignore[assignment]
-    sigma_0 = np.ones(100)  # type: ignore[assignment]
+    z = 1.0
+    sigma_0 = np.ones(100)
 
     phz = gaussian_phz(z, sigma_0, lower=0.5, upper=1.5)
 
+    assert isinstance(phz, np.ndarray)
     assert phz.shape == (100,)
     assert np.all(phz >= 0.5)
     assert np.all(phz <= 1.5)
@@ -188,7 +231,7 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
 
     # case: scalar redshift, scalar sigma_0
 
-    z = 1.0  # type: ignore[assignment]
+    z = 1.0
     sigma_0 = 0.0
 
     phz = gaussian_phz(z, sigma_0)
@@ -203,35 +246,38 @@ def test_gaussian_phz(rng: np.random.Generator) -> None:
 
     phz = gaussian_phz(z, sigma_0)
 
+    assert isinstance(phz, np.ndarray)
     assert phz.shape == (10,)
     np.testing.assert_array_equal(z, phz)
 
     # case: scalar redshift, array sigma_0
 
-    z = 1.0  # type: ignore[assignment]
-    sigma_0 = np.zeros(10)  # type: ignore[assignment]
+    z = 1.0
+    sigma_0 = np.zeros(10)
 
     phz = gaussian_phz(z, sigma_0)
 
+    assert isinstance(phz, np.ndarray)
     assert phz.shape == (10,)
     np.testing.assert_array_equal(z, phz)
 
     # case: array redshift, array sigma_0
 
     z = np.linspace(0, 1, 10)
-    sigma_0 = np.zeros((11, 1))  # type: ignore[assignment]
+    sigma_0 = np.zeros((11, 1))
 
     phz = gaussian_phz(z, sigma_0)
 
+    assert isinstance(phz, np.ndarray)
     assert phz.shape == (11, 10)
     np.testing.assert_array_equal(np.broadcast_to(z, (11, 10)), phz)
 
     # test resampling
 
     phz = gaussian_phz(np.array(0.0), np.array(1.0), lower=np.array([0]))
-    assert isinstance(phz, float)  # type: ignore[unreachable]
+    assert isinstance(phz, float)
 
-    phz = gaussian_phz(np.array(0.0), np.array(1.0), upper=np.array([1]))  # type: ignore[unreachable]
+    phz = gaussian_phz(np.array(0.0), np.array(1.0), upper=np.array([1]))
     assert isinstance(phz, float)
 
     # test error
diff --git a/tests/test_lensing.py b/tests/test_lensing.py
index 96bc1944..356b2b78 100644
--- a/tests/test_lensing.py
+++ b/tests/test_lensing.py
@@ -1,5 +1,10 @@
+from __future__ import annotations
+
+import typing
+
 import healpix
 import numpy as np
+import numpy.typing as npt
 import pytest
 
 from glass.lensing import (
@@ -10,29 +15,36 @@
 )
 from glass.shells import RadialWindow
 
+if typing.TYPE_CHECKING:
+    from cosmology import Cosmology
+
 
 @pytest.fixture
-def shells():  # type: ignore[no-untyped-def]
+def shells() -> list[RadialWindow]:
     return [
-        RadialWindow([0.0, 1.0, 2.0], [0.0, 1.0, 0.0], 1.0),
-        RadialWindow([1.0, 2.0, 3.0], [0.0, 1.0, 0.0], 2.0),
-        RadialWindow([2.0, 3.0, 4.0], [0.0, 1.0, 0.0], 3.0),
-        RadialWindow([3.0, 4.0, 5.0], [0.0, 1.0, 0.0], 4.0),
-        RadialWindow([4.0, 5.0, 6.0], [0.0, 1.0, 0.0], 5.0),
+        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),
     ]
 
 
 @pytest.fixture
-def cosmo():  # type: ignore[no-untyped-def]
+def cosmo() -> Cosmology:
     class MockCosmology:
         @property
-        def omega_m(self):  # type: ignore[no-untyped-def]
+        def omega_m(self) -> float:
             return 0.3
 
-        def ef(self, z):  # type: ignore[no-untyped-def]
+        def ef(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
             return (self.omega_m * (1 + z) ** 3 + 1 - self.omega_m) ** 0.5
 
-        def xm(self, z, z2=None):  # type: ignore[no-untyped-def]
+        def xm(
+            self,
+            z: npt.NDArray[np.float64],
+            z2: npt.NDArray[np.float64] | None = None,
+        ) -> npt.NDArray[np.float64]:
             if z2 is None:
                 return np.array(z) * 1000
             return (np.array(z2) - np.array(z)) * 1000
@@ -41,37 +53,31 @@ def xm(self, z, z2=None):  # type: ignore[no-untyped-def]
 
 
 @pytest.mark.parametrize("usecomplex", [True, False])
-def test_deflect_nsew(usecomplex):  # type: ignore[no-untyped-def]
+def test_deflect_nsew(usecomplex: bool) -> None:  # noqa: FBT001
     d = 5.0
     r = np.radians(d)
 
-    if usecomplex:
-
-        def alpha(re, im):  # type: ignore[no-untyped-def]
-            return re + 1j * im
-    else:
-
-        def alpha(re, im):  # type: ignore[no-untyped-def]
-            return [re, im]
+    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))  # type: ignore[no-untyped-call]
+    lon, lat = 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))  # type: ignore[no-untyped-call]
+    lon, lat = 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))  # type: ignore[no-untyped-call]
+    lon, lat = 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))  # type: ignore[no-untyped-call]
+    lon, lat = deflect(0.0, 0.0, alpha(0, -r, usecomplex=usecomplex))
     np.testing.assert_allclose([lon, lat], [d, 0.0], atol=1e-15)
 
 
-def test_deflect_many(rng):  # type: ignore[no-untyped-def]
+def test_deflect_many(rng: np.random.Generator) -> None:
     n = 1000
     abs_alpha = rng.uniform(0, 2 * np.pi, size=n)
     arg_alpha = rng.uniform(-np.pi, np.pi, size=n)
@@ -89,7 +95,11 @@ def test_deflect_many(rng):  # type: ignore[no-untyped-def]
     np.testing.assert_allclose(dotp, np.cos(abs_alpha))
 
 
-def test_multi_plane_matrix(shells, cosmo, rng):  # type: ignore[no-untyped-def]
+def test_multi_plane_matrix(
+    shells: list[RadialWindow],
+    cosmo: Cosmology,
+    rng: np.random.Generator,
+) -> None:
     mat = multi_plane_matrix(shells, cosmo)
 
     np.testing.assert_array_equal(mat, np.tril(mat))
@@ -101,12 +111,17 @@ def test_multi_plane_matrix(shells, cosmo, rng):  # type: ignore[no-untyped-def]
     kappas = []
     for shell, delta in zip(shells, deltas):
         convergence.add_window(delta, shell)
-        kappas.append(convergence.kappa.copy())  # type: ignore[union-attr]
+        if convergence.kappa is not None:
+            kappas.append(convergence.kappa.copy())
 
     np.testing.assert_allclose(mat @ deltas, kappas)
 
 
-def test_multi_plane_weights(shells, cosmo, rng):  # type: ignore[no-untyped-def]
+def test_multi_plane_weights(
+    shells: list[RadialWindow],
+    cosmo: Cosmology,
+    rng: np.random.Generator,
+) -> None:
     w_in = np.eye(len(shells))
     w_out = multi_plane_weights(w_in, shells, cosmo)
 
@@ -125,4 +140,4 @@ def test_multi_plane_weights(shells, cosmo, rng):  # type: ignore[no-untyped-def
 
     wmat = multi_plane_weights(weights, shells, cosmo)
 
-    np.testing.assert_allclose(np.einsum("ij,ik", wmat, deltas), kappa)  # type: ignore[arg-type]
+    np.testing.assert_allclose(np.einsum("ij,ik", wmat, deltas), kappa)
diff --git a/tests/test_points.py b/tests/test_points.py
index 584e824d..d76a8204 100644
--- a/tests/test_points.py
+++ b/tests/test_points.py
@@ -1,39 +1,62 @@
+from __future__ import annotations
+
+import typing
+
 import numpy as np
+import numpy.typing as npt
 
 from glass.points import position_weights, positions_from_delta, uniform_positions
 
-
-def catpos(pos):  # type: ignore[no-untyped-def]
-    lon, lat, cnt = [], [], 0  # type: ignore[var-annotated]
+if typing.TYPE_CHECKING:
+    import collections.abc
+
+
+def catpos(
+    pos: collections.abc.Generator[
+        tuple[
+            npt.NDArray[np.float64],
+            npt.NDArray[np.float64],
+            int | npt.NDArray[np.int_],
+        ]
+    ],
+) -> tuple[
+    npt.NDArray[np.float64],
+    npt.NDArray[np.float64],
+    int | npt.NDArray[np.int_],
+]:
+    lon = np.empty(0)
+    lat = np.empty(0)
+    cnt: int | npt.NDArray[np.int_] = 0
     for lo, la, co in pos:
-        lon = np.concatenate([lon, lo])  # type: ignore[assignment]
-        lat = np.concatenate([lat, la])  # type: ignore[assignment]
+        lon = np.concatenate([lon, lo])
+        lat = np.concatenate([lat, la])
         cnt = cnt + co
     return lon, lat, cnt
 
 
-def test_positions_from_delta():  # type: ignore[no-untyped-def]
+def test_positions_from_delta() -> None:
     # case: single-dimensional input
 
-    ngal = 1e-3
+    ngal: float | npt.NDArray[np.float64] = 1e-3
     delta = np.zeros(12)
     bias = 0.8
     vis = np.ones(12)
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))  # type: ignore[no-untyped-call]
+    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))
 
     assert isinstance(cnt, int)
     assert lon.shape == lat.shape == (cnt,)
 
     # case: multi-dimensional ngal
 
-    ngal = [1e-3, 2e-3]  # type: ignore[assignment]
+    ngal = np.array([1e-3, 2e-3])
     delta = np.zeros(12)
     bias = 0.8
     vis = np.ones(12)
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))  # type: ignore[no-untyped-call]
+    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))
 
+    assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (2,)
     assert lon.shape == (cnt.sum(),)
     assert lat.shape == (cnt.sum(),)
@@ -45,62 +68,66 @@ def test_positions_from_delta():  # type: ignore[no-untyped-def]
     bias = 0.8
     vis = np.ones(12)
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))  # type: ignore[no-untyped-call]
+    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))
 
+    assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (3, 2)
     assert lon.shape == (cnt.sum(),)
     assert lat.shape == (cnt.sum(),)
 
     # case: multi-dimensional broadcasting
 
-    ngal = [1e-3, 2e-3]  # type: ignore[assignment]
+    ngal = np.array([1e-3, 2e-3])
     delta = np.zeros((3, 1, 12))
     bias = 0.8
     vis = np.ones(12)
 
-    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))  # type: ignore[no-untyped-call]
+    lon, lat, cnt = catpos(positions_from_delta(ngal, delta, bias, vis))
 
+    assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (3, 2)
     assert lon.shape == (cnt.sum(),)
     assert lat.shape == (cnt.sum(),)
 
 
-def test_uniform_positions():  # type: ignore[no-untyped-def]
+def test_uniform_positions() -> None:
     # case: scalar input
 
-    ngal = 1e-3
+    ngal: float | npt.NDArray[np.float64] = 1e-3
 
-    lon, lat, cnt = catpos(uniform_positions(ngal))  # type: ignore[no-untyped-call]
+    lon, lat, cnt = catpos(uniform_positions(ngal))
 
     assert isinstance(cnt, int)
     assert lon.shape == lat.shape == (cnt,)
 
     # case: 1-D array input
 
-    ngal = [1e-3, 2e-3, 3e-3]  # type: ignore[assignment]
+    ngal = np.array([1e-3, 2e-3, 3e-3])
 
-    lon, lat, cnt = catpos(uniform_positions(ngal))  # type: ignore[no-untyped-call]
+    lon, lat, cnt = catpos(uniform_positions(ngal))
 
+    assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (3,)
     assert lon.shape == lat.shape == (cnt.sum(),)
 
     # case: 2-D array input
 
-    ngal = [[1e-3, 2e-3], [3e-3, 4e-3], [5e-3, 6e-3]]  # type: ignore[assignment]
+    ngal = np.array([[1e-3, 2e-3], [3e-3, 4e-3], [5e-3, 6e-3]])
 
-    lon, lat, cnt = catpos(uniform_positions(ngal))  # type: ignore[no-untyped-call]
+    lon, lat, cnt = catpos(uniform_positions(ngal))
 
+    assert isinstance(cnt, np.ndarray)
     assert cnt.shape == (3, 2)
     assert lon.shape == lat.shape == (cnt.sum(),)
 
 
-def test_position_weights(rng):  # type: ignore[no-untyped-def]
+def test_position_weights(rng: np.random.Generator) -> None:
     for bshape in None, (), (100,), (100, 1):
         for cshape in (100,), (100, 50), (100, 3, 2):
             counts = rng.random(cshape)
             bias = None if bshape is None else rng.random(bshape)
 
-            weights = position_weights(counts, bias)  # type: ignore[no-untyped-call]
+            weights = 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 5f3bf3bc..e965fdce 100644
--- a/tests/test_shapes.py
+++ b/tests/test_shapes.py
@@ -9,31 +9,31 @@
 )
 
 
-def test_triaxial_axis_ratio(rng):  # type: ignore[no-untyped-def]
+def test_triaxial_axis_ratio(rng: np.random.Generator) -> None:
     # single axis ratio
 
-    q = triaxial_axis_ratio(0.8, 0.4)  # type: ignore[no-untyped-call]
+    q = triaxial_axis_ratio(0.8, 0.4)
     assert np.isscalar(q)
 
     # many axis ratios
 
-    q = triaxial_axis_ratio(0.8, 0.4, size=1000)  # type: ignore[no-untyped-call]
+    q = 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))  # type: ignore[no-untyped-call]
+    q = triaxial_axis_ratio(0.8, 0.4, size=(10, 10))
     assert np.shape(q) == (10, 10)
 
     # implicit size
 
-    q1 = triaxial_axis_ratio([0.8, 0.9], 0.4)  # type: ignore[no-untyped-call]
-    q2 = triaxial_axis_ratio(0.8, [0.4, 0.5])  # type: ignore[no-untyped-call]
+    q1 = triaxial_axis_ratio(np.array([0.8, 0.9]), 0.4)
+    q2 = 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([[0.6, 0.7], [0.8, 0.9]], [0.4, 0.5])  # type: ignore[no-untyped-call]
+    q = 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,50 +42,55 @@ def test_triaxial_axis_ratio(rng):  # type: ignore[no-untyped-def]
     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)  # type: ignore[no-untyped-call]
+    q = 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)  # type: ignore[no-untyped-call]
+    e = 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)  # type: ignore[no-untyped-call]
+    e = 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)  # type: ignore[no-untyped-call]
+    e = 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))  # type: ignore[no-untyped-call]
+    e = 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, [0.222, 0.333], 0.056)  # type: ignore[no-untyped-call]
-    e2 = ellipticity_ryden04(-1.85, 0.89, 0.222, [0.056, 0.067])  # type: ignore[no-untyped-call]
-    e3 = ellipticity_ryden04([-1.85, -2.85], 0.89, 0.222, 0.056)  # type: ignore[no-untyped-call]
-    e4 = ellipticity_ryden04(-1.85, [0.89, 1.001], 0.222, 0.056)  # type: ignore[no-untyped-call]
+    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)
     assert np.shape(e1) == np.shape(e2) == np.shape(e3) == np.shape(e4) == (2,)
 
     # broadcasting rule
 
-    e = ellipticity_ryden04([-1.9, -2.9], 0.9, [[0.2, 0.3], [0.4, 0.5]], 0.1)  # type: ignore[no-untyped-call]
+    e = ellipticity_ryden04(
+        np.array([-1.9, -2.9]),
+        0.9,
+        np.array([[0.2, 0.3], [0.4, 0.5]]),
+        0.1,
+    )
     assert np.shape(e) == (2, 2)
 
     # check that result is in the specified range
 
-    e = ellipticity_ryden04(0.0, 1.0, 0.222, 0.056, size=10)  # type: ignore[no-untyped-call]
+    e = 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)  # type: ignore[no-untyped-call]
+    e = ellipticity_ryden04(0.0, 1.0, 0.0, 1.0, size=10)
     assert np.all((e.real >= -1.0) & (e.real <= 1.0))
 
 
@@ -108,7 +113,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([n, n], [0.128, 0.256])
+    eps = ellipticity_gaussian(np.array([n, n]), np.array([0.128, 0.256]))
 
     assert eps.shape == (2 * n,)
 
@@ -138,7 +143,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([n, n], [0.128, 0.256])
+    eps = ellipticity_intnorm(np.array([n, n]), np.array([0.128, 0.256]))
 
     assert eps.shape == (2 * n,)
 
diff --git a/tests/test_shells.py b/tests/test_shells.py
index f468d3f6..5f15f5ed 100644
--- a/tests/test_shells.py
+++ b/tests/test_shells.py
@@ -4,8 +4,8 @@
 from glass.shells import RadialWindow, partition, restrict, tophat_windows
 
 
-def test_tophat_windows():  # type: ignore[no-untyped-def]
-    zb = [0.0, 0.1, 0.2, 0.5, 1.0, 2.0]
+def test_tophat_windows() -> None:
+    zb = np.array([0.0, 0.1, 0.2, 0.5, 1.0, 2.0])
     dz = 0.005
 
     ws = tophat_windows(zb, dz)
@@ -18,16 +18,19 @@ def test_tophat_windows():  # type: ignore[no-untyped-def]
         zn <= z0 + len(w.za) * dz <= zn + dz for w, z0, zn in zip(ws, zb, zb[1:])
     )
 
-    assert all(np.all(w.wa == 1) for w in ws)  # type: ignore[comparison-overlap]
+    assert all(np.all(w.wa == 1) for w in ws)
 
 
-def test_restrict():  # type: ignore[no-untyped-def]
+def test_restrict() -> None:
     # 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(za=[1.0, 2.0, 3.0, 4.0], wa=[0.0, 0.5, 0.5, 0.0], zeff=None)  # type: ignore[arg-type]
+    w = 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)
 
@@ -49,34 +52,14 @@ def test_restrict():  # type: ignore[no-untyped-def]
 
 
 @pytest.mark.parametrize("method", ["lstsq", "nnls", "restrict"])
-def test_partition(method):  # type: ignore[no-untyped-def]
+def test_partition(method: str) -> None:
     shells = [
-        RadialWindow(np.array([0.0, 1.0]), np.array([1.0, 0.0]), 0.0),  # type: ignore[arg-type]
-        RadialWindow(
-            np.array([0.0, 1.0, 2.0]),  # type: ignore[arg-type]
-            np.array([0.0, 1.0, 0.0]),  # type: ignore[arg-type]
-            0.5,
-        ),
-        RadialWindow(
-            np.array([1.0, 2.0, 3.0]),  # type: ignore[arg-type]
-            np.array([0.0, 1.0, 0.0]),  # type: ignore[arg-type]
-            1.5,
-        ),
-        RadialWindow(
-            np.array([2.0, 3.0, 4.0]),  # type: ignore[arg-type]
-            np.array([0.0, 1.0, 0.0]),  # type: ignore[arg-type]
-            2.5,
-        ),
-        RadialWindow(
-            np.array([3.0, 4.0, 5.0]),  # type: ignore[arg-type]
-            np.array([0.0, 1.0, 0.0]),  # type: ignore[arg-type]
-            3.5,
-        ),
-        RadialWindow(
-            np.array([4.0, 5.0]),  # type: ignore[arg-type]
-            np.array([0.0, 1.0]),  # type: ignore[arg-type]
-            5.0,
-        ),
+        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),
     ]
 
     z = np.linspace(0.0, 5.0, 1000)
@@ -87,6 +70,12 @@ def test_partition(method):  # type: ignore[no-untyped-def]
 
     part = partition(z, fz, shells, method=method)
 
-    assert part.shape == (len(shells), 3, 2)  # type: ignore[union-attr]
+    assert part.shape == (len(shells), 3, 2)
 
-    np.testing.assert_allclose(part.sum(axis=0), np.trapz(fz, z))  # type: ignore[attr-defined, union-attr]
+    np.testing.assert_allclose(
+        part.sum(axis=0),
+        np.trapz(  # type: ignore[attr-defined]
+            fz,
+            z,
+        ),
+    )