mass_assignment

Mass assignment and velocity moments on regular 3D grids, implemented in C++ with OpenMP.

Installation

pip install mass-assignment

Inputs

Inputs are NumPy arrays prepared by the caller:

• pos: particle positions, (N,3) or (3,N)
• vel: particle velocities, (N,3) or (3,N)
• mass: optional particle weights, (N,); unit weights when omitted

Positions and velocities may be float32 or float64. Helpers for applying the same coordinate, box, and boundary conventions are available from mass_assignment.inputs.

Usage

import numpy as np
import mass_assignment as ma

rng = np.random.default_rng(42)
pos = rng.uniform(0.0, 100.0, size=(100_000, 3))
vel = rng.normal(size=(100_000, 3))
mass = rng.uniform(0.5, 1.5, size=100_000)

grid = ma.Grid(nmesh=128, lbox=100.0, method="TSC", nthreads=0)

rho = grid.dens(pos, mass=mass)
vx, vy, vz = grid.velc(pos, vel, mass=mass)
m2 = grid.moment2(pos, vel, mass=mass)
m3 = grid.moment3(pos, vel, mass=mass)
m4 = grid.moment4(pos, vel, mass=mass)

dens returns mass per cell, not density contrast. For equal-volume cells:

delta = rho / rho.mean() - 1.0

The module-level functions dens, scalar, velc, moment2, moment3, and moment4 accept the same settings directly. The names sigma, skewness, and kurtosis are aliases for moments of order 2, 3, and 4; the latter two are not standardized skewness or excess kurtosis.

Moment functions return one array per independent Cartesian component. Use np.asarray(result) when a stacked array is wanted. The component labels are available through moment_components(order).

Grid settings

• method: "NGP", "CIC", "TSC", or "PCS"; integer orders 1–4 are also accepted.
• nmesh: a positive scalar or three axis sizes.
• lbox: a positive scalar or three box lengths, in the units of pos.
• bc: "periodic", "open", or one label per axis. Open boundaries discard the part of an assignment kernel lying outside the box.
• dtype: output precision, np.float32 by default or np.float64.
• nthreads: 0 uses the OpenMP default; a positive integer sets the count.

Mesh utilities

values = ma.mesh_to_ptcl(pos, rho, lbox=100.0, method="TSC")
gx, gy, gz = ma.mesh_diff(rho, lbox=100.0, order=4)

mesh_to_ptcl samples a float32 scalar mesh at particle positions. mesh_diff supports finite-difference orders 2, 4, 6, and 8. An open axis needs at least order + 1 cells.

Shot noise

Grid.shotnoise returns the leading Poisson term V * sum(w**2) / sum(w)**2, or V/N for equal weights:

shotnoise = grid.shotnoise(pos, mass=mass)

An explicit lbox= may be supplied when the physical volume differs from the coordinate normalization used for deposition.

Moment component order

m2: xx, xy, xz, yy, yz, zz
m3: xxx, xxy, xxz, xyy, xyz, xzz, yyy, yyz, yzz, zzz
m4: xxxx, xxxy, xxxz, xxyy, xxyz, xxzz,
    xyyy, xyyz, xyzz, xzzz,
    yyyy, yyyz, yyzz, yzzz, zzzz

License

MIT; see LICENSE.