diff --git a/.readthedocs.yml b/.readthedocs.yml
index 8076f03..7e47f55 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yml
@@ -1,5 +1,10 @@
 version: 2
 
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.11"
+
 python:
   install:
     - requirements: requirements.txt
diff --git a/examples/advanced/plot_s4_galaxies.py b/examples/advanced/plot_s4_galaxies.py
index 2227cdd..afc3266 100644
--- a/examples/advanced/plot_s4_galaxies.py
+++ b/examples/advanced/plot_s4_galaxies.py
@@ -2,9 +2,9 @@
 Stage IV Galaxy Survey
 ======================
 
-This example simulates a galaxy catalogue from a Stage IV Space Satellite Galaxy
-Survey such as *Euclid* and *Roman* combining the :doc:`/basic/plot_density` and
-:doc:`/basic/plot_lensing` examples with galaxy ellipticities and galaxy shears,
+This example simulates a galaxy catalogue from a Stage IV Space Satellite
+Galaxy Survey such as *Euclid* and *Roman* combining the :doc:`/basic/galaxies`
+and :doc:`/basic/lensing` examples with galaxy ellipticities and galaxy shears,
 as well as using some auxiliary functions.
 
 The focus in this example is mock catalogue generation using auxiliary functions
@@ -71,9 +71,12 @@
 # compute the angular matter power spectra of the shells with CAMB
 cls = glass.ext.camb.matter_cls(pars, lmax, ws)
 
-# compute Gaussian cls for lognormal fields for 3 correlated shells
+# angular discretisation with 3 correlated shells
 # putting nside here means that the HEALPix pixel window function is applied
-gls = glass.fields.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)
+cls = glass.fields.discretized_cls(cls, nside=nside, lmax=lmax, ncorr=3)
+
+# compute Gaussian cls for lognormal fields
+gls = glass.fields.lognormal_gls(cls)
 
 # generator for lognormal matter fields
 matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)
diff --git a/examples/advanced/plot_shears.py b/examples/advanced/plot_shears.py
index 8ed3781..2988e69 100644
--- a/examples/advanced/plot_shears.py
+++ b/examples/advanced/plot_shears.py
@@ -3,8 +3,8 @@
 ============
 
 This example simulates a galaxy catalogue with shears affected by weak lensing,
-combining the :doc:`/basic/plot_density` and :doc:`/basic/plot_lensing` examples
-with generators for the intrinsic galaxy ellipticity and the resulting shear.
+combining the :doc:`/basic/galaxies` and :doc:`/basic/lensing` examples with
+generators for the intrinsic galaxy ellipticity and the resulting shear.
 
 '''
 
@@ -31,6 +31,7 @@
 import glass.shapes
 import glass.lensing
 import glass.galaxies
+import glass.user
 from glass.core.constants import ARCMIN2_SPHERE
 
 
@@ -56,15 +57,18 @@
 ws = glass.shells.tophat_windows(zb)
 
 # load the angular matter power spectra previously computed with CAMB
-cls = np.load('../basic/cls.npy')
+cls = glass.user.load_cls('../basic/cls.npz')
+
+# angular discretisation with 3 correlated shells
+# putting nside here means that the HEALPix pixel window function is applied
+cls = glass.fields.discretized_cls(cls, nside=nside, lmax=lmax, ncorr=3)
 
 # %%
 # Matter
 # ------
 
-# compute Gaussian cls for lognormal fields for 3 correlated shells
-# putting nside here means that the HEALPix pixel window function is applied
-gls = glass.fields.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)
+# compute Gaussian cls for lognormal fields
+gls = glass.fields.lognormal_gls(cls)
 
 # generator for lognormal matter fields
 matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)
@@ -136,7 +140,7 @@
 
         # apply the shear fields to the ellipticities
         gal_she = glass.galaxies.galaxy_shear(gal_lon, gal_lat, gal_eps,
-                                            kappa_i, gamm1_i, gamm2_i)
+                                              kappa_i, gamm1_i, gamm2_i)
 
         # map the galaxy shears to a HEALPix map; this is opaque but works
         gal_pix = hp.ang2pix(nside, gal_lon, gal_lat, lonlat=True)
diff --git a/examples/basic/biased_lognormal.py b/examples/basic/biased_lognormal.py
new file mode 100644
index 0000000..7e4dc20
--- /dev/null
+++ b/examples/basic/biased_lognormal.py
@@ -0,0 +1,184 @@
+"""
+Biased lognormal fields
+=======================
+
+This example shows how to compute a biased lognormal field by rescaling the
+*alms* of the matter field with a biased angular power spectrum.
+
+**WARNING**: This model has a major caveat. See below.
+
+"""
+
+# %%
+# Setup
+# -----
+#
+# This example simulates a lognormal matter field and a biased tracer field.
+#
+# The precomputed angular matter power spectra from the :doc:`/basic/shells`
+# example are used, so the simulation is set up in the same way.
+
+import numpy as np
+import healpy as hp
+import matplotlib.pyplot as plt
+
+# use the CAMB cosmology that generated the matter power spectra
+import camb
+from cosmology import Cosmology
+
+# GLASS imports
+from glass.shells import distance_grid, partition, tophat_windows
+from glass.fields import (discretized_cls, lognormal_gls, generate_alms,
+                          biased_cls, rescaled_alm, alm_to_lognormal, getcl)
+from glass.user import load_cls
+
+
+# cosmology for the simulation
+h = 0.7
+Oc = 0.25
+Ob = 0.05
+
+# basic parameters of the simulation
+nside = lmax = 256
+
+# set up CAMB parameters for matter angular power spectrum
+pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,
+                       NonLinear=camb.model.NonLinear_both)
+
+# get the cosmology from CAMB
+cosmo = Cosmology.from_camb(pars)
+
+# shells of 200 Mpc in comoving distance spacing
+zgrid = distance_grid(cosmo, 0., 1., dx=200.)
+
+# uniform matter weight function
+shells = tophat_windows(zgrid)
+
+# load the previously-computed angular matter power spectra
+cls = load_cls("cls.npz")
+
+# apply angular discretisation to cls
+cls = discretized_cls(cls, nside=nside, lmax=lmax, ncorr=3)
+
+# %%
+# Matter
+# ------
+#
+# The matter fields in this example are lognormal.  To simulate both matter and
+# the biased tracer field, we require the Gaussian *alm* array for each shell.
+# We therefore use the :func:`~glass.fields.generate_alms` generator, instead
+# of generating the lognormal matter fields directly.
+
+# compute Gaussian angular power spectra for lognormal matter fields
+gls = lognormal_gls(cls)
+
+# generator for Gaussian random field alms
+alms = generate_alms(gls, ncorr=3)
+
+# %%
+# Biased tracer
+# -------------
+#
+# The biased tracer fields in this example have a single bias factor *bias*,
+# but this could be an array with a different value for each shell.  The
+# biasing works by first applying the bias factors to the original *cls* array
+# using the :func:`~glass.fields.biased_cls` function.  We then compute the
+# Gaussian angular power spectra for biased lognormal fields.
+#
+# This is just one example of how to obtained biased tracer fields.  The
+# Gaussian angular power spectra could be obtained in any number of ways.
+
+# this is the bias we want to implement
+bias = 1.8
+
+# apply linear bias factor to the cls
+cls_b = biased_cls(cls, bias)
+
+# compute Gaussian cls for the biased lognormal galaxy fields
+gls_b = lognormal_gls(cls_b)
+
+# %%
+# This example will produce an integrated density for the unbiased and biased
+# fields.  Here, we are making an example *dndz* distribution, and computing
+# its contribution to each matter shell.
+
+# localised redshift distribution
+# the actual density per arcmin2 does not matter here, it is never used
+z = np.linspace(0, 1, 101)
+dndz = np.exp(-(z - 0.5)**2/(0.1)**2)
+
+# partition the dndz function into a ngal value for each shell
+ngal = partition(z, dndz, shells)
+
+# %%
+# Simulation
+# ----------
+# The simulation generates the Gaussian *alm* for each matter shell and
+# computes its lognormal matter field.  The *alm* array is rescaled with the
+# biased *gls* to compute the biased lognormal tracer field.  Both are then
+# added to an integrated number density field, using the shell's contribution
+# *ngal* to the total redshift distribution.
+
+# initialise maps for biased and unbiased total number densities
+ntot = np.zeros(12 * nside**2)
+ntot_b = np.zeros(12 * nside**2)
+
+# main loop to simulate the matter fields iteratively
+for i, alm in enumerate(alms):
+
+    # first, compute the lognormal matter field
+    delta = alm_to_lognormal(alm, nside)
+
+    # now rescale the alm to the biased lognormal field for galaxies
+    alm_b = rescaled_alm(alm, getcl(gls_b, i), getcl(gls, i))
+
+    # compute the biased lognormal galaxies field
+    delta_b = alm_to_lognormal(alm_b, nside)
+
+    # add both fields to number density
+    ntot += ngal[i] * (1 + delta)
+    ntot_b += ngal[i] * (1 + delta_b)
+
+# %%
+# Analysis
+# --------
+# To make sure the biased tracer field works as expected, we compute the
+# density contrast of the unbiased and biased integrated fields, then compare
+# their bias to the linear bias factor we used above.
+
+# this is the mean number density over the entire distribution
+nbar = np.trapz(dndz, z)
+
+# compute the integrated density contrasts
+delta_tot = (ntot - nbar)/nbar
+delta_tot_b = (ntot_b - nbar)/nbar
+
+# get the angular power spectra of the number density maps
+sim_cls = hp.anafast([delta_tot, delta_tot_b],
+                     pol=False, lmax=lmax, use_pixel_weights=True)
+
+# plot the unbiased and biased cls
+l = np.arange(lmax+1)
+fig, ax = plt.subplots(2, 1, sharex=True, sharey=True, layout="constrained")
+ax[0].plot(l, (sim_cls[1]/sim_cls[0])**0.5, "-k", lw=2, label="simulated bias (biased $\\times$ biased)")
+ax[0].axhline(bias, c="r", ls="-", lw=1, label="bias factor")
+ax[0].legend(frameon=False)
+ax[0].set_ylabel(r"bias factor")
+ax[1].plot(l, sim_cls[2]/sim_cls[0], "-k", lw=2, label="simulated bias (biased $\\times$ unbiased)")
+ax[1].axhline(bias, c="r", ls="-", lw=1, label="bias factor")
+ax[1].legend(frameon=False)
+ax[1].set_ylabel(r"bias factor")
+ax[1].set_xscale("log")
+ax[1].set_xlabel(r"angular mode number $l$")
+ax[1].set_ylim(bias - 0.25, bias + 0.25)
+plt.show()
+
+# %%
+# The results show that the biased lognormal tracer is able to capture the
+# linear bias in the auto-correlation of the biased tracer reasonably well.
+# **However, there is a significant deviation from the linear model in the
+# cross-correlation between biased tracer and unbiased field.**  The "defect"
+# in the model is roughly a function of the amplitude of the spectrum, and
+# hence tendentially worse for lower redshifts, higher :math:`sigma_8` values,
+# and higher bias values.  The biased lognormal tracer model should therefore
+# only be used where this defect is acceptable.
diff --git a/examples/basic/cls.npy b/examples/basic/cls.npy
deleted file mode 100644
index f7b3db9..0000000
Binary files a/examples/basic/cls.npy and /dev/null differ
diff --git a/examples/basic/cls.npz b/examples/basic/cls.npz
new file mode 100644
index 0000000..efaf6f0
Binary files /dev/null and b/examples/basic/cls.npz differ
diff --git a/examples/basic/plot_density.py b/examples/basic/galaxies.py
similarity index 93%
rename from examples/basic/plot_density.py
rename to examples/basic/galaxies.py
index 190d84d..1ecc96c 100644
--- a/examples/basic/plot_density.py
+++ b/examples/basic/galaxies.py
@@ -28,6 +28,7 @@
 import glass.fields
 import glass.points
 import glass.galaxies
+import glass.user
 
 
 # cosmology for the simulation
@@ -52,14 +53,17 @@
 ws = glass.shells.tophat_windows(zb)
 
 # load the angular matter power spectra previously computed with CAMB
-cls = np.load('cls.npy')
+cls = glass.user.load_cls('cls.npz')
+
+# angular discretisation with 3 correlated shells
+cls = glass.fields.discretized_cls(cls, nside=nside, lmax=lmax, ncorr=3)
 
 # %%
 # Matter
 # ------
 
-# compute Gaussian cls for lognormal fields for 3 correlated shells
-gls = glass.fields.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)
+# compute Gaussian cls for lognormal fields
+gls = glass.fields.lognormal_gls(cls)
 
 # generator for lognormal matter fields
 matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)
diff --git a/examples/basic/plot_lensing.py b/examples/basic/lensing.py
similarity index 94%
rename from examples/basic/plot_lensing.py
rename to examples/basic/lensing.py
index 7489639..4190fe0 100644
--- a/examples/basic/plot_lensing.py
+++ b/examples/basic/lensing.py
@@ -36,6 +36,7 @@
 import glass.fields
 import glass.lensing
 import glass.galaxies
+import glass.user
 
 
 # cosmology for the simulation
@@ -60,15 +61,18 @@
 ws = glass.shells.tophat_windows(zb)
 
 # load the angular matter power spectra previously computed with CAMB
-cls = np.load('cls.npy')
+cls = glass.user.load_cls('cls.npz')
+
+# angular discretisation with 3 correlated shells
+# putting nside here means that the HEALPix pixel window function is applied
+cls = glass.fields.discretized_cls(cls, nside=nside, lmax=lmax, ncorr=3)
 
 # %%
 # Matter
 # ------
 
-# compute Gaussian cls for lognormal fields for 3 correlated shells
-# putting nside here means that the HEALPix pixel window function is applied
-gls = glass.fields.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)
+# compute Gaussian cls for lognormal fields
+gls = glass.fields.lognormal_gls(cls)
 
 # generator for lognormal matter fields
 matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)
diff --git a/examples/basic/matter.py b/examples/basic/matter.py
index b6a52d2..8d4cfa5 100644
--- a/examples/basic/matter.py
+++ b/examples/basic/matter.py
@@ -26,6 +26,7 @@
 # these are the GLASS imports: matter and random fields
 import glass.shells
 import glass.fields
+import glass.user
 
 
 # cosmology for the simulation
@@ -48,10 +49,13 @@
 zb = glass.shells.distance_grid(cosmo, 0., 1., dx=200.)
 
 # load precomputed angular matter power spectra
-cls = np.load('cls.npy')
+cls = glass.user.load_cls('cls.npz')
 
-# compute Gaussian cls for lognormal fields with 3 correlated shells
-gls = glass.fields.lognormal_gls(cls, ncorr=3, nside=nside)
+# angular discretisation with 3 correlated shells
+cls = glass.fields.discretized_cls(cls, nside=nside, ncorr=3)
+
+# compute Gaussian cls for lognormal fields
+gls = glass.fields.lognormal_gls(cls)
 
 # this generator will yield the matter fields in each shell
 matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)
diff --git a/examples/basic/shells.py b/examples/basic/shells.py
index 5d4e6aa..630cddd 100644
--- a/examples/basic/shells.py
+++ b/examples/basic/shells.py
@@ -17,12 +17,12 @@
 # Here we define the shells for these examples, and use CAMB to compute the
 # angular matter power spectra for the shell definitions.
 
-import numpy as np
 import camb
 from cosmology import Cosmology
 
-import glass.shells
-import glass.ext.camb
+from glass.shells import distance_grid, tophat_windows
+from glass.ext.camb import matter_cls, camb_tophat_weight
+from glass.user import save_cls
 
 
 # cosmology for the simulation
@@ -41,18 +41,14 @@
 cosmo = Cosmology.from_camb(pars)
 
 # shells of 200 Mpc in comoving distance spacing
-zb = glass.shells.distance_grid(cosmo, 0., 1., dx=200.)
+zgrid = distance_grid(cosmo, 0., 1., dx=200.)
 
 # uniform matter weight function
 # CAMB requires linear ramp for low redshifts
-ws = glass.shells.tophat_windows(zb, weight=glass.ext.camb.camb_tophat_weight)
+shells = tophat_windows(zgrid, weight=camb_tophat_weight)
 
 # compute angular matter power spectra with CAMB
-cls = glass.ext.camb.matter_cls(pars, lmax, ws)
+cls = matter_cls(pars, lmax, shells)
 
-
-# %%
-# Save
-# ----
-# We save the shell definitions to file, for use in other examples.
-np.save('cls.npy', cls)
+# save the angular power spectra to file, to use in other examples
+save_cls("cls.npz", cls)