Quickstart

This walks through a complete, runnable prediction: build the inputs, choose a velocity PDF, and compute the redshift-space multipoles.

1. Inputs

You need a real-space \(\xi(r)\), pairwise velocity moments, and the constant \(aH\). Here we use simple analytic stand-ins so the example runs without data; in practice these come from a simulation measurement or theory (see Inputs & interpolation).

import numpy as np
from liulu import StreamingModel
from liulu.physics.real_space import TabulatedXi
from liulu.physics.pairwise_velocity import TabulatedMoments
from liulu.physics.velocity_pdf import NIGPDF, SkewTPDF, GaussianPDF

# tabulation grid in r [Mpc/h]
r = np.geomspace(0.5, 60.0, 40)

# real-space correlation function (power-law stand-in)
xi = 1.5 * (r / 8.0) ** -1.8

# pairwise velocity moments (stand-ins; real ones are measured/modelled)
v_r     = -300.0 * np.exp(-r / 15.0)          # mean radial infall (km/s)
sigma_r = 400.0 * np.ones_like(r)             # radial dispersion (km/s)
sigma_t = 300.0 * np.ones_like(r)             # tangential dispersion (km/s)
gamma_r = -0.3 * np.exp(-r / 10.0)            # radial skewness
kappa_r =  0.8 * np.exp(-r / 10.0)            # radial excess kurtosis
kappa_t =  0.4 * np.exp(-r / 10.0)            # tangential excess kurtosis

2. Build the physics components

xi_real = TabulatedXi(r_table=r, xi_table=xi)

moments = TabulatedMoments(
    r_table=r, v_r=v_r, sigma_r=sigma_r, sigma_t=sigma_t,
    gamma_r=gamma_r, kappa_r=kappa_r, kappa_t=kappa_t,
)

aH = 91.0   # a * H(z) in km/s/(Mpc/h); e.g. ~91 for LCDM at z=0.25

TabulatedMoments takes the user-friendly sigma_*/gamma_*/kappa_* and converts them to the canonical central moments c_ij internally (see Conventions).

3. Choose a PDF and run

model = StreamingModel(xi_real, NIGPDF(moments), aH=aH, y_max=120.0)

s = np.geomspace(1.0, 50.0, 19)
mp = model.multipoles(s, ell_max=4)   # dict {0: xi0(s), 2: xi2(s), 4: xi4(s)}

print(np.c_[s, s**2 * mp[0], s**2 * mp[2]])

multipoles returns a dict keyed by even \(\ell\); take mp[0], mp[2], mp[4] for the monopole, quadrupole, hexadecapole.

Comparing to a measured (binned) multipole? Bin-average the model.

A pair-count estimator (e.g. pycorr) reports the \(s^2\)-weighted average of \(\xi_\ell\) over each \(s\) bin, not its value at the bin centre. Point-evaluating the model at geometric bin centres over-predicts the steep monopole by \(\sim\)1.7% for typical log-spaced bins and is first-order wrong wherever \(\xi_\ell\) is curved (the Fingers-of-God quadrupole, zero crossings). Compare with

mp = model.multipoles_binned(s_edges, ell_max=4)   # same bin edges as the measurement

and tabulate the input \(\xi(r)\) at the pair-weighted bin separation (pycorr's return_sep, or liulu.numerics.binning.pair_weighted_centres). See Validation for what happens if you don't.

4. Swap the PDF

The PDF is a drop-in component — compare models by swapping the one line:

for pdf in (GaussianPDF(moments), SkewTPDF(moments), NIGPDF(moments)):
    xi2 = StreamingModel(xi_real, pdf, aH=aH, y_max=120.0).multipoles(s, ell_max=4)[2]
    print(type(pdf).__name__, (s**2 * xi2)[:5])

Other outputs

model.xi_s_smu(s, mu)            # xi^s on a (s, mu) grid
model.xi_s_2d(s_perp, s_par)     # xi^s on the (s_perp, s_par) plane
model.wedges(s, mu_edges=[(0,0.5),(0.5,1.0)])   # mu wedges

See the API reference for full signatures, and Validation for how the PDF choice plays out against \(N\)-body halos.