LSS projected basics (#140)

* Move files

* Move files

* Add galactic to RA/dec

* Update sky maps

* Add projected fields

* Remove old import

* Quick update

* Add IO

* Add imports

* Update imports

* Add basic file

csiborgtools/__init__.py

@@ -22,7 +22,7 @@ from .utils import (center_of_mass, delta2ncells, number_counts,
                     radec_to_cartesian, cartesian_to_radec,                     # noqa
                     thin_samples_by_acl, BIC_AIC, radec_to_galactic,            # noqa
                     heliocentric_to_cmb, calculate_acl, harmonic_evidence,      # noqa
-                    dict_samples_to_array)                                      # noqa
+                    dict_samples_to_array, galactic_to_radec)                   # noqa
 from .params import (paths_glamdring, simname2boxsize, simname2Omega_m,         # noqa
                      snap2redshift)                                             # noqa
 

csiborgtools/cmb_correlation/__init__.py

@@ -0,0 +1,16 @@
+# Copyright (C) 2024 Richard Stiskalek
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 3 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+from .io import read_CMB_temperature                                            # noqa

csiborgtools/cmb_correlation/io.py

@@ -0,0 +1,59 @@
+# Copyright (C) 2022 Richard Stiskalek
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 3 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+"""
+Various I/O functions for reading and writing data.
+"""
+from warnings import warn
+
+import healpy as hp
+from astropy.io import fits
+
+
+def read_CMB_temperature(fname=None, nside_out=None,
+                         convert_to_ring_ordering=True, normalize=False,
+                         verbose=True):
+    """
+    Read the CMB temperature map from a FITS file.
+    """
+    if fname is None:
+        warn("Using the glamdrnig path to the default temperature map.",
+             UserWarning)
+        fname = "/mnt/extraspace/rstiskalek/catalogs/CMB/COM_CMB_IQU-smica_2048_R3.00_full.fits"  # noqa
+
+    f = fits.open(fname)
+    if verbose:
+        print(f"Reading CMB temperature map from `{fname}`.")
+
+    skymap = f[1].data["I_STOKES"]
+    mask = f[1].data["TMASK"]
+    f.close()
+
+    if nside_out is not None:
+        if verbose:
+            print(f"Moving to nside = {nside_out}...")
+        skymap = hp.pixelfunc.ud_grade(skymap, nside_out, order_in="NESTED")
+        mask = hp.pixelfunc.ud_grade(mask, nside_out, order_in="NESTED")
+
+    if convert_to_ring_ordering:
+        if verbose:
+            print("Converting to RING ordering...")
+        skymap = hp.reorder(skymap, n2r=True)
+        mask = hp.reorder(mask, n2r=True)
+
+    if normalize:
+        skymap -= skymap.mean()
+        skymap /= skymap.std()
+
+    return skymap, mask

csiborgtools/field/__init__.py

@@ -19,8 +19,8 @@ try:
     from .density import (DensityField, PotentialField, TidalTensorField,       # noqa
                           VelocityField, radial_velocity, power_spectrum,       # noqa
                           overdensity_field)                                    # noqa
-    from .interp import (evaluate_cartesian_cic, evaluate_sky, evaluate_los,    # noqa
+    from .interp import (evaluate_cartesian_cic, evaluate_los, field2rsp,       # noqa
-                         field2rsp, fill_outside, make_sky,                     # noqa
+                         fill_outside, make_sky,                                # noqa
                          observer_peculiar_velocity, smoothen_field,            # noqa
                          field_at_distance)                                     # noqa
 except ImportError:

csiborgtools/field/interp.py

@@ -20,7 +20,7 @@ import numpy as np
 import smoothing_library as SL
 from numba import jit
 from scipy.interpolate import RegularGridInterpolator
-from tqdm import tqdm, trange
+from tqdm import tqdm
 
 from ..utils import periodic_wrap_grid, radec_to_cartesian
 from .utils import divide_nonzero, force_single_precision, nside2radec
@@ -303,56 +303,13 @@ def evaluate_los(*fields, sky_pos, boxsize, rmax, dr, smooth_scales=None,
 ###############################################################################
 
 
-def evaluate_sky(*fields, pos, mpc2box, smooth_scales=None, verbose=False):
+def make_sky(field, angpos, rmax, dr, boxsize, interpolation_method="cic",
-    """
+             return_full=False, verbose=True):
-    Evaluate a scalar field(s) at radial distance `Mpc / h`, right ascensions
-    [0, 360) deg and declinations [-90, 90] deg. Uses CIC interpolation.
-
-    Parameters
-    ----------
-    fields : (list of) 3-dimensional array of shape `(grid, grid, grid)`
-        Field to be interpolated.
-    pos : 2-dimensional array of shape `(n_samples, 3)`
-        Query spherical coordinates.
-    mpc2box : float
-        Conversion factor to multiply the radial distance by to get box units.
-    smooth_scales : (list of) float, optional
-        Smoothing scales in `Mpc / h`. If `None`, no smoothing is performed.
-    verbose : bool, optional
-        Smoothing verbosity flag.
-
-    Returns
-    -------
-    (list of) 1-dimensional array of shape `(n_samples, len(smooth_scales))`
-    """
-    # Make a copy of the positions to avoid modifying the input.
-    pos = np.copy(pos)
-
-    pos = force_single_precision(pos)
-    pos[:, 0] *= mpc2box
-
-    cart_pos = radec_to_cartesian(pos) + 0.5
-
-    if smooth_scales is not None:
-        if isinstance(smooth_scales, (int, float)):
-            smooth_scales = [smooth_scales]
-
-        if isinstance(smooth_scales, list):
-            smooth_scales = np.array(smooth_scales, dtype=np.float32)
-
-        smooth_scales *= mpc2box
-
-    return evaluate_cartesian_cic(*fields, pos=cart_pos,
-                                  smooth_scales=smooth_scales,
-                                  verbose=verbose)
-
-
-def make_sky(field, angpos, dist, boxsize, verbose=True):
     r"""
     Make a sky map of a scalar field. The observer is in the centre of the
     box the field is evaluated along directions `angpos` (RA [0, 360) deg,
-    dec [-90, 90] deg). Along each direction, the field is evaluated distances
+    dec [-90, 90] deg). The field is evaluated up to `rmax` with a linear
-    `dist` (`Mpc / h`) and summed. Uses CIC interpolation.
+    spacing of `dr` in `Mpc / h`.
 
     Parameters
     ----------
@@ -360,39 +317,38 @@ def make_sky(field, angpos, dist, boxsize, verbose=True):
         Field to be interpolated
     angpos : 2-dimensional arrays of shape `(ndir, 2)`
         Directions to evaluate the field.
-    dist : 1-dimensional array
+    rmax : float
-        Uniformly spaced radial distances to evaluate the field in `Mpc / h`.
+        Maximum radial distance in `Mpc / h`.
+    dr : float
+        Radial distance step in `Mpc / h`.
     boxsize : float
         Box size in `Mpc / h`.
+    interpolation_method : str, optional
+        Interpolation method. Must be one of `cic` or one of the methods of
+        `scipy.interpolate.RegularGridInterpolator`.
+    return_full : bool, optional
+        If `True`, return the full interpolated field instead of the average
+        field at each radial distance.
     verbose : bool, optional
         Verbosity flag.
 
     Returns
     -------
-    interp_field : 1-dimensional array of shape `(n_pos, )`.
+    interp_field : 1-dimensional array of shape `(n_pos, )`
     """
-    dx = dist[1] - dist[0]
+    rdist, finterp = evaluate_los(
-    assert np.allclose(dist[1:] - dist[:-1], dx)
+        field, sky_pos=angpos, boxsize=boxsize, rmax=rmax, dr=dr,
-    assert angpos.ndim == 2 and dist.ndim == 1
+        smooth_scales=None, verbose=verbose,
+        interpolation_method=interpolation_method)
 
-    # We loop over the angular directions, at each step evaluating a vector
+    if return_full:
-    # of distances. We pre-allocate arrays for speed.
+        return rdist, finterp
-    dir_loop = np.full((dist.size, 3), np.nan, dtype=np.float32)
 
-    ndir = angpos.shape[0]
+    finterp *= rdist**2
-    out = np.full(ndir, np.nan, dtype=np.float32)
+    finterp = np.trapz(finterp, x=rdist, axis=1)
-    for i in trange(ndir) if verbose else range(ndir):
+    finterp /= np.trapz(rdist**2, x=rdist)
-        dir_loop[:, 0] = dist
-        dir_loop[:, 1] = angpos[i, 0]
-        dir_loop[:, 2] = angpos[i, 1]
 
-        out[i] = np.sum(
+    return finterp
-            dist**2 * evaluate_sky(field, pos=dir_loop, mpc2box=1 / boxsize))
-
-    # Assuming the field is in h^2 Msun / kpc**3, we need to convert Mpc / h
-    # to kpc / h and multiply by the pixel area.
-    out *= dx * 1e9 * 4 * np.pi / len(angpos)
-    return out
 
 
 ###############################################################################

csiborgtools/read/paths.py

@@ -634,3 +634,24 @@ class Paths:
         fdir = join(self.postdir, "field_los")
         try_create_directory(fdir)
         return join(fdir, f"los_{catalogue}_{simnname}.hdf5")
+
+    def field_projected(self, simname, kind):
+        """
+        Path to the files containing the projected fields on the sky.
+
+        Parameters
+        ----------
+        simname : str
+            Simulation name.
+        kind : str
+            Field type.
+
+        Returns
+        -------
+        str
+        """
+        fdir = join(self.postdir, "field_projected")
+        try_create_directory(fdir)
+        return join(fdir, f"{simname}_{kind}_volume_weighted.hdf5")
+
+

csiborgtools/summary/__init__.py

@@ -13,6 +13,7 @@
 # with this program; if not, write to the Free Software Foundation, Inc.,
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
+from .cmb_correlation import read_projected_matter                              # noqa
 from .knn_summary import kNNCDFReader                                           # noqa
 from .nearest_neighbour_summary import NearestNeighbourReader                   # noqa
 from .overlap_summary import weighted_stats                                     # noqa

csiborgtools/summary/cmb_correlation.py

@@ -0,0 +1,80 @@
+# Copyright (C) 2024 Richard Stiskalek
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 3 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+"""
+Useful functions for getting summary data for CMB x MATTER cross-correlation.
+"""
+import numpy as np
+from h5py import File
+from healpy.sphtfunc import smoothing
+import healpy as hp
+from tqdm import tqdm
+
+
+def read_projected_matter(simname, paths, fwhm_deg=None, remove_monopole=False,
+                          remove_dipole=False, normalize=False):
+    """
+    Read the projected matter density field for a given simulation.
+
+    Parameters
+    ----------
+    simname : str
+        The name of the simulation.
+    paths : csiborgtools.read.Paths
+        Paths object.
+    fwhm_deg : float, optional
+        The full-width at half-maximum of the smoothing kernel in degrees.
+    remove_monopole : bool, optional
+        Whether to remove the monopole from the field.
+    remove_dipole : bool, optional
+        Whether to remove the dipole from the field.
+    normalize : bool, optional
+        Whether to apply standard normalization to the field.
+
+    Returns
+    -------
+    dist_ranges : 2-dimensional array of shape (n_dist_ranges, 2)
+        The distance ranges for the field.
+    data : 3-dimensional array of shape (n_sims, n_dist_ranges, npix)
+        The projected matter density field.
+    """
+    kind = "density"
+    nsims = paths.get_ics(simname)
+
+    fname = paths.field_projected(simname, kind)
+    with File(fname, "r") as f:
+        dist_ranges = f["dist_ranges"][...]
+
+        npix = len(f[f"nsim_{nsims[0]}/dist_range_0"])
+
+        data = np.zeros((len(nsims), len(dist_ranges), npix))
+        for i, nsim in enumerate(tqdm(nsims, desc="Simulations")):
+            for j in range(len(dist_ranges)):
+                skymap = f[f"nsim_{nsim}/dist_range_{j}"][...]
+
+                if fwhm_deg is not None:
+                    skymap = smoothing(skymap, fwhm=fwhm_deg * np.pi / 180.0)
+
+                if remove_monopole:
+                    hp.pixelfunc.remove_monopole(skymap, copy=False)
+                if remove_dipole:
+                    hp.pixelfunc.remove_dipole(skymap, copy=False)
+
+                if normalize:
+                    skymap -= np.mean(skymap)
+                    skymap /= np.std(skymap)
+
+                data[i, j] = skymap
+
+    return dist_ranges, data

csiborgtools/utils.py

@@ -163,6 +163,12 @@ def radec_to_galactic(ra, dec):
     return c.galactic.l.degree, c.galactic.b.degree
 
 
+def galactic_to_radec(l, b):  # noqa
+    """Convert galactic coordinates to right ascension and declination."""
+    c = SkyCoord(l=l*u.degree, b=b*u.degree, frame='galactic')
+    return c.icrs.ra.degree, c.icrs.dec.degree
+
+
 def radec_to_supergalactic(ra, dec):
     """Convert right ascension and declination to supergalactic coordinates."""
     c = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')

scripts/field_prop/field_projected.py

@@ -0,0 +1,127 @@
+# Copyright (C) 2023 Richard Stiskalek
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 3 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+"""
+A script to calculate the projected density field for a given simulation as a
+sky map. The script is not parallelized in any way. The generated fields are
+converted to galactic coordinates to match the CMB maps.
+"""
+from argparse import ArgumentParser
+from os import remove
+from os.path import exists, join
+
+import csiborgtools
+import numpy as np
+from h5py import File
+
+
+def get_field(simname, nsim, field_kind, MAS, grid):
+    """Get the appropriate field reader for the simulation."""
+    if simname == "csiborg1":
+        reader = csiborgtools.read.CSiBORG1Field(nsim)
+    elif "csiborg2" in simname:
+        kind = simname.split("_")[-1]
+        reader = csiborgtools.read.CSiBORG2Field(nsim, kind)
+    elif simname == "csiborg2X":
+        reader = csiborgtools.read.CSiBORG2XField(nsim)
+    elif "quijote" in simname:
+        reader = csiborgtools.read.QuijoteField(nsim)
+    else:
+        raise ValueError(f"Unknown simname: `{simname}`.")
+
+    if field_kind == "density":
+        return reader.density_field(MAS, grid)
+    else:
+        raise ValueError(f"Unknown field kind: `{field_kind}`.")
+
+
+def main(simname, nsims, field_kind, nside, dist_ranges, MAS, grid,
+         volume_weight, folder, normalize_to_overdensity=True):
+    boxsize = csiborgtools.simname2boxsize(simname)
+    Om0 = csiborgtools.simname2Omega_m(simname)
+    matter_density = Om0 * 277.53662724583074  # Msun / kpc^3
+
+    fname = join(folder, f"{simname}_{field_kind}.hdf5")
+    if volume_weight:
+        fname = fname.replace(".hdf5", "_volume_weighted.hdf5")
+    print(f"Writing to `{fname}`...")
+    if exists(fname):
+        remove(fname)
+
+    with File(fname, "w") as f:
+        f.create_dataset("dist_ranges", data=np.asarray(dist_ranges))
+        f.create_dataset("nsims", data=nsims)
+
+    # These are at first generated in RA/dec but we can assume it is galactic
+    # and convert it to RA/dec.
+    pixel_angpos = csiborgtools.field.nside2radec(nside)
+    RA, dec = csiborgtools.galactic_to_radec(*pixel_angpos.T)
+    pixel_angpos = np.vstack([RA, dec]).T
+    npix = len(pixel_angpos)
+
+    Rmax = np.asanyarray(dist_ranges).reshape(-1).max()
+    dr = 0.5 * boxsize / grid
+    print(f"{'R_max:':<20} {Rmax} Mpc / h", flush=True)
+    print(f"{'dr:':<20} {dr} Mpc / h", flush=True)
+
+    for nsim in nsims:
+        print(f"Interpolating at {npix} pixel for simulation {nsim}...",
+              flush=True)
+
+        field = get_field(simname, nsim, field_kind, MAS, grid)
+        rdist, finterp = csiborgtools.field.make_sky(
+            field, pixel_angpos, Rmax, dr, boxsize, return_full=True,
+            interpolation_method="linear")
+
+        with File(fname, "a") as f:
+            grp = f.create_group(f"nsim_{nsim}")
+
+            for n in range(len(dist_ranges)):
+                dmin, dmax = dist_ranges[n]
+                k_start = np.searchsorted(rdist, dmin)
+                k_end = np.searchsorted(rdist, dmax)
+
+                r = rdist[k_start:k_end + 1]
+                y = r**2 * finterp[:, k_start:k_end + 1]
+                skymap = np.trapz(y, r, axis=-1) / np.trapz(r**2, r)
+
+                if normalize_to_overdensity:
+                    skymap /= matter_density
+                    skymap -= 1
+
+                    grp.create_dataset(f"dist_range_{n}", data=skymap)
+
+
+if __name__ == "__main__":
+    parser = ArgumentParser()
+    parser.add_argument("--simname", type=str, help="Simulation name.")
+    args = parser.parse_args()
+
+    fdir = "/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_projected"
+    dx = 25
+    dist_ranges = [[0, n * dx] for n in range(1, 5)]
+    dist_ranges += [[n * dx, (n + 1) * dx] for n in range(0, 5)]
+
+    paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
+    nsims = paths.get_ics(args.simname)
+    print(f"{'Num. sims:':<20} {len(nsims)}", flush=True)
+
+    MAS = "SPH"
+    grid = 1024
+    nside = 128
+    field_kind = "density"
+    volume_weight = True
+
+    main(args.simname, nsims, field_kind, nside, dist_ranges, MAS, grid,
+         volume_weight, fdir)

scripts/field_prop/field_projected.sh

@@ -0,0 +1,26 @@
+nthreads=1
+memory=64
+on_login=${1}
+simname=${2}
+queue="berg"
+env="/mnt/users/rstiskalek/csiborgtools/venv_csiborg/bin/python"
+file="field_projected.py"
+
+
+if [ "$on_login" != "1" ] && [ "$on_login" != "0" ]
+then
+    echo "'on_login' (1) must be either 0 or 1."
+    exit 1
+fi
+
+pythoncm="$env $file --simname $simname"
+if [ $on_login -eq 1 ]; then
+    echo $pythoncm
+    $pythoncm
+else
+    cm="addqueue -q $queue -n $nthreads -m $memory $pythoncm"
+    echo "Submitting:"
+    echo $cm
+    echo
+    eval $cm
+fi

scripts/vrad_covariance.py

@@ -1,158 +0,0 @@
-# Copyright (C) 2024 Richard Stiskalek
-# This program is free software; you can redistribute it and/or modify it
-# under the terms of the GNU General Public License as published by the
-# Free Software Foundation; either version 3 of the License, or (at your
-# option) any later version.
-#
-# This program is distributed in the hope that it will be useful, but
-# WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
-# Public License for more details.
-#
-# You should have received a copy of the GNU General Public License along
-# with this program; if not, write to the Free Software Foundation, Inc.,
-# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
-"""
-NOTE: The script is far from finished or written well.
-"""
-from os.path import join
-import numpy as np
-
-import csiborgtools
-from h5py import File
-
-
-
-from sklearn.neighbors import NearestNeighbors
-from numba import jit
-from scipy.stats import binned_statistic
-
-
-###
-
-def find_indxs(rtree, r, radius):
-    """
-    Find the indices of points that are within a given radius of a given
-    point `r`.
-    """
-    if isinstance(r, (int, float)):
-        r = np.array(r)
-
-    return rtree.radius_neighbors(
-        r.reshape(-1, 1), radius=radius, return_distance=False, )[0]
-
-@jit(nopython=True)
-def dot_product_norm(x1_norm, x2_norm):
-    """Dot product of two normalised 1D vectors."""
-    return x1_norm[0] * x2_norm[0] + x1_norm[1] * x2_norm[1] + x1_norm[2] * x2_norm[2]  # noqa
-
-
-
-@jit(nopython=True)
-def get_angdist_vrad_product(i_comb, j_comb, pos_norm, vrad, r, rbin_i, rbin_j):
-    # TODO: Add itself?
-    len_i = len(i_comb)
-    len_j = len(j_comb)
-
-    cos_angdist_values = np.full(len_i * len_j, np.nan)
-    vrad_product = np.full(len_i * len_j, np.nan)
-    weights = np.full(len_i * len_j, np.nan)
-
-    k = 0
-    for i in range(len_i):
-        pos_norm_i = pos_norm[i_comb[i]]
-        vrad_i = vrad[i_comb[i]]
-        w_i = (rbin_i / r[i_comb[i]])**2
-        for j in range(len_j):
-            cos_angdist_values[k] = dot_product_norm(pos_norm_i, pos_norm[j_comb[j]])
-
-            # Product of the peculiar velocities
-            vrad_product[k] = vrad_i * vrad[j_comb[j]]
-            # Weight the product
-            w = w_i * (rbin_j / r[j_comb[j]])**2
-            vrad_product[k] *= w
-
-            weights[k] = w
-
-            k += 1
-
-    return cos_angdist_values, vrad_product, weights
-
-
-def main(out_summed_product, out_weights, ri, rj, angular_bins, pos, vel, observer, rmax):
-    # Centre the positions at the observer
-    pos = np.copy(pos) - observer
-    r = np.linalg.norm(pos, axis=1)
-    mask = r < rmax
-
-    # Select only the haloes within the radial range
-    pos = pos[mask]
-    vel = vel[mask]
-    r = r[mask]
-
-    # Create a KDTree for the radial positions
-    rtree = NearestNeighbors().fit(r.reshape(-1, 1))
-
-    # Calculate the radial velocity and the normalised position vector
-    pos_norm = pos / r[:, None]
-    vrad = np.sum(vel * pos_norm, axis=1)
-
-    # TODO: eventually here loop over the radii
-    # for ....
-    dr = 2.5
-    i_indxs = find_indxs(rtree, ri, radius=dr)
-    j_indxs = find_indxs(rtree, rj, radius=dr)
-
-    # Calculate the cosine of the angular distance and the product of the
-    # radial velocities for each pair of points.
-    cos_angdist, vrad_product, weights = get_angdist_vrad_product(
-        i_indxs, j_indxs, pos_norm, vrad, r, ri, rj)
-
-    out_summed_product += binned_statistic(
-        cos_angdist, vrad_product, bins=angular_bins, statistic="sum")[0]
-
-    out_weights += binned_statistic(
-        cos_angdist, weights, bins=angular_bins, statistic="sum")[0]
-
-    return out_summed_product, out_weights
-
-
-if __name__ == "__main__":
-    fdir = "/mnt/extraspace/rstiskalek/BBF/Quijote_C_ij"
-
-    rmax = 150
-    nradial = 20
-    nangular = 40
-    ncatalogue = 1
-
-    # NOTE check this
-    radial_bins = np.linspace(0, rmax, nradial + 1)
-    angular_bins = np.linspace(-1, 1, nangular + 1)
-
-    summed_product = np.zeros(nangular)
-    weights = np.zeros(nangular)
-
-    fiducial_observers = csiborgtools.read.fiducial_observers(1000, rmax)
-
-    ri = 100
-    rj = 120
-
-    # for i in trange(ncatalogue, desc="Catalogues"):
-    for i in [30]:
-        cat = csiborgtools.read.QuijoteCatalogue(i)
-
-        for j in range(len(fiducial_observers)):
-            # Loop over all the fiducial observers in this simulation
-            observer = fiducial_observers[j]
-            summed_product, weights = main(
-                summed_product, weights, ri, rj, angular_bins,
-                cat["cartesian_pos"], cat["cartesian_vel"], observer, rmax)
-
-
-    # TODO: Save
-    fname = join(fdir, "test.h5")
-    print(f"Saving to `{fname}`.")
-    with File(fname, 'w') as f:
-        f.create_dataset("summed_product", data=summed_product)
-        f.create_dataset("weights", data=weights)
-        f.create_dataset("angular_bins", data=angular_bins)