ENH(shells): radial windows for linear and cubic interpolation (#119)

Adds window functions `linear_windows()` and `cubic_windows()` that
correspond to linear and cubic spline interpolation, respectively.

Documentation is updated to give an overview of available window
functions.

Closes: #118 
Added: New `linear_windows()` and `cubic_windows()` window functions for
  shells.

docs/user/how-glass-works.rst

@@ -44,7 +44,6 @@ which are flat and non-overlapping.
         plt.fill_between(za, np.zeros_like(wa), wa, alpha=0.5)
         plt.annotate(f'shell {i+1}', (zeff, 0.5), ha='center', va='center')
 
-    plt.ylim(0., 2.)
     plt.xlabel('redshift $z$')
     plt.ylabel('window function $W(z)$')
     plt.tight_layout()
@@ -63,6 +62,48 @@ then simulates the (radially) continuous field :math:`F(z)` as the (radially)
 discretised fields :math:`F_i`.
 
 
+.. _user-window-functions:
+
+Window functions
+^^^^^^^^^^^^^^^^
+
+*GLASS* supports arbitrary window functions (although the computation of
+:ref:`line-of-sight integrals <user-los-integrals>` makes some assumptions).
+The following :ref:`window functions <reference-window-functions>` are
+included:
+
+.. plot::
+
+    from glass.shells import (redshift_grid, tophat_windows, linear_windows,
+                              cubic_windows)
+
+    plot_windows = [tophat_windows, linear_windows,
+                    cubic_windows]
+    nr = (len(plot_windows)+1)//2
+
+    fig, axes = plt.subplots(nr, 2, figsize=(8, nr*3), layout="constrained",
+                             squeeze=False, sharex=False, sharey=True)
+
+    zs = redshift_grid(0., 0.5, dz=0.1)
+    zt = np.linspace(0., 0.5, 200)
+
+    for ax in axes.flat:
+        ax.axis(False)
+    for windows, ax in zip(plot_windows, axes.flat):
+        ws = windows(zs)
+        wt = np.zeros_like(zt)
+        ax.axis(True)
+        ax.set_title(windows.__name__)
+        for i, (za, wa, zeff) in enumerate(ws):
+            wt += np.interp(zt, za, wa, left=0., right=0.)
+            ax.fill_between(za, np.zeros_like(wa), wa, alpha=0.5)
+        ax.plot(zt, wt, c="k", lw=2)
+    for ax in axes.flat:
+        ax.set_xlabel("redshift $z$")
+    for ax in axes[:, 0]:
+        ax.set_ylabel("window function $W(z)$")
+
+
 Angular discretisation
 ----------------------
 
@@ -89,6 +130,8 @@ difference between the two is whether or not the pixel window function is
 applied to the spherical harmonic expansion of the fields.
 
 
+.. _user-los-integrals:
+
 Line-of-sight integrals
 -----------------------
 

glass/shells.py

@@ -10,11 +10,15 @@
 matter shells, i.e. the radial discretisation of the light cone.
 
 
+.. _reference-window-functions:
+
 Window functions
 ----------------
 
 .. autoclass:: RadialWindow
 .. autofunction:: tophat_windows
+.. autofunction:: linear_windows
+.. autofunction:: cubic_windows
 
 
 Window function tools
@@ -151,6 +155,10 @@ def tophat_windows(zbins: ArrayLike1D, dz: float = 1e-3,
     ws : (N,) list of :class:`RadialWindow`
         List of window functions.
 
+    See Also
+    --------
+    :ref:`user-window-functions`
+
     '''
     if len(zbins) < 2:
         raise ValueError('zbins must have at least two entries')
@@ -173,6 +181,118 @@ def tophat_windows(zbins: ArrayLike1D, dz: float = 1e-3,
     return ws
 
 
+def linear_windows(zgrid: ArrayLike1D, dz: float = 1e-3,
+                   weight: Optional[WeightFunc] = None,
+                   ) -> List[RadialWindow]:
+    '''Linear interpolation window functions.
+
+    Uses the *N+2* given redshifts as nodes to construct *N* triangular
+    window functions between the first and last node.  These correspond
+    to linear interpolation of radial functions.  The redshift spacing
+    of the windows is approximately equal to ``dz``.
+
+    An optional weight function :math:`w(z)` can be given using
+    ``weight``; it is applied to the triangular windows.
+
+    The resulting windows functions are :class:`RadialWindow` instances.
+    Their effective redshifts correspond to the given nodes.
+
+    Parameters
+    ----------
+    zgrid : (N+2,) array_like
+        Redshift grid for the triangular window functions.
+    dz : float, optional
+        Approximate spacing of the redshift grid.
+    weight : callable, optional
+        If given, a weight function to be applied to the window
+        functions.
+
+    Returns
+    -------
+    ws : (N,) list of :class:`RadialWindow`
+        List of window functions.
+
+    See Also
+    --------
+    :ref:`user-window-functions`
+
+    '''
+    if len(zgrid) < 3:
+        raise ValueError('nodes must have at least 3 entries')
+    if zgrid[0] != 0:
+        warnings.warn('first triangular window does not start at z=0')
+
+    ws = []
+    for zmin, zmid, zmax in zip(zgrid, zgrid[1:], zgrid[2:]):
+        n = max(round((zmid - zmin)/dz), 2) - 1
+        m = max(round((zmax - zmid)/dz), 2)
+        z = np.concatenate([np.linspace(zmin, zmid, n, endpoint=False),
+                            np.linspace(zmid, zmax, m)])
+        w = np.concatenate([np.linspace(0., 1., n, endpoint=False),
+                            np.linspace(1., 0., m)])
+        if weight is not None:
+            w *= weight(z)
+        ws.append(RadialWindow(z, w, zmid))
+    return ws
+
+
+def cubic_windows(zgrid: ArrayLike1D, dz: float = 1e-3,
+                  weight: Optional[WeightFunc] = None,
+                  ) -> List[RadialWindow]:
+    '''Cubic interpolation window functions.
+
+    Uses the *N+2* given redshifts as nodes to construct *N* cubic
+    Hermite spline window functions between the first and last node.
+    These correspond to cubic spline interpolation of radial functions.
+    The redshift spacing of the windows is approximately equal to
+    ``dz``.
+
+    An optional weight function :math:`w(z)` can be given using
+    ``weight``; it is applied to the cubic spline windows.
+
+    The resulting windows functions are :class:`RadialWindow` instances.
+    Their effective redshifts correspond to the given nodes.
+
+    Parameters
+    ----------
+    zgrid : (N+2,) array_like
+        Redshift grid for the cubic spline window functions.
+    dz : float, optional
+        Approximate spacing of the redshift grid.
+    weight : callable, optional
+        If given, a weight function to be applied to the window
+        functions.
+
+    Returns
+    -------
+    ws : (N,) list of :class:`RadialWindow`
+        List of window functions.
+
+    See Also
+    --------
+    :ref:`user-window-functions`
+
+    '''
+    if len(zgrid) < 3:
+        raise ValueError('nodes must have at least 3 entries')
+    if zgrid[0] != 0:
+        warnings.warn('first cubic spline window does not start at z=0')
+
+    ws = []
+    for zmin, zmid, zmax in zip(zgrid, zgrid[1:], zgrid[2:]):
+        n = max(round((zmid - zmin)/dz), 2) - 1
+        m = max(round((zmax - zmid)/dz), 2)
+        z = np.concatenate([np.linspace(zmin, zmid, n, endpoint=False),
+                            np.linspace(zmid, zmax, m)])
+        u = np.linspace(0., 1., n, endpoint=False)
+        v = np.linspace(1., 0., m)
+        w = np.concatenate([u**2*(3-2*u), v**2*(3-2*v)])
+        if weight is not None:
+            w *= weight(z)
+        ws.append(RadialWindow(z, w, zmid))
+    return ws
+
+
 def restrict(z: ArrayLike1D, f: ArrayLike1D, w: RadialWindow
              ) -> Tuple[np.ndarray, np.ndarray]:
     '''Restrict a function to a redshift window.