Skip to content

Physics

Astrophysical content only — concrete models inherit the ABCs in liulu.physics.base. See Architecture.

Base classes

Real-space correlation

liulu.physics.base.RealSpaceCorrelation

Bases: ABC

Real-space two-point correlation function xi(r).

Source code in liulu/physics/base.py 
12
13
14
15
16
17
class RealSpaceCorrelation(ABC):
    """Real-space two-point correlation function xi(r)."""

    @abstractmethod
    def __call__(self, r: ndarray) -> ndarray:
        ...

Pairwise velocity moments

liulu.physics.base.PairwiseVelocityMoments

Bases: ABC

Radial-tangential pairwise velocity moments as functions of r.

Naming convention
  • m_ij(r): moment about the origin, i-th radial order, j-th tangential order.
  • c_ij(r): central moment, i-th radial order, j-th tangential order.

Odd tangential orders vanish by isotropy and are handled automatically.

Required for Gaussian PDF: m_10, c_20, c_02 Required for skew-t PDF: additionally c_30, c_12, c_40, c_04, c_22

Source code in liulu/physics/base.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
class PairwiseVelocityMoments(ABC):
    """Radial-tangential pairwise velocity moments as functions of r.

    Naming convention
    -----------------
    - ``m_ij(r)``: moment about the origin, i-th radial order, j-th tangential order.
    - ``c_ij(r)``: central moment, i-th radial order, j-th tangential order.

    Odd tangential orders vanish by isotropy and are handled automatically.

    Required for Gaussian PDF: m_10, c_20, c_02
    Required for skew-t PDF:   additionally c_30, c_12, c_40, c_04, c_22
    """

    # ── required (Gaussian-level) ─────────────────────────────────────

    @abstractmethod
    def m_10(self, r: ndarray) -> ndarray:
        """Mean radial pairwise velocity (km/s). Negative = infall."""
        ...

    @abstractmethod
    def c_20(self, r: ndarray) -> ndarray:
        """Radial velocity variance sigma_r^2 (km/s)^2."""
        ...

    @abstractmethod
    def c_02(self, r: ndarray) -> ndarray:
        """Tangential velocity variance per component sigma_t^2 (km/s)^2."""
        ...

    # ── optional (skew-t level, defaults assume Gaussian + independent) ──

    def c_30(self, r: ndarray) -> ndarray:
        """Third radial central moment <(v_r - <v_r>)^3>."""
        return np.zeros_like(np.atleast_1d(np.asarray(r, dtype=np.float64)))

    def c_12(self, r: ndarray) -> ndarray:
        """Cross central moment <(v_r - <v_r>) v_t^2>."""
        return np.zeros_like(np.atleast_1d(np.asarray(r, dtype=np.float64)))

    def c_40(self, r: ndarray) -> ndarray:
        """Fourth radial central moment."""
        return 3.0 * self.c_20(r) ** 2

    def c_04(self, r: ndarray) -> ndarray:
        """Fourth tangential central moment."""
        return 3.0 * self.c_02(r) ** 2

    def c_22(self, r: ndarray) -> ndarray:
        """Cross central moment <(v_r - <v_r>)^2 v_t^2>."""
        return self.c_20(r) * self.c_02(r)

    # ── dispatcher ────────────────────────────────────────────────────

    def get_moment(self, r: ndarray, r_order: int, t_order: int,
                   mode: str) -> ndarray:
        """Retrieve a moment by radial/tangential order.

        Parameters
        ----------
        r : ndarray
            Pair separation.
        r_order, t_order : int
            Radial and tangential orders.
        mode : str
            ``'m'`` for moments about the origin, ``'c'`` for central moments.
        """
        r = np.atleast_1d(np.asarray(r, dtype=np.float64))
        if t_order % 2 != 0:
            return np.zeros_like(r)
        if r_order == 0 and t_order == 0:
            return np.ones_like(r)
        if mode == "c" and r_order + t_order == 1:
            return np.zeros_like(r)
        return getattr(self, f"{mode}_{r_order}{t_order}")(r)

m_10 abstractmethod 

m_10(r: ndarray) -> ndarray

Mean radial pairwise velocity (km/s). Negative = infall.

Source code in liulu/physics/base.py 
36
37
38
39
@abstractmethod
def m_10(self, r: ndarray) -> ndarray:
    """Mean radial pairwise velocity (km/s). Negative = infall."""
    ...

c_20 abstractmethod 

c_20(r: ndarray) -> ndarray

Radial velocity variance sigma_r^2 (km/s)^2.

Source code in liulu/physics/base.py 
41
42
43
44
@abstractmethod
def c_20(self, r: ndarray) -> ndarray:
    """Radial velocity variance sigma_r^2 (km/s)^2."""
    ...

c_02 abstractmethod 

c_02(r: ndarray) -> ndarray

Tangential velocity variance per component sigma_t^2 (km/s)^2.

Source code in liulu/physics/base.py 
46
47
48
49
@abstractmethod
def c_02(self, r: ndarray) -> ndarray:
    """Tangential velocity variance per component sigma_t^2 (km/s)^2."""
    ...

c_30 

c_30(r: ndarray) -> ndarray

Third radial central moment <(v_r - )^3>.

Source code in liulu/physics/base.py 
53
54
55
def c_30(self, r: ndarray) -> ndarray:
    """Third radial central moment <(v_r - <v_r>)^3>."""
    return np.zeros_like(np.atleast_1d(np.asarray(r, dtype=np.float64)))

c_12 

c_12(r: ndarray) -> ndarray

Cross central moment <(v_r - ) v_t^2>.

Source code in liulu/physics/base.py 
57
58
59
def c_12(self, r: ndarray) -> ndarray:
    """Cross central moment <(v_r - <v_r>) v_t^2>."""
    return np.zeros_like(np.atleast_1d(np.asarray(r, dtype=np.float64)))

c_40 

c_40(r: ndarray) -> ndarray

Fourth radial central moment.

Source code in liulu/physics/base.py 
61
62
63
def c_40(self, r: ndarray) -> ndarray:
    """Fourth radial central moment."""
    return 3.0 * self.c_20(r) ** 2

c_04 

c_04(r: ndarray) -> ndarray

Fourth tangential central moment.

Source code in liulu/physics/base.py 
65
66
67
def c_04(self, r: ndarray) -> ndarray:
    """Fourth tangential central moment."""
    return 3.0 * self.c_02(r) ** 2

c_22 

c_22(r: ndarray) -> ndarray

Cross central moment <(v_r - )^2 v_t^2>.

Source code in liulu/physics/base.py 
69
70
71
def c_22(self, r: ndarray) -> ndarray:
    """Cross central moment <(v_r - <v_r>)^2 v_t^2>."""
    return self.c_20(r) * self.c_02(r)

get_moment 

get_moment(r: ndarray, r_order: int, t_order: int, mode: str) -> ndarray

Retrieve a moment by radial/tangential order.

Parameters:

Name Type Description Default
r ndarray

Pair separation.

 required
r_order int

Radial and tangential orders.

 required
t_order int

Radial and tangential orders.

 required
mode str

'm' for moments about the origin, 'c' for central moments.

 required
Source code in liulu/physics/base.py 
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
def get_moment(self, r: ndarray, r_order: int, t_order: int,
               mode: str) -> ndarray:
    """Retrieve a moment by radial/tangential order.

    Parameters
    ----------
    r : ndarray
        Pair separation.
    r_order, t_order : int
        Radial and tangential orders.
    mode : str
        ``'m'`` for moments about the origin, ``'c'`` for central moments.
    """
    r = np.atleast_1d(np.asarray(r, dtype=np.float64))
    if t_order % 2 != 0:
        return np.zeros_like(r)
    if r_order == 0 and t_order == 0:
        return np.ones_like(r)
    if mode == "c" and r_order + t_order == 1:
        return np.zeros_like(r)
    return getattr(self, f"{mode}_{r_order}{t_order}")(r)

Velocity PDF

liulu.physics.base.VelocityPDF

Bases: ABC

Line-of-sight pairwise velocity PDF P(v_los | r_perp, r_par).

Source code in liulu/physics/base.py 
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
class VelocityPDF(ABC):
    """Line-of-sight pairwise velocity PDF P(v_los | r_perp, r_par)."""

    @abstractmethod
    def __call__(self, v_los: ndarray, r_perp: ndarray,
                 r_par: ndarray) -> ndarray:
        """Evaluate the PDF.

        Parameters
        ----------
        v_los : ndarray
            Line-of-sight pairwise velocity (km/s).
        r_perp : ndarray
            Transverse pair separation (Mpc/h).
        r_par : ndarray
            Parallel pair separation (Mpc/h), non-negative.

        Returns
        -------
        ndarray
            PDF values, (km/s)^{-1}.
        """
        ...

__call__ abstractmethod 

__call__(v_los: ndarray, r_perp: ndarray, r_par: ndarray) -> ndarray

Evaluate the PDF.

Parameters:

Name Type Description Default
v_los ndarray

Line-of-sight pairwise velocity (km/s).

 required
r_perp ndarray

Transverse pair separation (Mpc/h).

 required
r_par ndarray

Parallel pair separation (Mpc/h), non-negative.

 required

Returns:

Type Description
ndarray

PDF values, (km/s)^{-1}.
Source code in liulu/physics/base.py 
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
@abstractmethod
def __call__(self, v_los: ndarray, r_perp: ndarray,
             r_par: ndarray) -> ndarray:
    """Evaluate the PDF.

    Parameters
    ----------
    v_los : ndarray
        Line-of-sight pairwise velocity (km/s).
    r_perp : ndarray
        Transverse pair separation (Mpc/h).
    r_par : ndarray
        Parallel pair separation (Mpc/h), non-negative.

    Returns
    -------
    ndarray
        PDF values, (km/s)^{-1}.
    """
    ...

Real-space correlation providers

TabulatedXi

liulu.physics.real_space.TabulatedXi

Bases: RealSpaceCorrelation

Interpolate xi(r) from tabulated data.

Parameters:

Name Type Description Default
r_table array_like

Separation values, Mpc/h.

 required
xi_table array_like

Correlation function values.

 required
interpolator callable

Interpolation factory with signature (r, xi) -> callable. Defaults to log-log interpolation (:func:numerics.loglog_interp).

 None
tail (zero, power)

Large-r extrapolation above the last tabulated point (ignored when a custom interpolator is given). "zero" (default) sets xi -> 0; "power" continues the power law of the trailing consistent-sign points -- for a table ending in the negative post-BAO regime this is a decaying negative tail (xi -> 0^-). A short table fed straight to the streaming model should use tail="power": the integral probes r well beyond the table, and a hard xi -> 0 truncation biases the large-separation quadrupole.

 "zero"
Notes

Interpolation is the hybrid :func:~liulu.numerics.interpolation.loglog_interp: linear in (log r, log xi) across intervals where both tabulated values are positive (exact for power laws), and linear in (log r, xi) across any interval touching xi <= 0 -- so tabulated negative values (BAO zero crossing, halo-exclusion dip) are honoured exactly, never bridged over or zeroed. The streaming integral needs xi bounded as r -> 0; the constant edge fill below r_min guarantees this for tabulated input.

Source code in liulu/physics/real_space.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
class TabulatedXi(RealSpaceCorrelation):
    """Interpolate xi(r) from tabulated data.

    Parameters
    ----------
    r_table : array_like
        Separation values, Mpc/h.
    xi_table : array_like
        Correlation function values.
    interpolator : callable, optional
        Interpolation factory with signature (r, xi) -> callable.
        Defaults to log-log interpolation (:func:`numerics.loglog_interp`).
    tail : {"zero", "power"}, optional
        Large-r extrapolation above the last tabulated point (ignored when a
        custom ``interpolator`` is given). ``"zero"`` (default) sets
        ``xi -> 0``; ``"power"`` continues the power law of the trailing
        consistent-sign points -- for a table ending in the negative
        post-BAO regime this is a decaying *negative* tail (``xi -> 0^-``).
        A short table fed straight to the streaming model should use
        ``tail="power"``: the integral probes r well beyond the table, and a
        hard ``xi -> 0`` truncation biases the large-separation quadrupole.

    Notes
    -----
    Interpolation is the hybrid :func:`~liulu.numerics.interpolation.loglog_interp`:
    linear in (log r, log xi) across intervals where both tabulated values are
    positive (exact for power laws), and linear in (log r, xi) across any
    interval touching ``xi <= 0`` -- so tabulated negative values (BAO zero
    crossing, halo-exclusion dip) are honoured exactly, never bridged over or
    zeroed. The streaming integral needs ``xi`` bounded as ``r -> 0``; the
    constant edge fill below ``r_min`` guarantees this for tabulated input.
    """

    def __init__(self, r_table, xi_table, interpolator=None, tail="zero"):
        if tail not in ("zero", "power"):
            raise ValueError(f"tail must be 'zero' or 'power', got {tail!r}.")
        self._r = np.asarray(r_table, dtype=np.float64)
        self._xi = np.asarray(xi_table, dtype=np.float64)

        if interpolator is not None:
            self._interp = interpolator(self._r, self._xi)
        else:
            from ..numerics.interpolation import loglog_interp
            # Hybrid power-law interpolation: log-log across positive-positive
            # intervals (near-exact for the steeply falling xi; avoids the
            # convex overestimate of log-linear, which biased the model
            # monopole ~1.4% high) and linear-in-xi across any interval with a
            # non-positive value, so measured negative xi is honoured. Constant
            # (edge) extrapolation below r_min keeps clustering on small
            # scales; above r_max, xi -> 0 (tail="zero") or a sign-aware
            # power-law continuation (tail="power").
            high = "power" if tail == "power" else 0.0
            self._interp = loglog_interp(
                self._r, self._xi, fill_value=("edge", high),
                name="TabulatedXi xi(r)")

    def __call__(self, r):
        return self._interp(r)

Linear-theory / Halofit

liulu.physics.real_space.LinearTheoryXi

Bases: RealSpaceCorrelation

xi(r) from linear theory power spectrum via Hankel transform.

Parameters:

Name Type Description Default
k array_like

Wavenumbers, h/Mpc.

 required
pk array_like

Linear power spectrum P(k), (Mpc/h)^3.

 required
hankel_transform callable

Transform function with signature (k, pk, r) -> xi. Defaults to numerics.hankel.pk_to_xi.

 None
Source code in liulu/physics/real_space.py 
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
class LinearTheoryXi(RealSpaceCorrelation):
    """xi(r) from linear theory power spectrum via Hankel transform.

    Parameters
    ----------
    k : array_like
        Wavenumbers, h/Mpc.
    pk : array_like
        Linear power spectrum P(k), (Mpc/h)^3.
    hankel_transform : callable, optional
        Transform function with signature (k, pk, r) -> xi.
        Defaults to numerics.hankel.pk_to_xi.
    """

    def __init__(self, k, pk, hankel_transform=None):
        self._k = np.asarray(k, dtype=np.float64)
        self._pk = np.asarray(pk, dtype=np.float64)

        if hankel_transform is not None:
            self._transform = hankel_transform
        else:
            from ..numerics.hankel import pk_to_xi
            self._transform = pk_to_xi

    def __call__(self, r):
        return self._transform(self._k, self._pk, r)

liulu.physics.real_space.HalofitXi

Bases: RealSpaceCorrelation

xi(r) from Halofit/HMCode nonlinear P(k) via Hankel transform.

Same interface as LinearTheoryXi but expects nonlinear P(k).

Source code in liulu/physics/real_space.py 
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
class HalofitXi(RealSpaceCorrelation):
    """xi(r) from Halofit/HMCode nonlinear P(k) via Hankel transform.

    Same interface as LinearTheoryXi but expects nonlinear P(k).
    """

    def __init__(self, k, pk_nl, hankel_transform=None):
        self._k = np.asarray(k, dtype=np.float64)
        self._pk = np.asarray(pk_nl, dtype=np.float64)

        if hankel_transform is not None:
            self._transform = hankel_transform
        else:
            from ..numerics.hankel import pk_to_xi
            self._transform = pk_to_xi

    def __call__(self, r):
        return self._transform(self._k, self._pk, r)

Pairwise velocity moments

TabulatedMoments

liulu.physics.pairwise_velocity.TabulatedMoments

Bases: PairwiseVelocityMoments

Interpolate pairwise velocity moments from tabulated data.

Parameters:

Name Type Description Default
r_table array_like

Separation values, Mpc/h.

 required
v_r array_like

Mean radial pairwise velocity, km/s (negative = infall).

 required
sigma_r array_like

Radial pairwise velocity dispersion, km/s.

 required
sigma_t array_like

Tangential pairwise velocity dispersion (per component), km/s.

 required
gamma_r array_like

Radial standardised skewness. Defaults to zero.

 None
kappa_r array_like

Radial excess kurtosis. Defaults to zero (Gaussian).

 None
kappa_t array_like

Tangential excess kurtosis. Defaults to zero (Gaussian).

 None
c_12_table array_like

Cross central moment <(v_r - ) v_t^2>. Defaults to zero.

 None
c_22_table array_like

Cross central moment <(v_r - )^2 v_t^2>. If None, defaults to sigma_r^2 * sigma_t^2 (independence assumption).

 None
interpolator callable

Interpolation factory with signature (r, y) -> callable.

 None
Source code in liulu/physics/pairwise_velocity.py 
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
class TabulatedMoments(PairwiseVelocityMoments):
    """Interpolate pairwise velocity moments from tabulated data.

    Parameters
    ----------
    r_table : array_like
        Separation values, Mpc/h.
    v_r : array_like
        Mean radial pairwise velocity, km/s (negative = infall).
    sigma_r : array_like
        Radial pairwise velocity dispersion, km/s.
    sigma_t : array_like
        Tangential pairwise velocity dispersion (per component), km/s.
    gamma_r : array_like, optional
        Radial standardised skewness. Defaults to zero.
    kappa_r : array_like, optional
        Radial excess kurtosis. Defaults to zero (Gaussian).
    kappa_t : array_like, optional
        Tangential excess kurtosis. Defaults to zero (Gaussian).
    c_12_table : array_like, optional
        Cross central moment <(v_r - <v_r>) v_t^2>. Defaults to zero.
    c_22_table : array_like, optional
        Cross central moment <(v_r - <v_r>)^2 v_t^2>. If None, defaults
        to sigma_r^2 * sigma_t^2 (independence assumption).
    interpolator : callable, optional
        Interpolation factory with signature (r, y) -> callable.
    """

    def __init__(self, r_table, v_r, sigma_r, sigma_t,
                 gamma_r=None, kappa_r=None, kappa_t=None,
                 c_12_table=None, c_22_table=None,
                 interpolator=None):
        from ..numerics.interpolation import log_linear_interp

        r = np.asarray(r_table, dtype=np.float64)
        if interpolator is not None:
            _interp = interpolator
        else:
            # Constant (edge) extrapolation: outside the table the dispersions
            # plateau to their boundary value rather than collapsing to 0
            # (which the PDF would read as "no pairs" and drop from the integral).
            def _interp(r_, y_):
                return log_linear_interp(r_, y_, fill_value="edge",
                                         name="TabulatedMoments")

        sigma_r = np.asarray(sigma_r, dtype=np.float64)
        sigma_t = np.asarray(sigma_t, dtype=np.float64)

        self._v_r = _interp(r, np.asarray(v_r, dtype=np.float64))
        self._sigma_r = _interp(r, sigma_r)
        self._sigma_t = _interp(r, sigma_t)

        if gamma_r is not None:
            _gamma_r = _interp(r, np.asarray(gamma_r, dtype=np.float64))
            self._c_30 = lambda r_: _gamma_r(r_) * self._sigma_r(r_) ** 3
        else:
            self._c_30 = None

        if kappa_r is not None:
            _kappa_r = _interp(r, np.asarray(kappa_r, dtype=np.float64))
            self._c_40 = lambda r_: (_kappa_r(r_) + 3.0) * self._sigma_r(r_) ** 4
        else:
            self._c_40 = None

        if kappa_t is not None:
            _kappa_t = _interp(r, np.asarray(kappa_t, dtype=np.float64))
            self._c_04 = lambda r_: (_kappa_t(r_) + 3.0) * self._sigma_t(r_) ** 4
        else:
            self._c_04 = None

        if c_12_table is not None:
            self._c_12 = _interp(r, np.asarray(c_12_table, dtype=np.float64))
        else:
            self._c_12 = None

        if c_22_table is not None:
            self._c_22 = _interp(r, np.asarray(c_22_table, dtype=np.float64))
        else:
            self._c_22 = None

    def m_10(self, r):
        return self._v_r(r)

    def c_20(self, r):
        return self._sigma_r(r) ** 2

    def c_02(self, r):
        return self._sigma_t(r) ** 2

    def c_30(self, r):
        if self._c_30 is not None:
            return self._c_30(r)
        return super().c_30(r)

    def c_12(self, r):
        if self._c_12 is not None:
            return self._c_12(r)
        return super().c_12(r)

    def c_40(self, r):
        if self._c_40 is not None:
            return self._c_40(r)
        return super().c_40(r)

    def c_04(self, r):
        if self._c_04 is not None:
            return self._c_04(r)
        return super().c_04(r)

    def c_22(self, r):
        if self._c_22 is not None:
            return self._c_22(r)
        return super().c_22(r)

Velocity PDFs

Gaussian

liulu.physics.velocity_pdf.GaussianPDF

Bases: VelocityPDF

Gaussian LOS pairwise velocity PDF.

Parameters:

Name Type Description Default
moments PairwiseVelocityMoments

Radial-tangential moment functions.

 required
Source code in liulu/physics/velocity_pdf.py 
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
class GaussianPDF(VelocityPDF):
    """Gaussian LOS pairwise velocity PDF.

    Parameters
    ----------
    moments : PairwiseVelocityMoments
        Radial-tangential moment functions.
    """

    def __init__(self, moments: PairwiseVelocityMoments):
        self._mean_los = project_to_los(moments, n=1, mode="m")
        self._c2_los = project_to_los(moments, n=2, mode="c")

    def __call__(self, v_los, r_perp, r_par):
        r, mu = separation_geometry(r_perp, r_par)
        mean = self._mean_los(r_perp, r_par, _r=r, _mu=mu)
        var = self._c2_los(r_perp, r_par, _r=r, _mu=mu)
        _require_positive_variance(var, "GaussianPDF")
        inv_sigma = 1.0 / np.sqrt(var)
        z = (v_los - mean) * inv_sigma
        return np.exp(-0.5 * z * z) * inv_sigma / np.sqrt(2.0 * np.pi)

Skew-t

liulu.physics.velocity_pdf.SkewTPDF

Bases: _FourMomentPDF

Skew-t LOS pairwise velocity PDF (Azzalini & Capitanio 2003).

Constructed from PairwiseVelocityMoments. Projects radial-tangential moments to the four LOS moments (mean, variance, skewness, kurtosis), then converts to skew-t parameters.

Parameters:

Name Type Description Default
moments PairwiseVelocityMoments

Radial-tangential moment functions.

 required
tabulate bool

If True (default), pre-compute (alpha, nu) on a (gamma1, gamma2) grid for fast vectorised evaluation. If False, use per-point fsolve (slow but exact).

 True
table_kw dict

Keyword arguments forwarded to _SkewTParamTable (e.g. gamma1_range, gamma2_range, n_g1, n_g2).

 {}
Source code in liulu/physics/velocity_pdf.py 
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
class SkewTPDF(_FourMomentPDF):
    """Skew-t LOS pairwise velocity PDF (Azzalini & Capitanio 2003).

    Constructed from PairwiseVelocityMoments. Projects radial-tangential
    moments to the four LOS moments (mean, variance, skewness, kurtosis),
    then converts to skew-t parameters.

    Parameters
    ----------
    moments : PairwiseVelocityMoments
        Radial-tangential moment functions.
    tabulate : bool
        If True (default), pre-compute (alpha, nu) on a (gamma1, gamma2)
        grid for fast vectorised evaluation.  If False, use per-point
        ``fsolve`` (slow but exact).
    table_kw : dict
        Keyword arguments forwarded to ``_SkewTParamTable`` (e.g.
        ``gamma1_range``, ``gamma2_range``, ``n_g1``, ``n_g2``).
    """

    def __init__(self, moments: PairwiseVelocityMoments, *,
                 tabulate: bool = True, **table_kw):
        super().__init__(moments)
        self._table = _SkewTParamTable(**table_kw) if tabulate else None

    def __call__(self, v_los, r_perp, r_par):
        mean, c2, gamma1, gamma2 = self._los_moments(r_perp, r_par)
        std = np.sqrt(c2)
        if self._table is not None:
            return self._eval_tabulated(v_los, mean, std, gamma1, gamma2)
        return self._eval_direct(v_los, mean, std, gamma1, gamma2)

    # ── fast path: tabulated ──────────────────────────────────────────

    def _eval_tabulated(self, v_los, mean, std, gamma1, gamma2):
        """Fully vectorised evaluation via pre-tabulated skew-t params."""
        v_los = np.asarray(v_los, dtype=np.float64)
        mean = np.atleast_1d(np.asarray(mean, dtype=np.float64))
        std = np.atleast_1d(np.asarray(std, dtype=np.float64))
        gamma1 = np.atleast_1d(np.asarray(gamma1, dtype=np.float64))
        gamma2 = np.atleast_1d(np.asarray(gamma2, dtype=np.float64))

        out_shape = np.broadcast_shapes(v_los.shape, mean.shape)

        alpha, nu = self._table(gamma1, gamma2)

        with np.errstate(invalid='ignore', divide='ignore', over='ignore'):
            delta = alpha / np.sqrt(1.0 + alpha ** 2)
            b = (np.sqrt(nu / np.pi)
                 * special.gamma((nu - 1.0) / 2.0)
                 / special.gamma(nu / 2.0))
            w = std / np.sqrt(nu / (nu - 2.0) - delta ** 2 * b ** 2)
            v_c = mean - w * delta * b
            result = _skewt_pdf(v_los, w, v_c, alpha, nu)

        # Gaussian fallback where the table had no valid solution
        bad = ~np.isfinite(result)
        self._count_fallback(bad)
        if np.any(bad):
            result = np.where(bad, self._gaussian(v_los, mean, std), result)

        return np.broadcast_to(result, out_shape)

    # ── slow path: direct fsolve per point ────────────────────────────

    def _eval_direct(self, v_los, mean, std, gamma1, gamma2):
        """Per-point evaluation using fsolve (exact, no interpolation)."""
        mean = np.atleast_1d(np.asarray(mean, dtype=np.float64))
        std = np.atleast_1d(np.asarray(std, dtype=np.float64))
        gamma1 = np.atleast_1d(np.asarray(gamma1, dtype=np.float64))
        gamma2 = np.atleast_1d(np.asarray(gamma2, dtype=np.float64))
        v_los = np.asarray(v_los, dtype=np.float64)

        if mean.size == 1:
            g1, g2 = gamma1.ravel()[0], gamma2.ravel()[0]
            if np.abs(g1) < 1e-8 and np.abs(g2) < 1e-4:
                inv_s = 1.0 / std.ravel()[0]
                z = (v_los - mean.ravel()[0]) * inv_s
                return np.exp(-0.5 * z * z) * inv_s / np.sqrt(2.0 * np.pi)
            w, v_c, alpha, nu = _moments_to_st_params(
                mean.ravel()[0], std.ravel()[0], g1, g2)
            return _skewt_pdf(v_los, w, v_c, alpha, nu)

        out_shape = np.broadcast_shapes(v_los.shape, mean.shape)
        v_flat = np.broadcast_to(v_los, out_shape).ravel()
        mean_flat = np.broadcast_to(mean, out_shape).ravel()
        std_flat = np.broadcast_to(std, out_shape).ravel()
        g1_flat = np.broadcast_to(gamma1, out_shape).ravel()
        g2_flat = np.broadcast_to(gamma2, out_shape).ravel()

        # Group by unique (mean, std, g1, g2) to amortise fsolve calls
        params = np.column_stack([mean_flat, std_flat, g1_flat, g2_flat])
        _, inverse, counts = np.unique(
            params, axis=0, return_inverse=True, return_counts=True,
        )
        result = np.empty_like(v_flat)
        for uid in range(len(counts)):
            mask = inverse == uid
            m, s, g1, g2 = params[mask][0]
            v_sub = v_flat[mask]
            if np.abs(g1) < 1e-8 and np.abs(g2) < 1e-4:
                inv_s = 1.0 / s
                z = (v_sub - m) * inv_s
                result[mask] = np.exp(-0.5 * z * z) * inv_s / np.sqrt(2.0 * np.pi)
            else:
                w, v_c, alpha, nu = _moments_to_st_params(m, s, g1, g2)
                result[mask] = _skewt_pdf(v_sub, w, v_c, alpha, nu)
        return result.reshape(out_shape)

Generalized hyperbolic

liulu.physics.velocity_pdf.GeneralizedHyperbolicPDF

Bases: _FourMomentPDF

Generalised hyperbolic LOS pairwise-velocity PDF (Kuruvilla & Porciani 2018). A normal variance-mean (Gaussian) mixture; strictly generalises the skew-t.

The four projected LOS moments (mean, variance, skewness, excess kurtosis) are matched to the GH parameters at class index lam:

  • lam == -0.5 (NIG) uses the analytic, fully vectorised inversion (fast; the recommended default -- see :class:NIGPDF);
  • any other scalar lam uses a vectorised damped-Newton shape solve;
  • lam may also be a callable lam(r_perp, r_par) returning a per-separation class index -- this realises the scale-dependent GHD of Kuruvilla & Porciani (2018), in which the fifth (shape) parameter drifts with pair separation. Calibrate such a lam(r) from a measured velocity PDF (e.g. via :func:calibrate_gh_lambda in scripts/calibrate_gh_lambda.py).

Points whose (skewness, excess kurtosis) fall outside the attainable region fall back to a Gaussian with the same first two moments.

Parameters:

Name Type Description Default
moments PairwiseVelocityMoments
required
lam float or callable

GH class index (default -0.5 = NIG), or lam(r_perp, r_par) -> ndarray for a scale-dependent shape.

 -0.5
Source code in liulu/physics/velocity_pdf.py 
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
class GeneralizedHyperbolicPDF(_FourMomentPDF):
    """Generalised hyperbolic LOS pairwise-velocity PDF (Kuruvilla & Porciani
    2018).  A normal variance-mean (Gaussian) mixture; strictly generalises the
    skew-t.

    The four projected LOS moments (mean, variance, skewness, excess kurtosis)
    are matched to the GH parameters at class index ``lam``:

    - ``lam == -0.5`` (NIG) uses the analytic, fully vectorised inversion
      (fast; the recommended default -- see :class:`NIGPDF`);
    - any other scalar ``lam`` uses a vectorised damped-Newton shape solve;
    - ``lam`` may also be a **callable** ``lam(r_perp, r_par)`` returning a
      per-separation class index -- this realises the *scale-dependent* GHD of
      Kuruvilla & Porciani (2018), in which the fifth (shape) parameter drifts
      with pair separation. Calibrate such a ``lam(r)`` from a measured velocity
      PDF (e.g. via :func:`calibrate_gh_lambda` in
      ``scripts/calibrate_gh_lambda.py``).

    Points whose (skewness, excess kurtosis) fall outside the attainable region
    fall back to a Gaussian with the same first two moments.

    Parameters
    ----------
    moments : PairwiseVelocityMoments
    lam : float or callable
        GH class index (default -0.5 = NIG), or ``lam(r_perp, r_par) -> ndarray``
        for a scale-dependent shape.
    """

    def __init__(self, moments: PairwiseVelocityMoments, *, lam=-0.5):
        super().__init__(moments)
        self.lam = lam if callable(lam) else float(lam)

    def __call__(self, v_los, r_perp, r_par):
        v_los = np.asarray(v_los, dtype=np.float64)
        mean, c2, g1, g2 = self._los_moments(r_perp, r_par)
        std = np.sqrt(c2)
        out_shape = np.broadcast_shapes(v_los.shape, mean.shape)

        if callable(self.lam):
            lam = np.broadcast_to(np.atleast_1d(np.asarray(self.lam(r_perp, r_par),
                                                           dtype=np.float64)), mean.shape)
        else:
            lam = self.lam

        if np.isscalar(lam) and abs(lam + 0.5) < 1e-12:
            alpha, beta, delta, mu, valid = _nig_params_from_moments(mean, c2, g1, g2)
        else:
            alpha, beta, delta, mu, valid = _gh_params_vectorized(lam, mean, c2, g1, g2)

        with np.errstate(invalid="ignore", divide="ignore",
                         over="ignore", under="ignore"):
            p = _gh_density(v_los, alpha, beta, delta, lam, mu)

        bad = np.broadcast_to((~valid) | ~np.isfinite(p) | (p < 0), out_shape)
        self._count_fallback(bad)
        if np.any(bad):
            p = np.where(bad, self._gaussian(v_los, mean, std), p)
        return np.broadcast_to(p, out_shape)

NIG

liulu.physics.velocity_pdf.NIGPDF

Bases: GeneralizedHyperbolicPDF

Normal-inverse-Gaussian LOS pairwise-velocity PDF (GH with lambda=-1/2).

Analytic, vectorised four-moment matching with Gaussian fallback; the recommended fast default among the generalised-hyperbolic family.

Source code in liulu/physics/velocity_pdf.py 
675
676
677
678
679
680
681
682
683
class NIGPDF(GeneralizedHyperbolicPDF):
    """Normal-inverse-Gaussian LOS pairwise-velocity PDF (GH with lambda=-1/2).

    Analytic, vectorised four-moment matching with Gaussian fallback; the
    recommended fast default among the generalised-hyperbolic family.
    """

    def __init__(self, moments: PairwiseVelocityMoments):
        super().__init__(moments, lam=-0.5)

Tabulated GH

liulu.physics.velocity_pdf.TabulatedGHPDF

Bases: _FourMomentPDF

Generalised-hyperbolic LOS velocity PDF with scale-dependent shape, its parameters tabulated over pair separation (Kuruvilla & Porciani 2018).

The expensive GH shape solve is done once on a coarse (r_perp, r_par) grid (the GH invariant shape (zeta, rho) is smooth in separation) and interpolated thereafter; the first two moments (location, scale) are taken from the exact local projected moments, so only the higher-order shape is interpolated. The class index lam may vary with separation -- pass a callable lam(r_perp, r_par) (e.g. the calibrated lam(r) of :func:calibrate_gh_lambda) to realise the scale-dependent GHD.

Outside the tabulation grid, or where the shape solve fails, it falls back to the analytic NIG (and then Gaussian), so it is always well defined.

Parameters:

Name Type Description Default
moments PairwiseVelocityMoments
required
lam float or callable

Class index, scalar or lam(r_perp, r_par) -> ndarray.

 -0.5
rp_max float

Extent of the tabulation grid in r_perp, r_par (Mpc/h).

 60.0
rpar_max float

Extent of the tabulation grid in r_perp, r_par (Mpc/h).

 60.0
n_rp int

Grid resolution.

 48
n_rpar int

Grid resolution.

 48
Source code in liulu/physics/velocity_pdf.py 
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
class TabulatedGHPDF(_FourMomentPDF):
    """Generalised-hyperbolic LOS velocity PDF with **scale-dependent** shape,
    its parameters tabulated over pair separation (Kuruvilla & Porciani 2018).

    The expensive GH shape solve is done once on a coarse ``(r_perp, r_par)``
    grid (the GH invariant shape ``(zeta, rho)`` is smooth in separation) and
    interpolated thereafter; the first two moments (location, scale) are taken
    from the *exact* local projected moments, so only the higher-order shape is
    interpolated. The class index ``lam`` may vary with separation -- pass a
    callable ``lam(r_perp, r_par)`` (e.g. the calibrated ``lam(r)`` of
    :func:`calibrate_gh_lambda`) to realise the scale-dependent GHD.

    Outside the tabulation grid, or where the shape solve fails, it falls back
    to the analytic NIG (and then Gaussian), so it is always well defined.

    Parameters
    ----------
    moments : PairwiseVelocityMoments
    lam : float or callable
        Class index, scalar or ``lam(r_perp, r_par) -> ndarray``.
    rp_max, rpar_max : float
        Extent of the tabulation grid in r_perp, r_par (Mpc/h).
    n_rp, n_rpar : int
        Grid resolution.
    """

    def __init__(self, moments: PairwiseVelocityMoments, *, lam=-0.5,
                 rp_max=60.0, rpar_max=150.0, n_rp=48, n_rpar=80):
        super().__init__(moments)
        self.lam = lam
        self._nig = NIGPDF(moments)

        from scipy.interpolate import RegularGridInterpolator
        rp_g = np.concatenate([[1e-3], np.geomspace(0.05, rp_max, n_rp - 1)])
        rpar_g = np.concatenate([[1e-3], np.geomspace(0.05, rpar_max, n_rpar - 1)])
        RP, RPAR = np.meshgrid(rp_g, rpar_g, indexing="ij")
        mean, c2, g1, g2 = self._los_moments(RP, RPAR)
        lam_g = (np.asarray(lam(RP, RPAR), float) if callable(lam)
                 else np.full(RP.shape, float(lam)))
        zeta = np.full(RP.shape, np.nan)
        rho = np.full(RP.shape, np.nan)
        gg1, gg2, ll = g1.ravel(), g2.ravel(), lam_g.ravel()
        c2r = c2.ravel()
        zf, rf = zeta.ravel(), rho.ravel()
        for i in range(gg1.size):
            if c2r[i] > 0 and gg2[i] > (5.0 / 3.0) * gg1[i] ** 2 + 1e-6:
                zf[i], rf[i] = _gh_shape_solve_scalar(ll[i], gg1[i], gg2[i])
        zeta = zf.reshape(RP.shape); rho = rf.reshape(RP.shape)
        valid = np.isfinite(zeta) & np.isfinite(rho)
        # interpolate the (smooth) shape; NaN -> 0 so fill marks invalid
        kw = dict(bounds_error=False, fill_value=np.nan)
        self._zeta_i = RegularGridInterpolator((rp_g, rpar_g), np.where(valid, zeta, np.nan), **kw)
        self._rho_i = RegularGridInterpolator((rp_g, rpar_g), np.where(valid, rho, np.nan), **kw)
        self._lam_i = RegularGridInterpolator((rp_g, rpar_g), lam_g, **kw)
        self._coverage = valid.mean()

    def __call__(self, v_los, r_perp, r_par):
        v_los = np.asarray(v_los, dtype=np.float64)
        rp = np.atleast_1d(np.asarray(r_perp, dtype=np.float64))
        rpar = np.atleast_1d(np.asarray(r_par, dtype=np.float64))
        mean, c2, g1, g2 = self._los_moments(rp, rpar)
        std = np.sqrt(c2)
        out_shape = np.broadcast_shapes(v_los.shape, mean.shape)

        pts = np.stack([np.broadcast_arrays(rp, rpar)[0].ravel(),
                        np.broadcast_arrays(rp, rpar)[1].ravel()], axis=-1)
        zeta = self._zeta_i(pts).reshape(mean.shape)
        rho = self._rho_i(pts).reshape(mean.shape)
        lam = self._lam_i(pts).reshape(mean.shape)
        good = np.isfinite(zeta) & np.isfinite(rho) & (c2 > 0)

        with np.errstate(all="ignore"):
            _, _, var0 = _gh_shape_skew_kurt(lam, zeta, rho)
            b = np.sqrt(c2 / var0)
            alpha = 1.0 / b
            beta = rho / b
            delta = zeta / np.sqrt(1.0 - rho ** 2) * b
            gam = np.sqrt(1.0 - rho ** 2) / b
            R1 = special.kv(lam + 1, zeta) / special.kv(lam, zeta)
            mu = mean - beta * (delta / gam) * R1
            p = _gh_density(v_los, alpha, beta, delta, lam, mu)
        good = good & np.isfinite(p) & (p >= 0)

        # fallback to analytic NIG where the GH table is invalid/out of range
        bad = np.broadcast_to(~good, out_shape)
        self._count_fallback(bad)
        if np.any(bad):
            p_nig = self._nig(v_los, rp, rpar)
            p = np.where(bad, p_nig, p)
        return np.broadcast_to(p, out_shape)

Scale-dependent

liulu.physics.velocity_pdf.scale_dependent_lambda 

scale_dependent_lambda(rp_grid, rpar_grid, lam_grid, fill=-0.5)

Build a callable lam(r_perp, r_par) interpolating a calibrated lambda(separation) table, for :class:TabulatedGHPDF. Outside the grid it returns fill (default -0.5 = NIG), so the model degrades gracefully to NIG on scales where no calibration exists.

Source code in liulu/physics/velocity_pdf.py 
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
def scale_dependent_lambda(rp_grid, rpar_grid, lam_grid, fill=-0.5):
    """Build a callable ``lam(r_perp, r_par)`` interpolating a calibrated
    lambda(separation) table, for :class:`TabulatedGHPDF`. Outside the grid it
    returns ``fill`` (default -0.5 = NIG), so the model degrades gracefully to
    NIG on scales where no calibration exists.
    """
    from scipy.interpolate import RegularGridInterpolator
    rgi = RegularGridInterpolator(
        (np.asarray(rp_grid, dtype=np.float64), np.asarray(rpar_grid, dtype=np.float64)),
        np.asarray(lam_grid, dtype=np.float64),
        bounds_error=False, fill_value=float(fill))

    def lam(r_perp, r_par):
        rp, rpar = np.broadcast_arrays(np.atleast_1d(np.asarray(r_perp, dtype=np.float64)),
                                       np.atleast_1d(np.asarray(r_par, dtype=np.float64)))
        out = rgi(np.stack([rp.ravel(), rpar.ravel()], axis=-1)).reshape(rp.shape)
        return np.where(np.isfinite(out), out, fill)

    return lam

LOS projection

liulu.physics.los_projection.project_to_los 

project_to_los(moments: PairwiseVelocityMoments, n: int, mode: str = 'c')

Project the n-th order r-t moment onto the line of sight.

Parameters:

Name Type Description Default
moments PairwiseVelocityMoments

Object providing get_moment(r, r_order, t_order, mode).

 required
n int

Order of the LOS moment.

 required
mode str

'c' for central moments, 'm' for moments about the origin.

 'c'

Returns:

Type Description
callable

los_moment(r_perp, r_par) -> ndarray returning the projected n-th LOS moment.
Source code in liulu/physics/los_projection.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
def project_to_los(moments: PairwiseVelocityMoments, n: int,
                   mode: str = "c"):
    """Project the n-th order r-t moment onto the line of sight.

    Parameters
    ----------
    moments : PairwiseVelocityMoments
        Object providing ``get_moment(r, r_order, t_order, mode)``.
    n : int
        Order of the LOS moment.
    mode : str
        ``'c'`` for central moments, ``'m'`` for moments about the origin.

    Returns
    -------
    callable
        ``los_moment(r_perp, r_par) -> ndarray`` returning the projected
        n-th LOS moment.
    """

    def los_moment(r_perp, r_par, *, _r=None, _mu=None):
        if _r is None:
            r, mu = separation_geometry(r_perp, r_par)
        else:
            r, mu = _r, _mu

        result = np.zeros_like(r)
        for k in range(n + 1):
            t_order = n - k
            coeff = comb(n, k, exact=True)
            m = moments.get_moment(r, r_order=k, t_order=t_order, mode=mode)
            result = result + coeff * mu ** k * (1.0 - mu ** 2) ** (t_order / 2.0) * m

        return result

    return los_moment

Cosmology

liulu.physics.cosmology.Cosmology

Bases: FlatLambdaCDM

Background cosmology based on :class:astropy.cosmology.FlatLambdaCDM.

Parameters:

Name Type Description Default
Om0 float

Present-day matter density parameter.

 required
H0 float

Hubble constant in km/s/Mpc.

 required
z float

Redshift of evaluation.

 required
gamma_growth float

Growth index gamma in f = Omega_m(z)**gamma (Linder 2005). Default 0.55 (GR/LCDM, ~0.2% accurate). Override for modified-gravity growth (e.g. ~0.42-0.46 for scalar-tensor / f(R)-like enhancement); note the true MG growth rate is generally scale-dependent, so a single gamma is itself an approximation.

 0.55
Source code in liulu/physics/cosmology.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
class Cosmology(FlatLambdaCDM):
    """Background cosmology based on :class:`astropy.cosmology.FlatLambdaCDM`.

    Parameters
    ----------
    Om0 : float
        Present-day matter density parameter.
    H0 : float
        Hubble constant in km/s/Mpc.
    z : float
        Redshift of evaluation.
    gamma_growth : float, optional
        Growth index gamma in ``f = Omega_m(z)**gamma`` (Linder 2005).
        Default 0.55 (GR/LCDM, ~0.2% accurate). Override for
        modified-gravity growth (e.g. ~0.42-0.46 for scalar-tensor /
        f(R)-like enhancement); note the true MG growth rate is generally
        scale-dependent, so a single gamma is itself an approximation.
    """

    def __init__(self, Om0: float, H0: float, z: float,
                 gamma_growth: float = 0.55, **kwargs):
        super().__init__(H0=H0, Om0=Om0, **kwargs)
        self.z = z
        self.gamma_growth = float(gamma_growth)

    @property
    def a(self) -> float:
        """Scale factor."""
        return 1.0 / (1.0 + self.z)

    def Hz(self) -> float:
        """Hubble parameter H(z) in km/s/Mpc."""
        return self.H(self.z).value

    def aH(self) -> float:
        """a * H(z) in km/s/(Mpc/h).

        For distances in Mpc/h:
            v[km/s] = a*H(z)[km/s/Mpc] * delta_r[Mpc/h] * (1/h)[Mpc/(Mpc/h)]

        So aH in km/s/(Mpc/h) = a * H(z) / h.
        """
        return self.a * self.Hz() / self.h

    def f(self) -> float:
        """Linear growth rate f = d ln D / d ln a.

        Uses the approximation f ~ Omega_m(z)**gamma_growth (Linder 2005);
        gamma_growth defaults to 0.55 (GR/LCDM, ~0.2% accurate) and can be
        overridden at construction for modified-gravity growth.
        """
        Omz = self.Om0 * (1 + self.z)**3 / self.efunc(self.z)**2
        return Omz**self.gamma_growth

    def D(self) -> float:
        """Linear growth factor D(z), normalised to D(0) = 1.

        Uses the Carroll, Press & Turner (1992) fitting formula
        (~1-2% accurate for LCDM; not valid for modified gravity).
        """
        def _growth_unnorm(Om, Ode):
            return (5.0 / 2.0) * Om / (
                Om**(4.0 / 7.0) - Ode
                + (1.0 + Om / 2.0) * (1.0 + Ode / 70.0)
            )

        Ez2 = self.efunc(self.z)**2
        Omz = self.Om0 * (1 + self.z)**3 / Ez2
        Odez = self.Ode0 / Ez2

        Dz = _growth_unnorm(Omz, Odez) / (1.0 + self.z)
        D0 = _growth_unnorm(self.Om0, self.Ode0)
        return Dz / D0

a property 

a: float

Scale factor.

Hz 

Hz() -> float

Hubble parameter H(z) in km/s/Mpc.

Source code in liulu/physics/cosmology.py 
40
41
42
def Hz(self) -> float:
    """Hubble parameter H(z) in km/s/Mpc."""
    return self.H(self.z).value

aH 

aH() -> float

a * H(z) in km/s/(Mpc/h).

For distances in Mpc/h: v[km/s] = a*H(z)[km/s/Mpc] * delta_r[Mpc/h] * (1/h)[Mpc/(Mpc/h)]

So aH in km/s/(Mpc/h) = a * H(z) / h.

Source code in liulu/physics/cosmology.py 
44
45
46
47
48
49
50
51
52
def aH(self) -> float:
    """a * H(z) in km/s/(Mpc/h).

    For distances in Mpc/h:
        v[km/s] = a*H(z)[km/s/Mpc] * delta_r[Mpc/h] * (1/h)[Mpc/(Mpc/h)]

    So aH in km/s/(Mpc/h) = a * H(z) / h.
    """
    return self.a * self.Hz() / self.h

f 

f() -> float

Linear growth rate f = d ln D / d ln a.

Uses the approximation f ~ Omega_m(z)**gamma_growth (Linder 2005); gamma_growth defaults to 0.55 (GR/LCDM, ~0.2% accurate) and can be overridden at construction for modified-gravity growth.

Source code in liulu/physics/cosmology.py 
54
55
56
57
58
59
60
61
62
def f(self) -> float:
    """Linear growth rate f = d ln D / d ln a.

    Uses the approximation f ~ Omega_m(z)**gamma_growth (Linder 2005);
    gamma_growth defaults to 0.55 (GR/LCDM, ~0.2% accurate) and can be
    overridden at construction for modified-gravity growth.
    """
    Omz = self.Om0 * (1 + self.z)**3 / self.efunc(self.z)**2
    return Omz**self.gamma_growth

D 

D() -> float

Linear growth factor D(z), normalised to D(0) = 1.

Uses the Carroll, Press & Turner (1992) fitting formula (~1-2% accurate for LCDM; not valid for modified gravity).

Source code in liulu/physics/cosmology.py 
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
def D(self) -> float:
    """Linear growth factor D(z), normalised to D(0) = 1.

    Uses the Carroll, Press & Turner (1992) fitting formula
    (~1-2% accurate for LCDM; not valid for modified gravity).
    """
    def _growth_unnorm(Om, Ode):
        return (5.0 / 2.0) * Om / (
            Om**(4.0 / 7.0) - Ode
            + (1.0 + Om / 2.0) * (1.0 + Ode / 70.0)
        )

    Ez2 = self.efunc(self.z)**2
    Omz = self.Om0 * (1 + self.z)**3 / Ez2
    Odez = self.Ode0 / Ez2

    Dz = _growth_unnorm(Omz, Odez) / (1.0 + self.z)
    D0 = _growth_unnorm(self.Om0, self.Ode0)
    return Dz / D0