pisces.physics.virialization.eddington.sample_eddington_velocities#

pisces.physics.virialization.eddington.sample_eddington_velocities(relative_energy_array: unyt_array, distribution_function_array: unyt_array, particle_relative_potentials: unyt_array, **kwargs) → unyt_array[source]#

Sample 3D Cartesian velocities for particles based on an isotropic Eddington distribution function.

This function uses rejection sampling to draw particle speeds from a spherically symmetric, ergodic distribution function \(f(\mathcal{E})\), and then assigns a random direction to each speed to produce 3D velocity vectors.

Parameters:

  • relative_energy_array (unyt_array) – 1D array containing the tabulated relative energies \(\mathcal{E}\) at which the distribution function \(f(\mathcal{E})\) is evaluated. Units must be compatible with \(\mathrm{km}^2\,\mathrm{s}^{-2}\). The length of this array must match the length of distribution_function_array.

  • distribution_function_array (unyt_array) – 1D array of the distribution function values \(f(\mathcal{E})\), corresponding to each entry in relative_energy_array. These values describe the phase-space density of particles as a function of energy in a spherically symmetric, isotropic system. Units must be \(\mathrm{M}_\odot\,\mathrm{s}^3\,\mathrm{km}^{-3}\,\mathrm{kpc}^{-3}\).

  • particle_relative_potentials (unyt_array) – 1D array of relative gravitational potentials \(\Psi(r)\) for each particle position, used to determine the maximum allowed velocity (escape velocity) at that location. Each value should be convertible to \(\mathrm{km}^2\,\mathrm{s}^{-2}\). The escape velocity for a particle is computed as \(v_\mathrm{esc} = \sqrt{2\Psi(r)}\).

  • **kwargs – Additional keyword arguments passed to the Cython sampling routine, such as show_progress to control progress bar visibility.

Returns:

output_velocities – A (3, N) array of sampled 3D Cartesian velocities with units of km/s. The velocities are sampled from the input Eddington distribution and distributed isotropically in direction.

Return type:

unyt_array

Notes

This function samples velocities for particles in a spherically symmetric, isotropic system governed by an ergodic distribution function \(f(\mathcal{E})\). For each particle, the escape velocity is computed from the relative potential using the expression \(v_\mathrm{esc}(r) = \sqrt{2\Psi(r)}\), where the relative potential is defined as \(\Psi(r) = -[\Phi(r) - \Phi_0]\). This gives the maximum speed a bound particle can have at that location.

The distribution function \(f(\mathcal{E})\) is provided as a one-dimensional array evaluated on a corresponding grid of relative energies \(\mathcal{E}\), and is interpolated internally using a cubic spline. Particle speeds are sampled via rejection sampling on the function \(f(\mathcal{E}) \cdot v^2\), which reflects the correct velocity-space weighting in spherical symmetry.

For each particle, trial velocities \(v\) are drawn uniformly in \([0, v_\mathrm{esc}]\) and converted to relative energy via \(\mathcal{E} = \Psi - \frac{1}{2}v^2\). A trial is accepted if the product \(f(\mathcal{E}) \cdot v^2\) is less than a uniform random deviate times a precomputed bounding value \(\mathrm{likelihood\_max}[i]\). This ensures unbiased sampling from the desired distribution.

Once a speed is accepted, a random isotropic direction is assigned by sampling spherical angles \(\theta\) and \(\phi\), and converting the result into Cartesian velocity components. The final output is a three-component array of velocities for each particle, sampled from \(f(\mathcal{E})\) and distributed isotropically in direction.

The core rejection sampling loop is implemented in Cython for performance and lives in the module pisces/particles/sampling/_eddington_sampling.pyx. All input arrays must be one-dimensional, consistently ordered, and convertible to the expected physical units.