Add your issue here

The current positions_from_delta() function carries out a large number of steps:

• Apply bias model to given delta map.
• Convert biased delta to mean number count in each map pixel.
• Apply a visibility map.
• Compute observed number counts from a distribution (e.g., Poisson).
• Accumulate the galaxies in the number count map into batches of a given size (e.g., 1_000_000)
• Sample random positions for each galaxy within its pixel.

This structure is very rigid, and makes it difficult to implement things such as #136.

We should split this up into individual functions, so that they can be mixed and matched with different implementations:

for i, delta in enumerate(matter):
    # mean number count in each pixel
    nbar = glass.nbar_from_density(arcmin2=ngal[i], nside=nside)
    # galaxy density from linear bias model, could be a complicated model
    gal = nbar * (1 + beff[i] * delta)
    # apply visibility map for this shell, could be a complicated model
    gal *= vmap[i]
    # Poisson-sample galaxy indices from map in batches
    for ipix in glass.batched_poisson(gal, rng=rng):
        # randomly sample galaxy positions within each pixel
        lon, lat = glass.random_healpix_positions(nside, ipix)

This needs a bit of experimentation for multi-dimensional arrays, but in principle seems feasible.