diff --git a/src/test_particle_sim_1d/collisions.py b/src/test_particle_sim_1d/collisions.py
new file mode 100644
index 0000000..32e5c8e
--- /dev/null
+++ b/src/test_particle_sim_1d/collisions.py
@@ -0,0 +1,191 @@
+"""
+collisions.py
+
+Implements Monte Carlo Collision (MCC) logic for a 1D3V particle simulation.
+
+This module handles the probabilistic nature of particle-neutral collisions.
+It calculates collision probabilities based on the particle's relative velocity
+to the background gas and updates particle velocities using elastic scattering
+dynamics.
+
+The collision probability is determined by:
+    P_coll = 1 - exp(-n_g * sigma(g) * g * dt)
+
+Functions
+---------
+apply_elastic_collisions : Main function to check for and apply collisions.
+scatter_isotropic_3d     : Updates velocities of colliding particles (Elastic).
+
+"""
+
+from __future__ import annotations
+
+from collections.abc import Callable
+
+import numpy as np
+
+from . import particles
+
+
+def apply_elastic_collisions(
+    species: particles.Species,
+    neutral_density: float,
+    neutral_temp_K: float,
+    neutral_mass: float,
+    cross_section_func: Callable[[np.ndarray], np.ndarray],
+    dt: float,
+    seed: int | None = None,
+) -> int:
+    """
+    Check for collisions between test particles and a background gas, then
+    update velocities using elastic scattering.
+
+    Parameters
+    ----------
+    species : Species
+        The particle species object containing position and velocity arrays.
+    neutral_density : float
+        Number density of the background gas [m^-3].
+    neutral_temp_K : float
+        Temperature of the background gas in Kelvin [K].
+    neutral_mass : float
+        Mass of a single background gas particle [kg].
+    cross_section_func : Callable
+        A function that takes relative velocity array `g` [m/s] and returns
+        cross-section array `sigma` [m^2].
+    dt : float
+        Time step [s].
+    seed : int, optional
+        Random seed for reproducibility.
+
+    Returns
+    -------
+    int
+        The total number of collisions that occurred in this step.
+    """
+    rng = np.random.default_rng(seed)
+
+    # 1. Access active particle data
+    # We only care about the first N particles that are alive
+    N = species.N
+    if N == 0:
+        return 0
+
+    vx = species.vx[:N]
+    vy = species.vy[:N]
+    vz = species.vz[:N]
+
+    # 2. Calculate EXACT relative velocity 'g'
+
+    # Sample candidate neutrals for all active particles
+    v_neutral_all = particles.sample_maxwellian(
+        n=N,
+        mass=neutral_mass,
+        temperature={"K": neutral_temp_K},
+        mean_velocity=0.0,
+        seed=seed,
+    )
+
+    # Calculate vector relative velocity for every particle
+    v_rel_x = vx - v_neutral_all[:, 0]
+    v_rel_y = vy - v_neutral_all[:, 1]
+    v_rel_z = vz - v_neutral_all[:, 2]
+
+    # Compute magnitude g
+    g = np.sqrt(v_rel_x**2 + v_rel_y**2 + v_rel_z**2)
+
+    # 3. Calculate Cross Section sigma(g)
+    sigma = cross_section_func(g)
+
+    # 4. Calculate Collision Probability (Eq. 2 in proposal)
+    # P = 1 - exp(-n * sigma * g * dt)
+    exponent = -neutral_density * sigma * g * dt
+    P_coll = 1.0 - np.exp(exponent)
+
+    # 5. Determine which particles collide (Vectorized Monte Carlo)
+    # Generate random numbers [0, 1) for every particle
+    R = rng.random(N)
+    collision_mask = P_coll > R
+    n_collisions = np.count_nonzero(collision_mask)
+
+    if n_collisions == 0:
+        return 0
+
+    # 6. Resolve Collisions (Update Velocities)
+    # Extract indices of colliding particles
+    idx = np.where(collision_mask)[0]
+
+    # Get velocities of colliding ions
+    v_ion = np.column_stack((vx[idx], vy[idx], vz[idx]))
+
+    # Reuse the specific neutrals that we calculated the collision probability against
+    v_neutral = v_neutral_all[idx]
+
+    # Perform Isotropic Elastic Scattering
+    v_ion_new = scatter_isotropic_3d(
+        v1=v_ion, v2=v_neutral, m1=species.m, m2=neutral_mass, rng=rng
+    )
+
+    # 7. Write new velocities back to the species object
+    species.vx[idx] = v_ion_new[:, 0]
+    species.vy[idx] = v_ion_new[:, 1]
+    species.vz[idx] = v_ion_new[:, 2]
+
+    return n_collisions
+
+
+def scatter_isotropic_3d(
+    v1: np.ndarray, v2: np.ndarray, m1: float, m2: float, rng: np.random.Generator
+) -> np.ndarray:
+    """
+    Perform elastic hard-sphere scattering in 3D.
+
+    Updates the velocity of species 1 (v1) assuming isotropic scattering
+    in the Center of Mass (CoM) frame.
+
+    Parameters
+    ----------
+    v1 : np.ndarray
+        Velocities of scattering particles (N, 3).
+    v2 : np.ndarray
+        Velocities of target particles (N, 3).
+    m1 : float
+        Mass of scattering particles.
+    m2 : float
+        Mass of target particles.
+    rng : np.random.Generator
+        Random number generator instance.
+
+    Returns
+    -------
+    np.ndarray
+        New velocities for species 1, shape (N, 3).
+    """
+    n_cols = len(v1)
+    total_mass = m1 + m2
+
+    # 1. Calculate Center of Mass Velocity
+    v_cm = (m1 * v1 + m2 * v2) / total_mass
+
+    # 2. Calculate Relative Velocity (v_rel = v1 - v2)
+    v_rel = v1 - v2
+    v_rel_mag = np.linalg.norm(v_rel, axis=1, keepdims=True)
+
+    # 3. Randomize direction on the unit sphere (Isotropic Scattering)
+    # Pick a random point on a sphere surface:
+    # cos(theta) is uniform [-1, 1], phi is uniform [0, 2pi]
+    cos_theta = 2.0 * rng.random(n_cols) - 1.0
+    sin_theta = np.sqrt(1.0 - cos_theta**2)
+    phi = 2.0 * np.pi * rng.random(n_cols)
+
+    # New direction vector in CoM frame
+    nx = sin_theta * np.cos(phi)
+    ny = sin_theta * np.sin(phi)
+    nz = cos_theta
+
+    # Shape into (N, 3)
+    n_vec = np.column_stack((nx, ny, nz))
+
+    # 4. Calculate new velocity for v1
+    # v1' = v_cm + (m2 / M_tot) * |v_rel| * n_vec
+    return v_cm + (m2 / total_mass) * v_rel_mag * n_vec
diff --git a/src/test_particle_sim_1d/diagnostics.py b/src/test_particle_sim_1d/diagnostics.py
index e69de29..68070c0 100644
--- a/src/test_particle_sim_1d/diagnostics.py
+++ b/src/test_particle_sim_1d/diagnostics.py
@@ -0,0 +1,274 @@
+"""
+diagnostics.py
+
+Tools for calculating physical quantities from particle data in a 1D3V simulation.
+
+This module provides functions to compute:
+1. Global quantities (scalar temperature, mean drift).
+2. Spatially resolved profiles (density n(z), temperature T(z)).
+3. Distribution functions (Energy Distribution Function f(E)).
+
+Functions
+---------
+compute_global_temperature    : Calculate scalar temperature of a species.
+compute_component_temperatures: Calculate temperature for each velocity component.
+compute_drift_velocity        : Calculate mean velocity vector.
+compute_density_profile       : Calculate number density on a 1D spatial grid.
+compute_temperature_profile   : Calculate temperature on a 1D spatial grid.
+compute_energy_distribution   : Calculate the Energy Distribution Function (IEDF/EEDF).
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from . import particles
+from .initialization import constants
+
+
+def compute_global_temperature(species: particles.Species) -> float:
+    """
+    Calculate the global scalar temperature of a particle species.
+
+    Temperature is derived from the variance of the velocity distribution
+    (kinetic energy relative to the mean drift).
+
+    T = (m / 3*kB) * <(v - <v>)^2>
+
+    Parameters
+    ----------
+    species : Species
+        Particle species to analyze.
+
+    Returns
+    -------
+    float
+        Temperature in Kelvin [K]. Returns 0.0 if no particles exist.
+    """
+    if species.N == 0:
+        return 0.0
+
+    # Extract active particle velocities
+    vx = species.vx[: species.N]
+    vy = species.vy[: species.N]
+    vz = species.vz[: species.N]
+
+    # Calculate thermal variances (v - v_drift)^2
+    var_x = np.var(vx)  # np.var calculates mean((x - mean(x))**2)
+    var_y = np.var(vy)
+    var_z = np.var(vz)
+
+    # Total velocity variance (sum of variances for independent dimensions)
+    total_variance = var_x + var_y + var_z
+
+    # T = (m * variance) / (k_B * degrees_of_freedom)
+    # Using 3 degrees of freedom (3V)
+    return (species.m * total_variance) / (3 * constants.kb)
+
+
+def compute_component_temperatures(
+    species: particles.Species,
+) -> tuple[float, float, float]:
+    """
+    Calculate the scalar temperature of a particle species for each velocity component.
+
+    Useful for checking anisotropy (e.g., T_x vs T_y).
+    T_i = (m / kB) * <(v_i - <v_i>)^2>  (1 degree of freedom per component)
+
+    Parameters
+    ----------
+    species : Species
+        Particle species to analyze.
+
+    Returns
+    -------
+    Tuple[float, float, float]
+        (Tx, Ty, Tz) in Kelvin [K].
+    """
+    if species.N == 0:
+        return 0.0, 0.0, 0.0
+
+    # Extract active particle velocities
+    vx = species.vx[: species.N]
+    vy = species.vy[: species.N]
+    vz = species.vz[: species.N]
+
+    # Calculate thermal variances
+    var_x = np.var(vx)
+    var_y = np.var(vy)
+    var_z = np.var(vz)
+
+    factor = species.m / constants.kb
+    return (factor * var_x, factor * var_y, factor * var_z)
+
+
+def compute_drift_velocity(species: particles.Species) -> np.ndarray:
+    """
+    Calculate the global mean drift velocity of a species.
+
+    Parameters
+    ----------
+    species : Species
+        Particle species.
+
+    Returns
+    -------
+    np.ndarray
+        Array [vx_mean, vy_mean, vz_mean] in [m/s].
+    """
+    if species.N == 0:
+        return np.zeros(3)
+
+    return np.array(
+        [
+            np.mean(species.vx[: species.N]),
+            np.mean(species.vy[: species.N]),
+            np.mean(species.vz[: species.N]),
+        ]
+    )
+
+
+def compute_density_profile(
+    species: particles.Species, z_grid: np.ndarray, area: float = 1.0
+) -> np.ndarray:
+    """
+    Calculate the 1D number density profile n(z).
+
+    Parameters
+    ----------
+    species : Species
+        Particle species.
+    z_grid : np.ndarray
+        Edges of the spatial bins (length M+1).
+    area : float, optional
+        Cross-sectional area of the simulation domain [m^2]. Default 1.0.
+
+    Returns
+    -------
+    np.ndarray
+        Number density [m^-3] at each bin center (length M).
+    """
+    if species.N == 0:
+        return np.zeros(len(z_grid) - 1)
+
+    z_p = species.z[: species.N]
+    weights = species.weight[: species.N]
+
+    # Histogram particles onto the grid
+    # "weights" allows for handling super-particles with different weights
+    counts, _ = np.histogram(z_p, bins=z_grid, weights=weights)
+
+    # Calculate bin volumes
+    dz = np.diff(z_grid)
+    volumes = dz * area
+
+    # n = N / Volume
+    return counts / volumes
+
+
+def compute_temperature_profile(
+    species: particles.Species, z_grid: np.ndarray
+) -> np.ndarray:
+    """
+    Calculate the 1D temperature profile T(z).
+
+    Parameters
+    ----------
+    species : Species
+        Particle species.
+    z_grid : np.ndarray
+        Edges of the spatial bins (length M+1).
+
+    Returns
+    -------
+    np.ndarray
+        Temperature [K] at each bin center (length M).
+        Returns 0.0 in bins with < 2 particles.
+    """
+    n_bins = len(z_grid) - 1
+    T_profile = np.zeros(n_bins)
+
+    if species.N == 0:
+        return T_profile
+
+    z_p = species.z[: species.N]
+    vx = species.vx[: species.N]
+    vy = species.vy[: species.N]
+    vz = species.vz[: species.N]
+
+    # Digitize finds which bin index each particle belongs to
+    bin_indices = np.digitize(z_p, z_grid) - 1
+
+    for i in range(n_bins):
+        # Mask for particles in the current bin
+        mask = bin_indices == i
+
+        # Need at least 2 particles to calculate variance/temperature sensibly
+        if np.count_nonzero(mask) < 2:
+            T_profile[i] = 0.0
+            continue
+
+        vx_bin = vx[mask]
+        vy_bin = vy[mask]
+        vz_bin = vz[mask]
+
+        # Calculate local variance (thermal energy)
+        # Using ddof=1 for sample variance if N is small, but np.var default is population
+        # Consistent with global: mean((v - <v>)^2)
+        var = np.var(vx_bin) + np.var(vy_bin) + np.var(vz_bin)
+
+        T_profile[i] = (species.m * var) / (3 * constants.kb)
+
+    return T_profile
+
+
+def compute_energy_distribution(
+    species: particles.Species, n_bins: int = 50, e_max: float | None = None
+) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Calculate the Energy Distribution Function (EDF) of the species.
+
+    Plots f(E) vs E, where E is kinetic energy.
+
+    Parameters
+    ----------
+    species : Species
+        Particle species.
+    n_bins : int, optional
+        Number of histogram bins, default 50.
+    e_max : float | None, optional
+        Maximum energy [eV] to include in histogram. If None, uses max particle energy.
+
+    Returns
+    -------
+    bin_centers : np.ndarray
+        Energy values [eV] (length n_bins).
+    pdf : np.ndarray
+        Probability Density Function value (normalized such that integral is 1).
+    """
+    if species.N == 0:
+        return np.array([]), np.array([])
+
+    vx = species.vx[: species.N]
+    vy = species.vy[: species.N]
+    vz = species.vz[: species.N]
+
+    # Calculate kinetic energy in Joules: 0.5 * m * v^2
+    v2 = vx**2 + vy**2 + vz**2
+    E_joules = 0.5 * species.m * v2
+
+    # Convert to eV
+    E_eV = E_joules / constants.q_e
+
+    # Determine range
+    if e_max is None:
+        e_max = np.max(E_eV)
+
+    # Compute Histogram
+    # density=True normalizes it so the integral over the range is 1
+    hist, bin_edges = np.histogram(E_eV, bins=n_bins, range=(0, e_max), density=True)
+
+    # Calculate bin centers for plotting
+    bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])
+
+    return bin_centers, hist
diff --git a/tests/test_physics_integration.py b/tests/test_physics_integration.py
new file mode 100644
index 0000000..eb51d58
--- /dev/null
+++ b/tests/test_physics_integration.py
@@ -0,0 +1,352 @@
+"""
+test_physics_integration.py
+
+Integration tests for the collision physics and diagnostics modules.
+Combines thermalization, isotropization, and density profile verification.
+
+Tests
+-----
+test_isothermal_relaxation   : Verifies T_ion decays to T_neutral.
+test_velocity_isotropization : Verifies energy redistribution (Tx -> Ty, Tz).
+test_distribution_relaxation : Verifies evolution to Maxwellian energy distribution.
+test_density_profile_reconstruction : Verifies density profile stability under collisions.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from test_particle_sim_1d import collisions, diagnostics
+from test_particle_sim_1d.initialization import constants
+from test_particle_sim_1d.particles import Species
+
+
+def constant_sigma(g: np.ndarray) -> np.ndarray:
+    """Helper: Constant cross-section of 1e-19 m^2."""
+    return np.full_like(g, 1.0e-19)
+
+
+def test_isothermal_relaxation():
+    """
+    Verify that hot ions relax to the neutral gas temperature over time.
+    Case: 0D Box (Fields=0), Hot Ions -> Cold Neutrals.
+    """
+    # --- 1. Simulation Parameters ---
+    n_particles = 20000
+    dt = 1e-8
+    n_steps = 1000
+
+    # Physics Parameters
+    T_ion_init_eV = 1.0  # Start hot (~11600 K)
+    T_neutral_K = 300.0  # Background is cool
+    neutral_density = 2e21
+
+    # --- 2. Initialize Species ---
+    ions = Species(
+        q=constants.q_e, m=constants.m_p, name="proton", capacity=n_particles
+    )
+    ions.initialize_maxwellian(
+        n=n_particles, z_min=0.0, z_max=0.1, temperature={"eV": T_ion_init_eV}, seed=42
+    )
+
+    T_start = diagnostics.compute_global_temperature(ions)
+    print(f"\n[Thermalization] Start T_ion: {T_start:.2f} K (Target: {T_neutral_K} K)")
+
+    assert T_start > 10000, "Initial temperature should be high."
+
+    # --- 3. Run Time Loop ---
+    temp_history = [T_start]
+
+    for _ in range(n_steps):
+        collisions.apply_elastic_collisions(
+            species=ions,
+            neutral_density=neutral_density,
+            neutral_temp_K=T_neutral_K,
+            neutral_mass=constants.m_p,
+            cross_section_func=constant_sigma,
+            dt=dt,
+            seed=None,
+        )
+        current_T = diagnostics.compute_global_temperature(ions)
+        temp_history.append(current_T)
+
+    # --- 4. Verify Results ---
+    T_final = temp_history[-1]
+    print(f"[Thermalization] End   T_ion: {T_final:.2f} K")
+
+    assert T_final < T_start, "Ions should have cooled down."
+
+    # Check convergence to neutral temp (within 20% margin for noise)
+    assert abs(T_final - T_neutral_K) < 0.2 * T_neutral_K, (
+        f"Final T {T_final:.1f} K not close to {T_neutral_K} K"
+    )
+
+    # Check monotonicity (compare start avg vs end avg)
+    avg_start = np.mean(temp_history[:10])
+    avg_end = np.mean(temp_history[-10:])
+    assert avg_end < avg_start
+
+
+def test_velocity_isotropization():
+    """
+    Verify that energy transfers from hot dimensions to cold dimensions.
+    Case: Hot in X, Cold in Y/Z -> Equipartition.
+    """
+    # --- 1. Parameters ---
+    n_particles = 10000
+    dt = 1e-8
+    n_steps = 500
+
+    T_neutral_K = 300.0
+    neutral_density = 2e21
+
+    # --- 2. Initialize Anisotropic Species ---
+    ions = Species(
+        q=constants.q_e, m=constants.m_p, name="proton", capacity=n_particles
+    )
+
+    # Manually create anisotropic distribution
+    rng = np.random.default_rng(42)
+    v_th_hot = 20000.0  # Hot in X
+    v_th_cold = 0.0  # Cold in Y, Z
+
+    vx = rng.normal(0.0, v_th_hot, n_particles)
+    vy = rng.normal(0.0, v_th_cold, n_particles)
+    vz = rng.normal(0.0, v_th_cold, n_particles)
+    z = np.zeros(n_particles)
+
+    ions.add_particles(z, np.column_stack((vx, vy, vz)))
+
+    # Use diagnostics module for component temps
+    Tx_start, Ty_start, Tz_start = diagnostics.compute_component_temperatures(ions)
+    print(
+        f"\n[Isotropization] Start Tx: {Tx_start:.0f}, Ty: {Ty_start:.0f}, Tz: {Tz_start:.0f} K"
+    )
+
+    assert Tx_start > 10000
+    assert Ty_start < 10
+
+    # --- 3. Run Collision Loop ---
+    for _ in range(n_steps):
+        collisions.apply_elastic_collisions(
+            species=ions,
+            neutral_density=neutral_density,
+            neutral_temp_K=T_neutral_K,
+            neutral_mass=constants.m_p,
+            cross_section_func=constant_sigma,
+            dt=dt,
+        )
+
+    # --- 4. Verify Results ---
+    Tx_final, Ty_final, Tz_final = diagnostics.compute_component_temperatures(ions)
+    print(
+        f"[Isotropization] End   Tx: {Tx_final:.0f}, Ty: {Ty_final:.0f}, Tz: {Tz_final:.0f} K"
+    )
+
+    # Y and Z should heat up, X should cool down
+    assert Ty_final > Ty_start
+    assert Tz_final > Tz_start
+    assert Tx_final < Tx_start
+
+    # Check for Equipartition (Isotropy)
+    T_mean = (Tx_final + Ty_final + Tz_final) / 3.0
+
+    assert abs(Tx_final - T_mean) < 0.2 * T_mean, "Tx did not converge to mean"
+    assert abs(Ty_final - T_mean) < 0.2 * T_mean, "Ty did not converge to mean"
+    assert abs(Tz_final - T_mean) < 0.2 * T_mean, "Tz did not converge to mean"
+
+
+def test_distribution_relaxation():
+    """
+    Verify that a mono-energetic beam relaxes to a Maxwellian distribution.
+    Tests: collisions, compute_temperature_profile, compute_energy_distribution.
+    """
+    # --- 1. Parameters ---
+    n_particles = 20000
+    dt = 1e-8
+    n_steps = 1000
+
+    T_neutral_K = 11600.0  # ~1 eV background
+    T_neutral_eV = 1.0
+    neutral_density = 2e21
+
+    # --- 2. Initialize Mono-energetic Species ---
+    ions = Species(
+        q=constants.q_e, m=constants.m_p, name="proton", capacity=n_particles
+    )
+
+    # Initialize uniform positions
+    ions.initialize_uniform(
+        n_particles, z_min=0.0, z_max=0.1, v0=[0.0, 0.0, 0.0], seed=42
+    )
+
+    # Overwrite velocities to be a "Cold Beam"
+    energy_init_eV = 5.0
+    v_mag = np.sqrt(2 * energy_init_eV * constants.q_e / constants.m_p)
+
+    rng = np.random.default_rng(42)
+    phi = rng.uniform(0, 2 * np.pi, n_particles)
+    cos_theta = rng.uniform(-1, 1, n_particles)
+    sin_theta = np.sqrt(1 - cos_theta**2)
+
+    ions.vx[:n_particles] = v_mag * sin_theta * np.cos(phi)
+    ions.vy[:n_particles] = v_mag * sin_theta * np.sin(phi)
+    ions.vz[:n_particles] = v_mag * cos_theta
+
+    print(f"\n[Distribution] Start: Mono-energetic at {energy_init_eV} eV")
+
+    # --- 3. Run Collision Loop ---
+    for _ in range(n_steps):
+        collisions.apply_elastic_collisions(
+            species=ions,
+            neutral_density=neutral_density,
+            neutral_temp_K=T_neutral_K,
+            neutral_mass=constants.m_p,
+            cross_section_func=constant_sigma,
+            dt=dt,
+        )
+
+    # --- 4. Verify Temperature Profile ---
+    z_grid = np.linspace(0.0, 0.1, 11)
+    T_profile = diagnostics.compute_temperature_profile(ions, z_grid)
+    avg_profile_T = np.mean(T_profile)
+    print(
+        f"[Distribution] End Avg Profile T: {avg_profile_T:.0f} K (Target: {T_neutral_K:.0f})"
+    )
+
+    assert abs(avg_profile_T - T_neutral_K) < 0.2 * T_neutral_K, (
+        "Temperature profile average did not converge"
+    )
+
+    # --- 5. Verify Energy Distribution (IEDF) ---
+    # Using 100 bins as requested
+    centers, pdf = diagnostics.compute_energy_distribution(ions, n_bins=100, e_max=5.0)
+
+    kT = T_neutral_eV
+    theory_pdf = (
+        2 * np.sqrt(centers / np.pi) * (1 / kT) ** (1.5) * np.exp(-centers / kT)
+    )
+
+    # --- FIXED L1 NORM CALCULATION ---
+    # 1. Calculate the bin width (Delta E)
+    bin_width = centers[1] - centers[0]
+
+    # 2. Calculate sum of absolute differences
+    diff = np.abs(pdf[1:] - theory_pdf[1:])
+    sum_diff = np.sum(diff)
+
+    # 3. Calculate Area Error (Integral difference)
+    # This value represents the total probability mass that is "misplaced"
+    l1_area_error = sum_diff * bin_width
+
+    print(f"[Distribution] L1 Area Error: {l1_area_error:.3f}")
+
+    # Threshold: 0.15 is generally safe for N=20000 particles and 100 bins.
+    assert l1_area_error < 0.15, (
+        f"L1 Area Error {l1_area_error:.3f} is too high; distribution shape is not Maxwellian."
+    )
+
+
+def test_density_profile_reconstruction():
+    """
+    Verify that compute_density_profile correctly measures a Step-Function density
+    AND that the Collision Module can act on this non-uniform distribution.
+
+    Case: Step Function Density (2:1 Contrast) + Collisional Heating (Cold -> Hot).
+    """
+    # --- 1. Parameters ---
+    n_particles = 30000
+    L = 0.1
+    half_L = L / 2
+    n_bins = 10
+    area = 0.01
+
+    # Collision Parameters
+    dt = 1e-8
+    n_steps = 200
+    T_neutral_K = 1000.0  # Hot neutrals to heat our cold ions
+    neutral_density = 2e21
+
+    # --- 2. Initialize Species ---
+    ions = Species(
+        q=constants.q_e, m=constants.m_p, name="proton", capacity=n_particles
+    )
+
+    # Create Step Function Distribution manually
+    rng = np.random.default_rng(42)
+    n_left = 20000
+    n_right = 10000
+
+    z_left = rng.uniform(0.0, half_L, n_left)
+    z_right = rng.uniform(half_L, L, n_right)
+    z_all = np.concatenate((z_left, z_right))
+
+    # Initialize velocities to ZERO (Cold)
+    v_zeros = np.zeros((n_particles, 3))
+
+    ions.add_particles(z_all, v_zeros)
+
+    # Check Initial Temperature (Should be 0)
+    T_init = diagnostics.compute_global_temperature(ions)
+    print(f"\n[Density+Collisions] Start T: {T_init:.1f} K")
+    assert T_init == 0.0, "Initial temperature should be 0"
+
+    # --- 3. Run Collision Loop (Heating) ---
+    # We do NOT update positions (z), only velocities via collisions.
+    # This ensures density profile should remain constant while T rises.
+
+    for _ in range(n_steps):
+        collisions.apply_elastic_collisions(
+            species=ions,
+            neutral_density=neutral_density,
+            neutral_temp_K=T_neutral_K,
+            neutral_mass=constants.m_p,
+            cross_section_func=constant_sigma,
+            dt=dt,
+        )
+
+    # --- 4. Verify Collision Physics (Heating) ---
+    T_final = diagnostics.compute_global_temperature(ions)
+    print(f"[Density+Collisions] End T:   {T_final:.1f} K (Expected > 50 K)")
+
+    # Particles should have heated up significantly from 0 K
+    assert T_final > 50.0, "Collisions failed to heat the particles"
+
+    # --- 5. Verify Density Profile (Diagnostics) ---
+    # The density profile must remain exactly the step function, verifying that
+    # collisions modified velocities without corrupting positions.
+
+    z_grid = np.linspace(0.0, L, n_bins + 1)
+    dz = L / n_bins
+    bin_volume = dz * area
+
+    density_profile = diagnostics.compute_density_profile(ions, z_grid, area=area)
+
+    expected_density_left = (n_left / 5.0) / bin_volume
+    expected_density_right = (n_right / 5.0) / bin_volume
+
+    avg_density_left = np.mean(density_profile[:5])
+    avg_density_right = np.mean(density_profile[5:])
+
+    print(f"[Density+Collisions] Density Left:  {avg_density_left:.2e}")
+    print(f"[Density+Collisions] Density Right: {avg_density_right:.2e}")
+
+    # Verify Accuracy
+    assert abs(avg_density_left - expected_density_left) < 0.05 * expected_density_left
+    assert (
+        abs(avg_density_right - expected_density_right) < 0.05 * expected_density_right
+    )
+
+    # Verify Contrast (2:1 ratio)
+    assert avg_density_left > 1.8 * avg_density_right, "Density contrast not preserved"
+
+
+if __name__ == "__main__":
+    try:
+        test_isothermal_relaxation()
+        test_velocity_isotropization()
+        test_distribution_relaxation()
+        test_density_profile_reconstruction()
+        print("\nAll Physics Tests Passed!")
+    except AssertionError as e:
+        print(f"\nTest Failed: {e}")