Skip to content

halocat.velocity

velocity

Pairwise velocity moments via the pairvel (PairVel.jl) wrapper.

PAIRVEL_PYTHON module-attribute 

PAIRVEL_PYTHON = '/cosma8/data/dp203/dc-ruan1/pairvel/python'

MOMENT_KEYS module-attribute 

MOMENT_KEYS = ('npairs', 'm10', 'c20', 'c02', 'c12', 'c30', 'c40', 'c04', 'c22')

measure_velocity_moments 

measure_velocity_moments(data: dict, box_size: float, mass_bins: list[tuple[float, float]], r_edges: ndarray, mass_key: str = 'Mtot') -> dict

Compute pairwise velocity central moments per mass-bin pair.

Returns a dict keyed by 'M{i}_M{j}' (i <= j) with the standard PairVel fields plus the bin centres and mass-bin metadata.

Source code in halocat/velocity.py 
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
def measure_velocity_moments(data: dict, box_size: float,
                             mass_bins: list[tuple[float, float]],
                             r_edges: np.ndarray,
                             mass_key: str = "Mtot") -> dict:
    """Compute pairwise velocity central moments per mass-bin pair.

    Returns a dict keyed by 'M{i}_M{j}' (i <= j) with the standard PairVel
    fields plus the bin centres and mass-bin metadata.
    """
    _ensure_pairvel_on_path()
    from pairvel import compute_moments

    r_edges = np.asarray(r_edges, dtype=np.float64)
    r_centres = 0.5 * (r_edges[:-1] + r_edges[1:])

    selections = []
    for lo, hi in mass_bins:
        sel = _select_mass_bin(data, lo, hi, mass_key=mass_key)
        pos = np.array([data["x"][sel], data["y"][sel], data["z"][sel]])
        vel = np.array([data["vx"][sel], data["vy"][sel], data["vz"][sel]])
        selections.append((pos, vel))

    out: dict = {}
    nbins = len(mass_bins)
    for i in range(nbins):
        for j in range(i, nbins):
            posL, velL = selections[i]
            posR, velR = selections[j]
            n1 = posL.shape[1]
            n2 = posR.shape[1]
            entry: dict = {
                "r": r_centres,
                "log10M1": np.array(mass_bins[i], dtype=np.float64),
                "log10M2": np.array(mass_bins[j], dtype=np.float64),
                "n1": int(n1),
                "n2": int(n2),
            }
            if n1 == 0 or n2 == 0:
                for k in MOMENT_KEYS:
                    entry[k] = np.full(len(r_centres),
                                       0.0 if k == "npairs" else np.nan)
            else:
                if i == j:
                    res = compute_moments(posL, velL, r_edges, box_size)
                else:
                    res = compute_moments(
                        posL, velL, r_edges, box_size,
                        positions_right=posR, velocities_right=velR,
                    )
                for k in MOMENT_KEYS:
                    entry[k] = np.asarray(res[k])
            out[f"M{i}_M{j}"] = entry

    out["_meta"] = {
        "r_edges": r_edges,
        "box_size": float(box_size),
        "mass_bins": np.array(mass_bins, dtype=np.float64),
    }
    return out

write_velocity_moments 

write_velocity_moments(vel_dict: dict, outpath: str, attrs: dict | None = None) -> None
Source code in halocat/velocity.py 
 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
def write_velocity_moments(vel_dict: dict, outpath: str,
                           attrs: dict | None = None) -> None:
    os.makedirs(os.path.dirname(outpath), exist_ok=True)
    meta = vel_dict.get("_meta", {})
    with h5py.File(outpath, "w") as f:
        if "r_edges" in meta:
            f.create_dataset("r_edges", data=meta["r_edges"])
        if "mass_bins" in meta:
            f.create_dataset("mass_bins", data=meta["mass_bins"])
        f.attrs["box_size"] = meta.get("box_size", 0.0)
        f.attrs["moment_keys"] = list(MOMENT_KEYS)
        if attrs:
            for k, v in attrs.items():
                f.attrs[k] = v
        for key, val in vel_dict.items():
            if key == "_meta":
                continue
            g = f.create_group(key)
            g.create_dataset("r", data=val["r"])
            for k in MOMENT_KEYS:
                g.create_dataset(k, data=val[k])
            g.attrs["log10M1"] = val["log10M1"]
            g.attrs["log10M2"] = val["log10M2"]
            g.attrs["n1"] = val["n1"]
            g.attrs["n2"] = val["n2"]