Annotated: electrosim.simulation.physics#

This page provides function-by-function, line-by-line explanations and cross-references.

minimum_image_displacement#

from electrosim.config import (
    K_COULOMB,
    COLOR_POSITIVE,
    COLOR_NEGATIVE,
    COLOR_NEUTRAL,
    NEUTRAL_CHARGE_EPS,
)

def minimum_image_displacement(p_i: np.ndarray, p_j: np.ndarray, world_size_m: np.ndarray) -> np.ndarray:
    """Compute displacement from i to j using the minimum distance convention in a 2D torus."""
    delta = p_j - p_i
    for axis in (0, 1):
        L = world_size_m[axis]
        if delta[axis] > 0.5 * L:
            delta[axis] -= L
        elif delta[axis] < -0.5 * L:
            delta[axis] += L
    return delta

delta = p_j - p_i: raw displacement before periodicity. We adopt j − i to point from i to j (consistent throughout code).
Loop over axes applies minimum-image: if displacement exceeds half-box, wrap by one box length the other way. See math/coulomb_periodic.
Returns the 2D vector with the shortest equivalent displacement on the torus.

electric_force_pair#

def electric_force_pair(
    p_i: np.ndarray,
    q_i: float,
    r_i: float,
    p_j: np.ndarray,
    q_j: float,
    r_j: float,
    world_size_m: np.ndarray,
    softening_fraction: float,
) -> np.ndarray:
    r_vec = minimum_image_displacement(p_j, p_i, world_size_m)
    r2 = float(np.dot(r_vec, r_vec))
    contact = r_i + r_j
    epsilon = softening_fraction * contact
    den = (r2 + epsilon * epsilon) ** 1.5
    if den == 0.0:
        return np.zeros(2, dtype=float)
    coef = K_COULOMB * q_i * q_j / den
    return coef * r_vec

Minimum-image uses (p_j, p_i) to point from i towards j (the force on i is along +r_vec if sign is positive). This matches the vector form of Coulomb force.
Softening radius epsilon based on contact distance (see math/plummer_softening).
Denominator implements \((r^2 + \epsilon^2)^{3/2}\) for the vector form \(\propto \mathbf{r}/r^3\).
Returns a 2D force vector.

compute_accelerations (Numba path shape and algorithm)#

At a high level, we pack particle arrays (structure-of-arrays) and call Numba kernels.

Key choices:

Skip fixed, massless, and neutral particles early.
Minimum-image per-axis inside tight loops.
Contact-based softening.
Multiply by inverse mass to produce accelerations.

See kernels: _compute_accelerations_numba_serial, _compute_accelerations_numba_parallel.

total_potential_energy#

def total_potential_energy(particles: List["Particle"], world_size_m: np.ndarray) -> float:
    E = 0.0
    N = len(particles)
    for i in range(N):
        pi = particles[i]
        for j in range(i + 1, N):
            pj = particles[j]
            r_vec = minimum_image_displacement(pi.pos_m, pj.pos_m, world_size_m)
            r = float(np.hypot(r_vec[0], r_vec[1]))
            r_eff = max(r, 1e-6)
            if r_eff == 0.0:
                continue
            E += K_COULOMB * pi.charge_c * pj.charge_c / r_eff
    return E

Pairwise sum over i<j avoids double counting.
Distance uses minimum-image.
r_eff = max(r, 1e-6) guards singularities; note the modeling choice: no softening applied here. See math/potential_energy_modeling.

rk4_integrate (stages and state policy)#

Packs current positions/velocities.
Defines helpers to temporarily set state and compute accelerations.
Computes k1..k4 for both position and velocity components.
Restores original state before committing final updates.
Writes back only for non-fixed particles.

Mathematical details in math/rk4_derivation.

resolve_collisions (merge and elastic phases)#

First pass: merge opposite-charge overlaps.
- New mass m1+m2, charge q1+q2, radius sqrt(r1^2+r2^2).
- Momentum conservation if not fixed; fixed result if any fixed.
- Update color and merge histories; wrap position; delete j; reindex ids.
Second pass: elastic collisions for remaining overlapping pairs.
- Compute normal, penetration correction along normal; handle fixed/infinite mass.
- If separating, skip; else impulse with restitution e=1.

See derivation in math/elastic_collision_impulse.