`progenax.profiles`

progenax.profiles¶

Auto-generated by scripts/build_api_reference.py from progenax source. Re-run after public-API changes; the script is idempotent.

Module path: progenax/profiles/

Public symbols: 15

Contents¶

• PlummerProfile

• KingProfile

• solve_king_profile

• MichieProfile

• solve_michie_profile

• LIMEPYProfile

• solve_limepy_profile

• solve_multimass_limepy

• find_alpha_for_masses

• EFFProfile

• UniformSphereProfile

• ProfileName

• make_profile

• sample_density_profile

• compute_profile_potential

profiles.PlummerProfile¶

class

📇 model card · ∇ gradient-verified — 2 audit cases

PlummerProfile(r_h: float = 1.0)

Plummer (1911) spherical density profile.

ρ(r) = (3M / 4πa³) × (1 + r²/a²)^(-5/2)

Attributes

References. Plummer (1911) MNRAS 71, 460

Examples. >>> profile = PlummerProfile(r_h=1.0) # 1 pc

masses = jnp.ones(100) key = jax.random.PRNGKey(42) positions = profile.sample_positions(masses, key)

Source: src/progenax/profiles/plummer.py#L16

profiles.KingProfile¶

class

📇 model card · ∇ gradient-verified — 4 audit cases

KingProfile(W0: ArrayLike, r_c: ArrayLike, r_t: ArrayLike, xi_grid: Float[Array, 'n_points'], psi_grid: Float[Array, 'n_points'], n_grid: int = 1000)

King (1966) spherical density profile.

Implements SpatialProfile protocol for IC assembly.

The CDF is precomputed at initialization for efficient sampling.

Attributes

References. King (1966), AJ, 71, 64

Examples. # Recommended: Use from_W0_rc for self-consistent model

profile = KingProfile.from_W0_rc(W0=7.0, r_c=1.0)

Or manually with pre-computed ODE solution¶

xi_grid, psi_grid, _ = solve_king_profile(W0=7.0) profile = KingProfile(W0=7.0, r_c=1.0, r_t=10.0, ... xi_grid=xi_grid, psi_grid=psi_grid) masses = jnp.ones(100) key = jax.random.PRNGKey(42) positions = profile.sample_positions(masses, key)

Source: src/progenax/profiles/king.py#L330

profiles.solve_king_profile¶

function

solve_king_profile(W0: ArrayLike, xi_max: float = 300.0, n_points: int = 2000)

Solve King’s Poisson equation numerically using diffrax.

Integrates from xi=0 (center) outward until psi(xi) -> 0 (tidal radius).

Boundary conditions (King 1966, Eq. 10): psi(0) = W0 (central potential) d psi/d xi|_0 = 0 (symmetry at center)

Args

Returns: 3-tuple (xi_grid, psi_clamped, psi_raw): - xi_grid: dimensionless radii (solution.ts). - psi_clamped: max(psi, 0), truncated at the tidal radius. This is the physical potential used for density / CDF / mu / virial. - psi_raw: the UNCLAMPED ODE solution solution.ys[:, 0] (goes negative beyond the zero-crossing). Feed psi_raw to _find_tidal_radius so the zero-crossing linear interpolation carries d(xi_t)/dW0 from the diffrax solve -- the clamp would otherwise set psi=0 at the crossing node, killing the gradient (audit Task 1.2b, RESOLVED). The forward value of xi_t is identical either way; only the gradient differs.

References. King (1966), AJ, 71, 64 Binney & Tremaine (2008), Section 4.3.2

Note. JIT-compatible when n_points (and xi_max) are static: they set the linspace size and are closed over, so jax.jit(solve_king_profile)(W0) traces fine (W0 may be a tracer). Uses Tsit5 (Runge-Kutta 5th order) from diffrax for robustness.

Source: src/progenax/profiles/king.py#L193

profiles.MichieProfile¶

class

📇 model card · ∇ gradient-verified — 4 audit cases

MichieProfile(W0, r_c, r_a, r_t, xi_grid, psi_grid, n_grid: int = 1000)

Michie-King anisotropic spherical density profile (SpatialProfile).

The self-consistent density of the Michie-King model (Michie 1963 anisotropy + King 1966 cutoff). More centrally-radial and more extended than the isotropic King model; r_a -> infinity recovers King. Construct with from_W0_rc.

Attributes

References. Michie (1963), MNRAS 125, 127; King (1966), AJ 71, 64.

Source: src/progenax/profiles/michie.py#L152

profiles.solve_michie_profile¶

function

📇 model card

solve_michie_profile(W0: ArrayLike, ra_hat: ArrayLike, xi_max: float = 800.0, n_points: int = 3000)

Solve the Michie-King Poisson equation from the centre outward to psi -> 0.

Args

Returns: 3-tuple (xi_grid, psi_clamped, psi_raw). psi_clamped is psi >= 0, truncated at the tidal radius (the physical potential used for density / CDF / mu). psi_raw is the UNCLAMPED ODE solution (negative past the crossing) -- feed it to _find_tidal_radius so d(xi_t)/dW0 flows (audit Task 1.2b; see solve_king_profile). The forward value of xi_t is the same either way; only the gradient differs.

References. Michie (1963), MNRAS 125, 127 (Eq. 5.8); King (1966), AJ 71, 64.

Source: src/progenax/profiles/michie.py#L82

profiles.LIMEPYProfile¶

class

LIMEPYProfile(W0, g, r_c, r_t, xi_grid, psi_grid, r_a=None, n_grid: int = 1000)

General-g LIMEPY (lowered-isothermal) spherical density profile.

Implements the SpatialProfile protocol. Generalizes KingProfile with a continuous truncation parameter g (g=0 Woolley, g=1 King, g=2 Wilson) and an optional Michie/Osipkov-Merritt radial anisotropy radius r_a. Corners: g=1 isotropic == KingProfile; g=1 anisotropic == MichieProfile.

The CDF is precomputed at construction for inverse-transform position sampling.

Attributes

Source: src/progenax/profiles/limepy.py#L321

profiles.solve_limepy_profile¶

function

solve_limepy_profile(W0: ArrayLike, g: ArrayLike, ra_hat: ArrayLike | None = None, xi_max: float = 300.0, n_points: int = 2000) -> Tuple[Float[Array, 'n_points'], Float[Array, 'n_points'], Float[Array, 'n_points']]

Solve the general-g (optionally anisotropic) LIMEPY Poisson equation (diffrax).

Integrates W(xi) from the centre (W=W0, dW/dxi=0) outward to the truncation radius where W -> 0. Generalizes solve_king_profile with the continuous truncation parameter g (g=0 Woolley, g=1 King, g=2 Wilson) and an optional Michie/OM radial anisotropy radius. At g=1: isotropic reproduces King, anisotropic reproduces the Michie-King model.

JIT/grad-safe in (W0, g, ra_hat): n_points and xi_max are static; W0, g, ra_hat may be tracers and the ODE -> W(xi) path carries their gradients.

Args

Returns: 3-tuple (xi_grid, psi_clamped, psi_raw) (mirrors solve_king_profile): - xi_grid: dimensionless radii. - psi_clamped: max(psi, 0) -- the physical potential W(xi) >= 0 used for density / CDF / mu / virial. - psi_raw: the UNCLAMPED ODE solution (negative past the zero-crossing). Feed psi_raw to _find_tidal_radius so the crossing interpolation carries d(xi_t)/dW0 -- the clamp would zero psi at the crossing node and kill that gradient (audit Task 1.2b pattern). The forward xi_t value is identical either way; only the gradient differs.

Source: src/progenax/profiles/limepy.py#L212

profiles.solve_multimass_limepy¶

function

solve_multimass_limepy(alpha_j: Float[Array, 'n_comp'], m_j: Float[Array, 'n_comp'], W0: ArrayLike, g: ArrayLike, delta: ArrayLike, xi_max: float = 300.0, n_points: int = 2000, ra_hat: float | None = None, eta: float = 0.0, aniso_method: str = 'table') -> Tuple[Float[Array, 'n_points'], Float[Array, 'n_points'], Float[Array, 'n_points'], Float[Array, 'n_comp n_points']]

Mass-segregation convenience over solve_multicomponent_limepy (Engine A).

The Gieles & Zocchi (2015) equipartition parametrization: per-component rescaling rescale_j = mu_j^(2 delta), mu_j = m_j / bar_m, bar_m = sum_j m_j alpha_j (central density-weighted mean mass, Eq. 26), and per-component anisotropy radius hat_r_{a,j} = ra_hat * mu_j^eta (eta=0 = mass-independent, the paper default). At delta=0 every rescaling is 1 -- identical to solve_limepy_profile. ra_hat=None is the (fast) isotropic case. aniso_method (“table” default, “quadrature” oracle) is passed through to solve_multicomponent_limepy; ignored when ra_hat is None.

JIT/grad-safe in (alpha_j, m_j, W0, g, delta, ra_hat, eta); n_points, xi_max, aniso_method static.

Returns (xi_grid, psi_grid, psi_raw, rho_j_grid) as solve_multicomponent_limepy.

Source: src/progenax/profiles/limepy_multimass.py#L359

profiles.find_alpha_for_masses¶

function

find_alpha_for_masses(m_j: Float[Array, 'n_comp'], M_j: Float[Array, 'n_comp'], W0: ArrayLike, g: ArrayLike, delta: ArrayLike, n_iter: int = 30, xi_max: float = 300.0, n_points: int = 2000, ra_hat=None, eta: float = 0.0, aniso_method: str = 'table', tol: float = 1e-06) -> Tuple[Float[Array, 'n_comp'], Float[Array, '']]

Find the central density fractions alpha_j that reproduce target masses M_j (Layer B).

Fixed-point iteration on the realized mass fractions, with the stabilized update of Gieles & Zocchi (2015, Section 4.1) -- NOT Gunn & Griffin’s linear M_j/M_j’ (which diverges for wide mass functions):

alpha_j <- alpha_j sqrt(f_j / f_j'),   renormalize sum_j alpha_j = 1,

f_j = M_j / sum M (target), f_j’ = realized fraction. Starts from alpha_j = f_j.

Solved by a hand-rolled jax.custom_vjp: the forward is an adaptive jax.lax.while_loop that iterates the sqrt-update until the residual max_j |f_j’ - f_j| < tol (or the n_iter safety cap), and the backward is the EXACT fixed-point gradient via a reverse-mode implicit VJP of the sqrt-map residual R(alpha, theta) = alpha - sqrt-map(alpha) (n x n Jacobian by vmapped vjp, lstsq solve, -vjp_theta). This is flat-in-n_iter and ~3x faster per value_and_grad than the old unrolled lax.scan, with gradients matching central finite differences to <1e-5. Two solvers are dispatched on ra_hat is None: the isotropic solver keeps (ra_hat, eta) out of the differentiated set so it can take the fast isotropic density path; the anisotropic solver also differentiates (ra_hat, eta).

The iteration deliberately uses the SAME aniso_method as the final solve (default “table”) so the converged alpha_j are self-consistent with the model actually built; the residual remains a reported diagnostic. Pass aniso_method=“quadrature” for the exact oracle path.

Args

Returns: (alpha_j, residual): converged central density fractions (sum to 1, positive) and the final fractional residual max_j |f_j’ - f_j| (reported, never branched on).

Source: src/progenax/profiles/limepy_multimass.py#L648

profiles.EFFProfile¶

class

📇 model card · ∇ gradient-verified — 4 audit cases

EFFProfile(a: float = 1.0, gamma: float = 3.0, r_t: float = 10.0, n_grid: int = 1000)

EFF (Elson-Fall-Freeman 1987) truncated power-law profile.

3-D volume density: rho(r) = rho_0 * (1 + r^2/a^2)^{-gamma/2} for r <= r_t = 0 for r > r_t

Provenance note: Elson, Fall & Freeman (1987) Eq. 1 defines this functional form as the projected surface brightness mu(r), with gamma the surface (projected) slope (their median ~2.6). Here we adopt the same form as the 3-D volume density -- the standard N-body / IC-code convention -- so this gamma is a 3-D density slope, offset by ~1 from EFF87’s surface slope (Abel projection of r^-gamma gives a surface slope gamma-1). At gamma=5 the form reduces exactly to Plummer. No closed-form DF; velocities are assigned via Eddington inversion in kinematics.EFFVelocityDF.

The CDF is precomputed at initialization for efficient sampling.

Attributes

References. Elson, Fall & Freeman (1987), ApJ, 323, 54 (Eq. 1 = surface brightness, used here as the 3-D volume density; see docs bibliography note).

Source: src/progenax/profiles/eff.py#L20

profiles.UniformSphereProfile¶

class

UniformSphereProfile(R: float = 1.0)

Uniform density sphere profile.

ρ(r) = ρ₀ for r ≤ R, 0 otherwise

This is the CW04 ‘3D0’ distribution with Q ≈ 0.79.

Attributes

References. Cartwright & Whitworth (2004) MNRAS 348, 589 - Table 1, ‘3D0’

Examples. >>> profile = UniformSphereProfile(R=1.0) # 1 pc

masses = jnp.ones(100) key = jax.random.PRNGKey(42) positions = profile.sample_positions(masses, key)

Source: src/progenax/profiles/uniform.py#L15

profiles.ProfileName¶

function

ProfileName(*args, **kwargs)

(no docstring)

profiles.make_profile¶

function

make_profile(name: Literal['plummer', 'king', 'eff'], R_half: float, **kwargs) -> Union[progenax.profiles.plummer.PlummerProfile, progenax.profiles.king.KingProfile, progenax.profiles.eff.EFFProfile]

Factory function for creating profile instances.

This is the primary entry point for creating density profiles. It provides a uniform interface where R_half is the external scale parameter, with profile-specific shape parameters passed via kwargs.

Args

Returns: Profile instance (PlummerProfile, KingProfile, or EFFProfile)

Raises: ValueError: If profile name is not recognized

Examples. >>> # Plummer profile with 1 pc half-mass radius

profile = make_profile(“plummer”, R_half=1.0)

King profile with W0=7 (globular cluster typical)¶

profile = make_profile(“king”, R_half=1.0, W0=7.0)

EFF profile with gamma=3 (young cluster typical)¶

profile = make_profile(“eff”, R_half=1.0, gamma=3.0, r_t=15.0)

Notes. For King and EFF profiles, the mapping from R_half to internal parameters is simplified: R_half is used directly as r_c (King) or a (EFF). For precise half-mass radius control, the user should compute the appropriate r_c or a value externally.

Source: src/progenax/profiles/api.py#L41

profiles.sample_density_profile¶

function

sample_density_profile(key: PRNGKeyArray, N_stars: int, profile: Literal['plummer', 'king', 'eff'], R_half: float, **kwargs) -> Float[Array, 'N 3']

Sample N_stars positions from the chosen density profile.

This function creates a profile instance and samples positions from its density distribution using inverse CDF sampling.

Args

Returns: positions: Array of shape (N_stars, 3) in length units [pc] Positions are sampled from the density profile’s radial distribution with isotropic angular distribution.

Examples. >>> import jax

key = jax.random.PRNGKey(42)

Sample 1000 positions from Plummer profile¶

positions = sample_density_profile(key, 1000, “plummer”, R_half=1.0) positions.shape (1000, 3)

Sample from King profile with specific concentration¶

positions = sample_density_profile(key, 1000, “king”, R_half=1.0, W0=7.0)

Notes. - All profiles sample radii via inverse CDF for efficiency

• Angular distribution is isotropic (uniform on sphere)

• The masses argument required by profile.sample_positions() is filled with ones internally (mass values don’t affect spatial sampling)

Source: src/progenax/profiles/api.py#L122

profiles.compute_profile_potential¶

function

compute_profile_potential(positions: Float[Array, 'N 3'], profile: Literal['plummer', 'king', 'eff'], M_total: float, R_half: float, G: float, **kwargs) -> Float[Array, 'N']

Compute analytic gravitational potential at given positions.

This function computes Phi(r) for the specified density profile. The potential is normalized to the total mass M_total and gravitational constant G.

Args

Returns: phi: SPECIFIC gravitational potential (per unit mass) at each position, shape (N,). All values are negative (bound system). Units: [length units]^2 / [time units]^2 (e.g. pc^2/Myr^2 in STELLAR units) — consistent with specific kinetic energy 0.5v^2, so E = 0.5v^2 + phi is a specific orbital energy.

Examples. >>> from jaxstro.units import STELLAR

import jax.numpy as jnp

positions = jnp.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0]]) phi = compute_profile_potential( ... positions, “plummer”, M_total=1000.0, R_half=1.0, G=STELLAR.G ... )

Notes. Plummer potential (analytic): Phi(r) = -G * M_total / sqrt(r^2 + a^2) where a = R_half * sqrt(2^(2/3) - 1) is the scale radius.

King potential (exact relative potential): Phi(r) = -sigma^2 * psi(r), with psi the dimensionless King ODE potential (psi(r_t) = 0) and sigma^2 = G M / (9 r_c mu(W0)) the self-consistent velocity scale (matches KingVelocityDF).

EFF potential (true spherical potential): Phi(r) = -G [ M(<r)/r + 4 pi int_r^rt rho s ds ], using the exact enclosed mass from the profile’s density grid (interior monopole + outer-shell term).

The potential is computed per-particle and vectorized for efficiency.

Source: src/progenax/profiles/api.py#L176