#!/usr/bin/env python3
"""Free Dirac FFT packets: energy sectors, current and coherent motion.

Run with --help. Python standard library + NumPy only; no network or plotting.
Natural units c=hbar=1, m>0. Motion is along z with fixed spin up in the
Dirac basis: the two-component column here embeds as (psi_1,0,psi_2,0).
Thus H(k)=k*sigma_x+m*sigma_z, rho=psi^dagger psi and j_z=psi^dagger sigma_x psi.
This is an invariant free 1D sector, not a complete arbitrary 3D spinor.

An even periodic grid z=-L/2+j*dz uses NumPy FFT(norm='ortho',axis=0):
dz*sum_z rho = dz*sum_k psihat^dagger psihat. The coefficients are discrete
box coefficients, not continuum Fourier amplitudes. U=cos(Et)I-i H sin(Et)/E
evolves every retained Fourier mode EXACTLY; there is no integration step.
The sinc implementation has the correct U=I limit when E=0, but energy
projectors there are undefined. The CLI uses m>0 for a gap and H inverse.

P+/-=(I+/-H/E)/2 split energy sectors. They do NOT delete fixed components.
Negative-sector norm is positive; neither its presence nor interference is a
calculation of vacuum pair creation. Alongside each coherent spinor, CSVs give
the sector-dephased density/current from |psi+><psi+|+|psi-><psi-|. This is a
statistical mixture, not a spinor obtained by adding probabilities.

Position needs particular care on a periodic box. Raw box moments are always
labelled as such. A line-position value is supplied only with a conservative
travel-distance guard, small edge density, and agreement with the independent
continuum Heisenberg displacement computed from retained modes. These are
numerical diagnostics, not a compact-support theorem for projected packets.
The free oscillatory term is representation/observable dependent and is not
a literal microscopic electron trajectory. Fixed-momentum internal tests do
not have localized position expectations.

CSV outputs: profiles, diagnostics, independent box/cutoff/time-difference
refinement, and a near-rest narrow-versus-broad momentum dephasing study.
Every run includes fixed deterministic self-check fixtures. Custom output
parameters do not silently replace those fixtures; inspect their diagnostics.

Physics owners: /relativistic-qm/zitterbewegung/,
/relativistic-qm/hamiltonian-form/, /relativistic-qm/dirac-current/,
/relativistic-qm/dirac-negative-energy-solutions/.
References: Bjorken & Drell, Relativistic Quantum Mechanics (1964);
Thaller, The Dirac Equation (1992). Tested: Python 3.12.14, NumPy 2.3.5.
"""
from __future__ import annotations

import argparse
import csv
from dataclasses import dataclass
from fractions import Fraction
import hashlib
import json
import math
from pathlib import Path
import platform
import sys
import tempfile

import numpy as np

ALPHA = np.array([[0, 1], [1, 0]], dtype=complex)
BETA = np.diag([1, -1]).astype(complex)
IDENTITY = np.eye(2, dtype=complex)
PREPARATIONS = ("positive_only", "negative_only", "upper_component_only", "coherent_mixed")


@dataclass(frozen=True)
class Grid:
    length: float
    points: int
    mass: float

    def __post_init__(self):
        if not math.isfinite(self.length) or self.length <= 0:
            raise ValueError("length must be finite and positive")
        if self.points < 8 or self.points % 2:
            raise ValueError("points must be even and >= 8")
        if not math.isfinite(self.mass) or self.mass <= 0:
            raise ValueError("mass must be positive: the notebook uses gapped energy projectors")

    @property
    def dz(self):
        return self.length / self.points

    @property
    def z(self):
        return -self.length / 2 + self.dz * np.arange(self.points)

    @property
    def k(self):
        return 2 * np.pi * np.fft.fftfreq(self.points, d=self.dz)


def fft(psi):
    return np.fft.fft(psi, axis=0, norm="ortho")


def ifft(coefficients):
    return np.fft.ifft(coefficients, axis=0, norm="ortho")


def act(matrices, vectors):
    return np.einsum("nij,nj->ni", matrices, vectors)


def hamiltonian(k, mass):
    return np.asarray(k)[:, None, None] * ALPHA + mass * BETA


def energy(k, mass):
    return np.hypot(k, mass)


def projectors(k, mass):
    e = energy(k, mass)
    if np.any(e == 0):
        raise ValueError("positive/negative spectral splitting is undefined at E=0")
    h_over_e = hamiltonian(k, mass) / e[:, None, None]
    return (IDENTITY + h_over_e) / 2, (IDENTITY - h_over_e) / 2


def unitary(k, mass, time):
    e = energy(k, mass)
    # sinc(x)=sin(pi*x)/(pi*x), including x=0: no 0/0 at the massless apex.
    sine_over_e = time * np.sinc(e * time / np.pi)
    return np.cos(e * time)[:, None, None] * IDENTITY - 1j * sine_over_e[:, None, None] * hamiltonian(k, mass)


def evolve(psi, grid, time):
    return ifft(act(unitary(grid.k, grid.mass, time), fft(psi)))


def norm_squared(psi, grid):
    return float(grid.dz * np.vdot(psi, psi).real)


def relative_error(actual, expected):
    return float(np.linalg.norm(actual - expected) / max(np.linalg.norm(expected), 1e-300))


def normalize(psi, grid):
    size = norm_squared(psi, grid)
    if not math.isfinite(size) or size < 1e-28:
        raise ValueError("preparation has zero or unresolved norm")
    return psi / math.sqrt(size)


def prepare(grid, kind, center=-6.0, width=2.0, momentum=0.9):
    y = grid.z - center
    envelope = (2 * np.pi * width ** 2) ** (-0.25) * np.exp(-y ** 2 / (4 * width ** 2) + 1j * momentum * y)
    spin = np.array([1, 1], dtype=complex) / math.sqrt(2) if kind == "coherent_mixed" else np.array([1, 0], dtype=complex)
    if kind == "negative_only":
        spin = np.array([0, 1], dtype=complex)
    psi = envelope[:, None] * spin
    if kind in ("positive_only", "negative_only"):
        plus, minus = projectors(grid.k, grid.mass)
        psi = ifft(act(plus if kind == "positive_only" else minus, fft(psi)))
    elif kind not in PREPARATIONS:
        raise ValueError(f"unknown preparation: {kind}")
    return normalize(psi, grid)


def sector_coefficients(psi, grid):
    plus, minus = projectors(grid.k, grid.mass)
    coefficients = fft(psi)
    return act(plus, coefficients), act(minus, coefficients)


def density_current(psi):
    density = np.sum(np.abs(psi) ** 2, axis=1)
    current = 2 * np.real(np.conj(psi[:, 0]) * psi[:, 1])
    return density, current


def observable_expectation(coefficients, matrices, grid):
    return complex(grid.dz * np.einsum("ni,ni->", coefficients.conj(), act(matrices, coefficients)))


def heisenberg_operators(k, mass, time):
    """Independent ordered velocity/displacement formula from the owner.

    A=alpha-k H^-1 anticommutes with H.  Keep A left of exp(-2iHt).
    This is a continuum position-operator identity, not an identity for the
    discontinuous multiplication-by-z observable on a periodic grid.
    """
    e = energy(k, mass)
    if np.any(e == 0):
        raise ValueError("H inverse is undefined at E=0")
    h_inv = hamiltonian(k, mass) / e[:, None, None] ** 2
    drift = np.asarray(k)[:, None, None] * h_inv
    oscillatory = ALPHA - drift
    twice_time = unitary(k, mass, 2 * time)
    velocity = drift + oscillatory @ twice_time
    displacement = drift * time + 0.5j * oscillatory @ h_inv @ (twice_time - IDENTITY)
    return velocity, displacement


def packet_state(initial, grid, time, center, width):
    """Diagnostics of a unit-normalized initial packet and its dephased mixture."""
    if abs(norm_squared(initial, grid) - 1) > 1e-10:
        raise ValueError("packet_state requires an initial packet normalized to one")
    coefficients = fft(initial)
    initial_density, _ = density_current(initial)
    initial_center = float(grid.dz * np.sum(grid.z * initial_density))
    psi = evolve(initial, grid, time)
    density, current = density_current(psi)
    a, b = sector_coefficients(initial, grid)
    plus = ifft(a * np.exp(-1j * energy(grid.k, grid.mass) * time)[:, None])
    minus = ifft(b * np.exp(1j * energy(grid.k, grid.mass) * time)[:, None])
    plus_density, plus_current = density_current(plus)
    minus_density, minus_current = density_current(minus)
    norm = float(grid.dz * np.sum(density))
    raw_center = float(grid.dz * np.sum(grid.z * density) / norm)
    raw_variance = float(grid.dz * np.sum((grid.z - raw_center) ** 2 * density) / norm)
    edge_fraction = float(grid.dz * np.sum(density[np.abs(grid.z) >= 0.4 * grid.length]) / norm)
    velocity_op, displacement_op = heisenberg_operators(grid.k, grid.mass, time)
    predicted = initial_center + observable_expectation(coefficients, displacement_op, grid).real
    discrepancy = abs(raw_center - predicted)
    guard = abs(center) + 8 * width + abs(time) < grid.length / 2
    line_valid = bool(guard and edge_fraction < 1e-10 and discrepancy < 2e-8)
    h2_weights = energy(grid.k, grid.mass) ** 2 * np.sum(np.abs(coefficients) ** 2, axis=1)
    row = {
        "time": float(time), "norm": norm,
        "positive_sector_norm": norm_squared(a, grid), "negative_sector_norm": norm_squared(b, grid),
        "upper_component_norm": float(grid.dz * np.sum(np.abs(psi[:, 0]) ** 2)),
        "lower_component_norm": float(grid.dz * np.sum(np.abs(psi[:, 1]) ** 2)),
        "energy_expectation": observable_expectation(fft(psi), hamiltonian(grid.k, grid.mass), grid).real,
        "mean_current_velocity": float(grid.dz * np.sum(current)),
        "sector_dephased_mean_velocity": float(grid.dz * np.sum(plus_current + minus_current)),
        "interference_mean_velocity": float(grid.dz * np.sum(current - plus_current - minus_current)),
        "heisenberg_mean_velocity": observable_expectation(coefficients, velocity_op, grid).real,
        "raw_periodic_box_centroid": raw_center,
        "raw_periodic_box_standard_deviation": math.sqrt(max(0, raw_variance)),
        "continuum_heisenberg_centroid_prediction": predicted,
        "box_versus_heisenberg_centroid_error": discrepancy,
        "conservative_travel_guard_passed": bool(guard),
        "line_position_diagnostics_passed": line_valid,
        "line_position_if_diagnostics_pass": raw_center if line_valid else None,
        "edge_probability_outer_ten_percent_each_side": edge_fraction,
        "h_squared_weight_above_80_percent_nyquist": float(np.sum(h2_weights[
            np.abs(grid.k) > 0.8 * np.pi / grid.dz]) / np.sum(h2_weights)),
        "max_pointwise_current_bound_violation": float(np.max(np.maximum(0, np.abs(current) - density))),
    }
    return psi, density, current, plus_density + minus_density, plus_current + minus_current, row


class Checks:
    def __init__(self):
        self.rows = []

    def bound(self, name, value, tolerance):
        self.rows.append({"name": name, "value": float(value), "absolute_tolerance": tolerance,
                          "passed": bool(math.isfinite(value) and abs(value) <= tolerance)})

    def condition(self, name, condition, evidence=None):
        self.rows.append({"name": name, "passed": bool(condition), "evidence": evidence})

    @property
    def passed(self):
        return all(row["passed"] for row in self.rows)


def exact_fixture(checks):
    f = Fraction
    def mul(a, b):
        return [[sum(a[i][r] * b[r][j] for r in range(2)) for j in range(2)] for i in range(2)]
    h = [[f(4, 5), f(3, 5)], [f(3, 5), f(-4, 5)]]
    alpha = [[f(0), f(1)], [f(1), f(0)]]
    ident = [[f(1), f(0)], [f(0), f(1)]]
    plus = [[(ident[i][j] + h[i][j]) / 2 for j in range(2)] for i in range(2)]
    minus = [[(ident[i][j] - h[i][j]) / 2 for j in range(2)] for i in range(2)]
    checks.condition("exact rational H squared", mul(h, h) == ident, "m=4/5, k=3/5, E=1")
    for name, p, sign in (("plus", plus, 1), ("minus", minus, -1)):
        checks.condition(f"exact rational {name} projector", mul(p, p) == p)
        expected = [[sign * f(3, 5) * p[i][j] for j in range(2)] for i in range(2)]
        checks.condition(f"exact rational {name} projected velocity", mul(mul(p, alpha), p) == expected)
    checks.condition("exact rational orthogonal sectors", mul(plus, minus) == [[0, 0], [0, 0]])


def structural_checks(checks):
    exact_fixture(checks)
    k = np.array([-1e6, -7, -0.9, -1e-10, 0, 1e-10, 0.9, 7, 1e6])
    for mass in (0.01, 1, 100):
        h, e = hamiltonian(k, mass), energy(k, mass)
        plus, minus = projectors(k, mass)
        checks.bound(f"m={mass}: H Hermitian", float(np.max(np.abs(h - h.conj().transpose(0, 2, 1)))), 0)
        checks.bound(f"m={mass}: scaled H squared", float(np.max(np.abs(h @ h / e[:, None, None] ** 2 - IDENTITY))), 8e-15)
        eigenvalues, eigenvectors = np.linalg.eigh(h)
        numerical_plus = eigenvectors[:, :, 1, None] * eigenvectors[:, None, :, 1].conj()
        checks.bound(f"m={mass}: independent eigensolver projectors", float(np.max(np.abs(plus - numerical_plus))), 8e-15)
        for name, p, sign in (("plus", plus, 1), ("minus", minus, -1)):
            checks.bound(f"m={mass}: {name} idempotence", float(np.max(np.abs(p @ p - p))), 8e-15)
            checks.bound(f"m={mass}: {name} projected velocity", float(np.max(np.abs(
                p @ ALPHA @ p - sign * (k / e)[:, None, None] * p))), 8e-15)
        checks.bound(f"m={mass}: sector orthogonality", float(np.max(np.abs(plus @ minus))), 8e-15)
        for time in (0, 0.17, -2.3):
            u = unitary(k, mass, time)
            checks.bound(f"m={mass},t={time}: unitarity", float(np.max(np.abs(u.conj().transpose(0, 2, 1) @ u - IDENTITY))), 2e-14)
            phase = plus * np.exp(-1j * e * time)[:, None, None] + minus * np.exp(1j * e * time)[:, None, None]
            checks.bound(f"m={mass},t={time}: independent spectral evolution", float(np.max(np.abs(u - phase))), 2e-14)
    checks.bound("massless E=0 exact evolution limit", float(np.max(np.abs(unitary(np.array([0.0]), 0, 7) - IDENTITY))), 0)
    try:
        projectors(np.array([0.0]), 0)
    except ValueError:
        checks.condition("massless zero mode refuses undefined sector split", True)
    else:
        checks.condition("massless zero mode refuses undefined sector split", False)
    grid = Grid(128, 1024, 1)
    summaries = {}
    for kind in PREPARATIONS:
        initial = prepare(grid, kind)
        a, b = sector_coefficients(initial, grid)
        e = energy(grid.k, grid.mass)
        initial_energy = observable_expectation(fft(initial), hamiltonian(grid.k, grid.mass), grid).real
        norm_error = sector_error = continuity_error = current_error = velocity_error = position_error = phase_error = mixture_error = 0.0
        drift = float(grid.dz * np.sum(grid.k / e * (
            np.sum(np.abs(a) ** 2, axis=1) - np.sum(np.abs(b) ** 2, axis=1))))
        max_interference = 0.0
        for time in (0, 0.3, 2, 7, 24):
            psi, rho, current, dephased_density, dephased_current, row = packet_state(initial, grid, time, -6, 2)
            independent = ifft(a * np.exp(-1j * e * time)[:, None] + b * np.exp(1j * e * time)[:, None])
            phase_error = max(phase_error, relative_error(psi, independent))
            norm_error = max(norm_error, abs(row["norm"] - 1), abs(row["energy_expectation"] - initial_energy))
            at, bt = sector_coefficients(psi, grid)
            sector_error = max(sector_error, abs(norm_squared(at, grid) - norm_squared(a, grid)),
                               abs(norm_squared(bt, grid) - norm_squared(b, grid)))
            dot = ifft(-1j * act(hamiltonian(grid.k, grid.mass), fft(psi)))
            rho_dot = 2 * np.real(np.sum(np.conj(psi) * dot, axis=1))
            current_divergence = ifft(1j * grid.k[:, None] * fft(current[:, None]))[:, 0].real
            continuity_error = max(continuity_error, float(np.max(np.abs(rho_dot + current_divergence))))
            current_error = max(current_error, row["max_pointwise_current_bound_violation"])
            velocity_error = max(velocity_error, abs(row["mean_current_velocity"] - row["heisenberg_mean_velocity"]))
            position_error = max(position_error, row["box_versus_heisenberg_centroid_error"])
            max_interference = max(max_interference, abs(row["interference_mean_velocity"]))
            mixture_error = max(mixture_error, abs(grid.dz * np.sum(dephased_density) - 1),
                abs(grid.dz * np.sum(dephased_current) - drift),
                float(np.max(np.maximum(0, np.abs(dephased_current) - dephased_density))))
            checks.condition(f"{kind},t={time}: position diagnostics controlled", row["line_position_diagnostics_passed"],
                             row["edge_probability_outer_ten_percent_each_side"])
        for label, value, tolerance in (("norm and energy conservation", norm_error, 2e-13),
                ("fixed sector populations", sector_error, 2e-13), ("independent branch evolution", phase_error, 2e-13),
                ("local continuity", continuity_error, 2e-11), ("pointwise abs(j)<=rho", current_error, 2e-14),
                ("Heisenberg velocity", velocity_error, 2e-13), ("Heisenberg position with wrap control", position_error, 2e-9),
                ("dephased mixture norm, current bound and constant drift", mixture_error, 2e-13)):
            checks.bound(f"{kind}: {label}", value, tolerance)
        if kind in ("positive_only", "negative_only"):
            checks.bound(f"{kind}: no energy-sector interference velocity", max_interference, 2e-13)
        if kind == "positive_only":
            lower = float(grid.dz * np.sum(np.abs(initial[:, 1]) ** 2))
            checks.condition("positive energy is not upper-component-only", lower > 0.01 and norm_squared(b, grid) < 1e-25,
                             {"lower_component_weight": lower, "negative_sector_weight": norm_squared(b, grid)})
        if kind == "upper_component_only":
            checks.condition("upper-component-only is not positive energy", norm_squared(b, grid) > 0.01,
                             {"negative_sector_weight": norm_squared(b, grid)})
        summaries[kind] = {"positive_norm": norm_squared(a, grid), "negative_norm": norm_squared(b, grid),
                           "maximum_interference_velocity": max_interference}
    # An exactly solvable internal k=0 fixture; no localized plane-wave X is claimed.
    spin = np.array([[1, 1]], dtype=complex) / math.sqrt(2)
    for time in (0, 0.2, 0.7, 2.1):
        v, delta = heisenberg_operators(np.array([0.0]), 1, time)
        checks.bound(f"rest internal velocity cos(2mt),t={time}", float((spin.conj() * act(v, spin)).sum().real) - math.cos(2 * time), 2e-14)
        checks.bound(f"rest internal displacement matrix,t={time}", float((spin.conj() * act(delta, spin)).sum().real) - math.sin(2 * time) / 2, 2e-14)
    return summaries


def refinement_studies(checks):
    rows = []
    for kind in ("positive_only", "coherent_mixed"):
        reference_grid = Grid(128, 2048, 1)
        reference = evolve(prepare(reference_grid, kind), reference_grid, 24)
        errors = []
        for points in (64, 128, 256, 512, 1024):
            grid = Grid(128, points, 1)
            result = evolve(prepare(grid, kind), grid, 24)
            error = relative_error(result, reference[::2048 // points])
            errors.append(error)
            rows.append({"study": "cutoff_fixed_box", "preparation": kind, "length": 128, "points": points,
                         "dz": grid.dz, "nyquist_abs_k": np.pi / grid.dz, "time": 24, "relative_error": error,
                         "region": "entire same periodic box"})
        checks.bound(f"{kind}: refined cutoff error", errors[-1], 2e-11)
        checks.condition(f"{kind}: cutoff refinement detects unresolved coarse grid",
                         errors[0] > 1e-4 and errors[1] < errors[0] / 30, errors)
        reference_grid = Grid(512, 4096, 1)
        reference = evolve(prepare(reference_grid, kind), reference_grid, 24)
        reference_mask = np.abs(reference_grid.z) <= 12
        errors = []
        for length in (32, 64, 128, 256):
            grid = Grid(length, int(length / 0.125), 1)
            mask = np.abs(grid.z) <= 12
            assert np.array_equal(grid.z[mask], reference_grid.z[reference_mask])
            initial = prepare(grid, kind)
            result = evolve(initial, grid, 24)
            error = relative_error(result[mask], reference[reference_mask])
            errors.append(error)
            state = packet_state(initial, grid, 24, -6, 2)[-1]
            rows.append({"study": "box_fixed_cutoff", "preparation": kind, "length": length, "points": grid.points,
                         "dz": grid.dz, "nyquist_abs_k": np.pi / grid.dz, "time": 24, "relative_error": error,
                         "region": "fixed -12<=z<=12", "edge_probability": state["edge_probability_outer_ten_percent_each_side"],
                         "line_position_diagnostics_passed": state["line_position_diagnostics_passed"],
                         "box_versus_heisenberg_centroid_error": state["box_versus_heisenberg_centroid_error"]})
            if length == 32:
                checks.condition(f"{kind}: small-box position not mislabelled as line position",
                                 not state["line_position_diagnostics_passed"], state["box_versus_heisenberg_centroid_error"])
        checks.bound(f"{kind}: refined box error", errors[-1], 2e-11)
        checks.condition(f"{kind}: box refinement detects finite-box effects",
                         errors[0] > 1e-6 and errors[1] < errors[0] / 30, errors)
    grid = Grid(128, 1024, 1)
    initial = prepare(grid, "coherent_mixed")
    time = 3.2
    psi = evolve(initial, grid, time)
    hpsi = ifft(act(hamiltonian(grid.k, grid.mass), fft(psi)))
    errors = []
    for dt in (0.08, 0.04, 0.02, 0.01):
        numerical_derivative = (evolve(initial, grid, time + dt) - evolve(initial, grid, time - dt)) / (2 * dt)
        residual = float(np.linalg.norm(1j * numerical_derivative - hpsi) / np.linalg.norm(hpsi))
        errors.append(residual)
        rows.append({"study": "time_difference_not_integration", "preparation": "coherent_mixed", "length": 128,
                     "points": 1024, "time": time, "diagnostic_dt": dt, "relative_dirac_equation_residual": residual})
    ratios = [a / b for a, b in zip(errors, errors[1:])]
    checks.condition("time-difference Dirac residual is quadratic", all(3.8 < r < 4.2 for r in ratios), ratios)
    checks.bound("finest time-difference residual", errors[-1], 5e-5)
    return rows


def dephasing_study(checks):
    """Measure coherent spectral envelope, not oscillation extrema sampled coarsely."""
    grid = Grid(512, 4096, 1)
    rows, endpoints = [], {}
    for width in (2.0, 8.0):
        initial = prepare(grid, "coherent_mixed", center=0, width=width, momentum=0)
        a, b = sector_coefficients(initial, grid)
        e = energy(grid.k, grid.mass)
        coefficients = np.einsum("ni,ij,nj->n", a.conj(), ALPHA, b) * grid.dz
        c0 = complex(np.sum(coefficients))
        weights_a = np.sum(np.abs(a) ** 2, axis=1)
        weights_b = np.sum(np.abs(b) ** 2, axis=1)
        drift = float(grid.dz * np.sum(grid.k / e * (weights_a - weights_b)))
        for time in np.linspace(0, 128, 257):
            coherence = complex(np.sum(coefficients * np.exp(2j * e * time)))
            envelope = abs(coherence) / abs(c0)
            rows.append({"position_width": width, "momentum_width": 1 / (2 * width), "time": float(time),
                         "coherence_real": coherence.real, "coherence_imag": coherence.imag,
                         "relative_coherence_envelope": envelope, "sector_dephased_mean_velocity": drift,
                         "coherent_mean_velocity": drift + 2 * coherence.real})
        endpoints[width] = envelope
        for time in (0, 2, 24, 64, 128):
            state = packet_state(initial, grid, time, 0, width)[-1]
            expected = drift + 2 * np.sum(coefficients * np.exp(2j * e * time)).real
            checks.bound(f"width={width},t={time}: spectral coherence versus real-space current",
                         state["mean_current_velocity"] - expected, 2e-13)
            checks.bound(f"width={width},t={time}: controlled packet Heisenberg centroid",
                         state["box_versus_heisenberg_centroid_error"], 2e-8)
            checks.condition(f"width={width},t={time}: line-position guard", state["line_position_diagnostics_passed"])
    checks.condition("broader momentum packet dephases earlier", endpoints[2.0] < 0.4 and endpoints[8.0] > 0.75, endpoints)
    return rows


def write_csv(path, rows):
    fields = list(dict.fromkeys(key for row in rows for key in row))
    with path.open("x", encoding="utf8", newline="") as stream:
        writer = csv.DictWriter(stream, fieldnames=fields)
        writer.writeheader()
        writer.writerows(rows)


def write_profiles(path, grid, args):
    diagnostics = []
    columns = ["preparation", "time", "z", "upper_real", "upper_imag", "lower_real", "lower_imag",
               "density", "current", "sector_dephased_density", "sector_dephased_current", "interference_density", "interference_current"]
    with path.open("x", encoding="utf8", newline="") as stream:
        writer = csv.writer(stream)
        writer.writerow(columns)
        for kind in PREPARATIONS:
            initial = prepare(grid, kind, args.center, args.width, args.momentum)
            for time in np.linspace(0, args.final_time, args.time_samples):
                psi, rho, current, dephased_rho, dephased_j, row = packet_state(initial, grid, float(time), args.center, args.width)
                diagnostics.append({"preparation": kind, **row})
                for i, z in enumerate(grid.z):
                    writer.writerow([kind, time, z, psi[i, 0].real, psi[i, 0].imag, psi[i, 1].real, psi[i, 1].imag,
                                     rho[i], current[i], dephased_rho[i], dephased_j[i], rho[i] - dephased_rho[i], current[i] - dephased_j[i]])
    return diagnostics


def parser():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument("--output-dir", type=Path, help="new CSV/JSON directory; default: new temporary directory")
    p.add_argument("--json-output", type=Path, help="override summary JSON path")
    p.add_argument("--profiles-csv", type=Path, help="override profiles CSV path")
    p.add_argument("--mass", type=float, default=1.0, help="m>0 (default 1)")
    p.add_argument("--length", type=float, default=128.0, help="periodic box length (default 128)")
    p.add_argument("--points", type=int, default=1024, help="even spatial sample count (default 1024)")
    p.add_argument("--width", type=float, default=2.0, help="unprojected Gaussian intensity standard deviation (default 2)")
    p.add_argument("--center", type=float, default=-6.0, help="Gaussian center (default -6)")
    p.add_argument("--momentum", type=float, default=0.9, help="Gaussian carrier momentum (default 0.9)")
    p.add_argument("--final-time", type=float, default=24.0, help="last exact output time, not an integration step (default 24)")
    p.add_argument("--time-samples", type=int, default=5, help="number of output times (default 5)")
    return p


def main(argv=None):
    p = parser()
    args = p.parse_args(argv)
    try:
        grid = Grid(args.length, args.points, args.mass)
        if not all(math.isfinite(value) for value in (args.width, args.center, args.momentum, args.final_time)):
            raise ValueError("all numerical arguments must be finite")
        if args.width < 2 * grid.dz or args.time_samples < 2:
            raise ValueError("width must span at least two grid intervals and time-samples must be >=2")
        if abs(args.center) + 4 * args.width >= grid.length / 2:
            raise ValueError("center plus four Gaussian widths must lie inside the box")
        if abs(args.momentum) >= np.pi / grid.dz:
            raise ValueError("carrier momentum must lie strictly inside the Nyquist cutoff")
    except ValueError as error:
        p.error(str(error))
    output = args.output_dir or Path(tempfile.mkdtemp(prefix="dirac-packet-"))
    paths = {"profiles": args.profiles_csv or output / "profiles.csv", "diagnostics": output / "diagnostics.csv",
             "refinement": output / "refinement.csv", "dephasing": output / "dephasing.csv",
             "summary": args.json_output or output / "summary.json"}
    if len({path.resolve() for path in paths.values()}) != len(paths):
        p.error("output paths must be distinct")
    for path in paths.values():
        if path.exists() or path.is_symlink():
            p.error(f"existing outputs are never overwritten: {path}")
        if ".git" in [part.lower() for part in path.resolve().parts]:
            p.error(f"output cannot be inside .git: {path}")
    for path in paths.values():
        path.parent.mkdir(parents=True, exist_ok=True)
    checks = Checks()
    sector_report = structural_checks(checks)
    refinement = refinement_studies(checks)
    dephasing = dephasing_study(checks)
    diagnostic_rows = write_profiles(paths["profiles"], grid, args)
    write_csv(paths["diagnostics"], diagnostic_rows)
    write_csv(paths["refinement"], refinement)
    write_csv(paths["dephasing"], dephasing)
    report = {
        "program": "dirac-wave-packets.py", "passed": checks.passed,
        "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
        "environment": {"python": platform.python_version(), "numpy": np.__version__, "platform": platform.platform()},
        "randomness": "none", "units": "c=hbar=1; m>0; 1D spin-up z-motion sector",
        "basis": "(upper,lower) embeds as (upper,0,lower,0); alpha=sigma_x, beta=sigma_z",
        "method": "unitary spatial FFT and exact finite-mode cosine/sine evolution",
        "time_step_error": "not applicable; dt refinement tests a numerical derivative only",
        "parameters": {key: getattr(args, key) for key in ("mass", "length", "points", "width", "center", "momentum", "final_time", "time_samples")},
        "interpretation": ["Energy sectors are momentum-dependent operator projectors, not fixed components.",
            "Both sector populations are positive Hilbert-space norms, not pair-production probabilities.",
            "Sector dephasing is a statistical mixture; coherent fields retain local interference.",
            "Raw periodic-box moments are not unwrapped line positions after wraparound.",
            "Line moments require all reported diagnostics; the travel guard is conservative and is not a compact-support theorem.",
            "Zitterbewegung uses the ordinary Dirac position representation and is not a literal electron trajectory.",
            "Fixed self-check fixtures are separate from configurable profile parameters; inspect the latter's edge/cutoff diagnostics."],
        "sector_fixtures": sector_report, "checks": checks.rows,
        "output_files": {key: str(path) for key, path in paths.items()},
        "output_sha256": {key: hashlib.sha256(path.read_bytes()).hexdigest() for key, path in paths.items() if key != "summary"},
    }
    with paths["summary"].open("x", encoding="utf8") as stream:
        json.dump(report, stream, indent=2, allow_nan=False)
        stream.write("\n")
    print(json.dumps({"passed": checks.passed, "checks": len(checks.rows),
                      "failed": [row for row in checks.rows if not row["passed"]], "summary_json": str(paths["summary"])}, indent=2))
    return 0 if checks.passed else 1


if __name__ == "__main__":
    sys.exit(main())