#!/usr/bin/env python3
"""Dirac step scattering: channel choice, spinor matching and conserved flux.

Run with --help. Standard library + NumPy only; offline and deterministic.
Natural units hbar=c=1, H=-i*sigma_x*d/dx+m*sigma_z+V(x), one fixed spin
channel at normal incidence. V is potential ENERGY, not voltage. The step
is V_L for x<0 and V_R for x>0; the conserved total energy is E. We assume
E-V_L>m and impose incidence only from the left.

For a propagating channel, kinetic energy eps=E-V obeys eps^2=k^2+m^2.
Outgoing means positive current/velocity k/eps on the right. In the Klein
region eps<-m this requires k<0. Real modes are normalized to unit ABSOLUTE
Dirac probability current; incoming/reflected/transmitted flux magnitudes
then give R=|r|^2,T=|tau|^2. The decaying gap mode has k=+i*kappa and zero
current, and uses unit density at the interface instead of flux normalization.
Its coefficient is not a transmitted-flux amplitude. Exact thresholds are
excluded and explicitly recorded: zero-current channels cannot be unit-flux
normalized. A stable spinor chart is used on each side of the mass gap.

The algorithm solves a 2x2 continuity system from explicit spinors. Independent
closed amplitudes, full two-port S-matrix unitarity, direct current evaluation,
Hermitian eigensolver projectors, exact rational fixtures, threshold approaches,
finite-difference stationary residuals, and nonrelativistic/massless limits
provide checks. No time propagation, spatial mesh scattering solve, or repeated
transfer-matrix product is used. CSV positions only sample analytic solutions.

This is a prescribed, abrupt, stationary ONE-BODY Dirac step. R+T=1 refers to
positive Dirac probability current with the stated boundary condition. It is
not the signed scalar Klein--Gordon current convention and is NOT a vacuum
pair-production calculation. A pair count requires quantized in/out modes,
initial occupations and a specified source preparation. No such count is
exported. Smooth interfaces, transverse momentum, recoil, finite barriers,
radiation and backreaction are outside this experiment.

Outputs: step sweep, spatial spinor profiles, threshold approaches, derivative
refinement and limiting-case CSVs; JSON records versions, checks and hashes.
Existing outputs are never overwritten; default is a new temporary directory.
Owner: /relativistic-qm/klein-paradox/.
References: Calogeracos & Dombey, Contemporary Physics 40, 313 (1999);
Gavrilov & Gitman, Physical Review D 93, 045002 (2016).
Tested environment: Python 3.12.14, NumPy 2.3.5.
"""
from __future__ import annotations

import argparse
import csv
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]], complex)
BETA = np.diag([1, -1]).astype(complex)


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

    def condition(self, name, passed, **details):
        self.rows.append({"name": name, "passed": bool(passed), **details})

    def close(self, name, actual, expected, tolerance=4e-12):
        absolute = float(np.linalg.norm(np.asarray(actual) - np.asarray(expected)))
        scale = max(1.0, float(np.linalg.norm(expected)))
        self.condition(name, absolute <= tolerance * scale, absolute_error=absolute,
                       scaled_error=absolute / scale, tolerance=tolerance)

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


def current(spinor):
    return float(np.vdot(spinor, ALPHA @ spinor).real)


def density(spinor):
    return float(np.vdot(spinor, spinor).real)


def channel(epsilon, mass, direction=1):
    """Select by velocity; direction is used only for real propagating modes."""
    threshold_distance = abs(abs(epsilon) - mass)
    if threshold_distance <= 1e-12 * max(1.0, abs(epsilon), mass):
        raise ValueError("threshold excluded: zero-current/asymptotically unresolved channel")
    if abs(epsilon) > mass:
        momentum = direction * math.copysign(1.0, epsilon) * math.sqrt((abs(epsilon) - mass) * (abs(epsilon) + mass))
        region = "positive_continuum" if epsilon > 0 else "Klein_negative_continuum"
    else:
        momentum = 1j * math.sqrt((mass - epsilon) * (mass + epsilon))
        region = "evanescent_gap"
    # Both charts solve (k*sigma_x+m*sigma_z)w=epsilon*w; choose the one
    # that avoids the small epsilon+m denominator near the negative threshold.
    spinor = np.array([epsilon + mass, momentum], complex) if epsilon >= 0 else np.array([momentum, epsilon - mass], complex)
    spinor /= np.linalg.norm(spinor)
    spinor *= spinor[0].conjugate() / abs(spinor[0])
    normalization = "unit_interface_density"
    if region != "evanescent_gap":
        spinor /= math.sqrt(abs(current(spinor)))
        normalization = "unit_absolute_probability_current"
    return {"epsilon": epsilon, "k": momentum, "spinor": spinor, "region": region,
            "normalization": normalization, "current": current(spinor), "density": density(spinor)}


def match_step(energy, mass, right_potential, left_potential=0.0):
    left_energy = energy - left_potential
    if left_energy <= mass:
        raise ValueError("left incident channel must have positive kinetic energy above the mass gap")
    incoming = channel(left_energy, mass, 1)
    reflected = channel(left_energy, mass, -1)
    transmitted = channel(energy - right_potential, mass, 1)
    matrix = np.column_stack((reflected["spinor"], -transmitted["spinor"]))
    reflection, transmission = np.linalg.solve(matrix, -incoming["spinor"])
    is_gap = transmitted["region"] == "evanescent_gap"
    return {"energy": energy, "mass": mass, "left_potential": left_potential, "right_potential": right_potential,
            "incoming": incoming, "reflected": reflected, "transmitted": transmitted,
            "r": reflection, "tau": transmission, "R": float(abs(reflection) ** 2),
            "T": 0.0 if is_gap else float(abs(transmission) ** 2),
            "matching_condition_number": float(np.linalg.cond(matrix))}


def wave(solution, x, side):
    if side == "left":
        inc, ref = solution["incoming"], solution["reflected"]
        return inc["spinor"] * np.exp(1j * inc["k"] * x) + solution["r"] * ref["spinor"] * np.exp(1j * ref["k"] * x)
    out = solution["transmitted"]
    return solution["tau"] * out["spinor"] * np.exp(1j * out["k"] * x)


def closed_coefficients(solution):
    left, right, mass = solution["incoming"], solution["transmitted"], solution["mass"]
    a = left["k"] / (left["epsilon"] + mass)
    b = right["k"] / (right["epsilon"] + mass)
    reflection = (a - b) / (a + b)
    chart_transmission = 2 * a / (a + b)
    # Convert the owner's upper-component-one chart to each declared normalization.
    transmission = chart_transmission * left["spinor"][0] / right["spinor"][0]
    return reflection, transmission, a, b


def scattering_matrix(energy, mass, right_potential, left_potential=0.0):
    left_in = channel(energy - left_potential, mass, 1)["spinor"]
    left_out = channel(energy - left_potential, mass, -1)["spinor"]
    right_out = channel(energy - right_potential, mass, 1)
    if right_out["region"] == "evanescent_gap":
        raise ValueError("the evanescent gap has no right asymptotic propagating port")
    right_in = channel(energy - right_potential, mass, -1)["spinor"]
    return np.linalg.solve(np.column_stack((left_out, -right_out["spinor"])), np.column_stack((-left_in, right_in)))


def check_solution(checks, solution, label):
    inc, ref, out = (solution[key] for key in ("incoming", "reflected", "transmitted"))
    mass = solution["mass"]
    checks.close(label + " continuity", wave(solution, 0, "left"), wave(solution, 0, "right"))
    checks.close(label + " incident current", inc["current"], 1)
    checks.close(label + " reflected unit-mode current", ref["current"], -1)
    checks.close(label + " outgoing or gap current", out["current"], 0 if out["region"] == "evanescent_gap" else 1)
    for name, mode in (("incident", inc), ("reflected", ref), ("right", out)):
        h = mode["k"] * ALPHA + mass * BETA
        checks.close(label + " " + name + " eigen-equation", h @ mode["spinor"], mode["epsilon"] * mode["spinor"])
        checks.condition(label + " " + name + " positive density", mode["density"] > 0)
        if mode["region"] != "evanescent_gap":
            checks.close(label + " " + name + " current/group-velocity ratio", mode["current"] / mode["density"], mode["k"] / mode["epsilon"])
            values, vectors = np.linalg.eigh(h)
            vector = vectors[:, int(np.argmin(abs(values - mode["epsilon"])))]
            normalized = mode["spinor"] / math.sqrt(mode["density"])
            checks.close(label + " " + name + " independent eigensolver projector", np.outer(normalized, normalized.conj()), np.outer(vector, vector.conj()))
    r, tau, _, _ = closed_coefficients(solution)
    checks.close(label + " independently closed r", solution["r"], r)
    checks.close(label + " independently closed normalized tau", solution["tau"], tau)
    checks.close(label + " conserved R+T", solution["R"] + solution["T"], 1)
    checks.condition(label + " nonnegative bounded R,T", min(solution["R"], solution["T"]) >= 0 and max(solution["R"], solution["T"]) <= 1 + 5e-12)
    for x in (-3.1, -.73, -.01):
        checks.close(f"{label} left current at {x}", current(wave(solution, x, "left")), 1 - solution["R"])
    for x in (.01, .73, 3.1):
        checks.close(f"{label} right current at {x}", current(wave(solution, x, "right")), solution["T"])
    if out["region"] != "evanescent_gap":
        s = scattering_matrix(solution["energy"], mass, solution["right_potential"], solution["left_potential"])
        checks.close(label + " two-port S unitarity", s.conj().T @ s, np.eye(2))
        checks.close(label + " two-port left-incident column", s[:, 0], [solution["r"], solution["tau"]])
    else:
        checks.close(label + " gap density decay", density(wave(solution, 1.1, "right")),
                     density(wave(solution, .1, "right")) * math.exp(-2 * out["k"].imag))


def solution_row(solution, label):
    out = solution["transmitted"]
    return {"case": label, "energy": solution["energy"], "mass": solution["mass"],
            "step_height": solution["right_potential"] - solution["left_potential"], "epsilon_right": out["epsilon"],
            "region": out["region"], "excluded": False, "reason": "", "k_right_real": complex(out["k"]).real,
            "k_right_imag": complex(out["k"]).imag, "right_mode_current": out["current"],
            "right_mode_density": out["density"], "velocity_if_propagating": out["current"] / out["density"] if out["region"] != "evanescent_gap" else None,
            "r_real": float(solution["r"].real), "r_imag": float(solution["r"].imag),
            "tau_real": float(solution["tau"].real), "tau_imag": float(solution["tau"].imag),
            "tau_normalization": out["normalization"], "R": solution["R"], "T": solution["T"],
            "flux_balance_error": abs(solution["R"] + solution["T"] - 1),
            "matching_condition_number": solution["matching_condition_number"]}


class RationalComplex:
    def __init__(self, real=0, imag=0):
        if isinstance(real, RationalComplex):
            self.real, self.imag = real.real, real.imag
        else:
            self.real, self.imag = Fraction(real), Fraction(imag)

    def __add__(self, other):
        other = RationalComplex(other)
        return RationalComplex(self.real + other.real, self.imag + other.imag)

    __radd__ = __add__

    def __neg__(self):
        return RationalComplex(-self.real, -self.imag)

    def __sub__(self, other):
        return self + -RationalComplex(other)

    def __mul__(self, other):
        other = RationalComplex(other)
        return RationalComplex(self.real * other.real - self.imag * other.imag, self.real * other.imag + self.imag * other.real)

    __rmul__ = __mul__

    def conjugate(self):
        return RationalComplex(self.real, -self.imag)

    def __truediv__(self, other):
        other = RationalComplex(other)
        norm = other.real ** 2 + other.imag ** 2
        numerator = self * other.conjugate()
        return RationalComplex(numerator.real / norm, numerator.imag / norm)

    def __eq__(self, other):
        other = RationalComplex(other)
        return self.real == other.real and self.imag == other.imag


def exact_fixtures(checks):
    q, f = RationalComplex, Fraction
    a = q(f(1, 2))
    fixtures = [("positive", q(f(1, 3)), q(f(1, 5)), q(f(6, 5)), f(1, 25), f(24, 25)),
                ("gap", q(0, 1), q(f(-3, 5), f(-4, 5)), q(f(2, 5), f(-4, 5)), f(1), f(0)),
                ("Klein_outgoing", q(3), q(f(-5, 7)), q(f(2, 7)), f(25, 49), f(24, 49)),
                ("Klein_incoming_counterexample", q(-3), q(f(-7, 5)), q(f(-2, 5)), f(49, 25), f(-24, 25))]
    count = 0
    for label, b, expected_r, expected_t, expected_R, expected_T in fixtures:
        r, t = (a - b) / (a + b), 2 * a / (a + b)
        reflection = (r.conjugate() * r).real
        transmission = b.real / a.real * (t.conjugate() * t).real
        identities = [("r", r, expected_r), ("t upper-one chart", t, expected_t),
                      ("continuity upper", q(1) + r, t), ("continuity lower", a * (q(1) - r), b * t),
                      ("R", reflection, expected_R), ("signed T", transmission, expected_T),
                      ("signed current balance", reflection + transmission, f(1))]
        for name, actual, expected in identities:
            checks.condition("exact " + label + " " + name, actual == expected, arithmetic="Fraction complex")
            count += 1
    return count


def fixed_studies(checks):
    threshold_rows, limit_rows, derivative_rows = [], [], []
    for step in (0, 5 / 12, 5 / 3, 35 / 12, 10):
        solution = match_step(5 / 3, 1, step)
        check_solution(checks, solution, "fixed V=" + str(step))
        shifted = match_step(5 / 3 + 128, 1, step + 128, 128)
        checks.close("energy-origin invariance V=" + str(step), [shifted["r"], shifted["tau"]], [solution["r"], solution["tau"]])
    example = match_step(5 / 3, 1, 35 / 12)
    checks.close("owner exact Klein R,T", [example["R"], example["T"]], [25 / 49, 24 / 49])
    checks.close("owner exact outgoing momentum and velocity", [example["transmitted"]["k"], example["transmitted"]["current"] / example["transmitted"]["density"]], [-.75, .6])
    for epsilon in (1.25, 2.0, -1.25, -2.0):
        mode = channel(epsilon, 1)
        k = mode["k"]
        step = 1e-3
        dispersion = lambda p: math.copysign(math.hypot(p, 1), epsilon)
        numerical_velocity = (-dispersion(k + 2 * step) + 8 * dispersion(k + step) - 8 * dispersion(k - step) + dispersion(k - 2 * step)) / (12 * step)
        checks.close(f"independent dispersion derivative eps={epsilon}", numerical_velocity, mode["current"] / mode["density"], tolerance=2e-10)
    for threshold in (-1, 1):
        for side in (-1, 1):
            previous = None
            for delta in (1e-2, 1e-4, 1e-6, 1e-8):
                epsilon = threshold + side * delta
                solution = match_step(2, 1, 2 - epsilon)
                check_solution(checks, solution, f"threshold {threshold} side={side} delta={delta}")
                row = solution_row(solution, "threshold_approach")
                row.update({"threshold_epsilon": threshold, "approach_side": side, "distance": delta})
                threshold_rows.append(row)
                if solution["T"] > 0 and previous is not None:
                    checks.condition(f"threshold transmission decreases {threshold} {delta}", solution["T"] < previous)
                previous = solution["T"]
            if previous:
                checks.condition(f"zero transmission threshold limit {threshold}", previous < .001, final_T=previous)
    for mass, energy, step in ((0, 2, 0), (0, 2, .5), (0, 2, 3), (0, 2, 5)):
        solution = match_step(energy, mass, step)
        check_solution(checks, solution, "massless V=" + str(step))
        checks.close("massless perfect transmission V=" + str(step), [solution["R"], solution["T"]], [0, 1])
        limit_rows.append({"study": "massless_exact", **solution_row(solution, "massless")})
    for mass in (1, .5, .2, .05, .01, 0):
        solution = match_step(2, mass, 5)
        limit_rows.append({"study": "massless_approach", **solution_row(solution, "massless_approach")})
    a = .5
    high_step_limit = 4 * a / (1 + a) ** 2
    previous = float("inf")
    for step in (1e2, 1e4, 1e6):
        solution = match_step(5 / 3, 1, step)
        error = abs(solution["T"] - high_step_limit)
        checks.condition(f"high-step limit V={step}", error < previous, error=error)
        previous = error
        limit_rows.append({"study": "high_step_limit", **solution_row(solution, "high_step"), "limiting_T": high_step_limit, "limit_error": error})
    # Restore c in dimensionless matching via the rest ENERGY M=m*c^2.
    # Both flux-normalization factors acquire the same c, so R,T are unchanged.
    kinetic, potential = .2, .1
    sch_r = (math.sqrt(kinetic) - math.sqrt(kinetic - potential)) / (math.sqrt(kinetic) + math.sqrt(kinetic - potential))
    previous = None
    for speed in (2, 4, 8, 16, 32):
        rest_energy = speed ** 2  # fixed physical m=1
        solution = match_step(rest_energy + kinetic, rest_energy, potential)
        error = abs(solution["R"] - sch_r ** 2)
        order = None if previous is None else math.log2(previous / error)
        limit_rows.append({"study": "nonrelativistic_fixed_kinetic_energy", "physical_mass": 1, "c": speed,
                           "rest_energy": rest_energy, "kinetic_energy": kinetic, "step_height": potential,
                           "R": solution["R"], "T": solution["T"], "Schrodinger_R": sch_r ** 2, "limit_error": error, "observed_order": order})
        forbidden = match_step(rest_energy + kinetic, rest_energy, .3)
        checks.close(f"NR forbidden semi-infinite gap c={speed}", [forbidden["R"], forbidden["T"]], [1, 0])
        previous = error
    checks.condition("NR correction approaches c^-2", abs(order - 2) < .01, observed_order=order)
    # Centered spatial derivative away from the interface: error is O(dx^2).
    # No derivative matching is imposed at the discontinuous potential itself.
    for step in (.5, 2, 3.3):
        solution = match_step(2, 1, step)
        for side, position, potential in (("left", -1.3, 0), ("right", 1.3, step)):
            previous = None
            for dx in (.08, .04, .02, .01):
                psi = wave(solution, position, side)
                derivative = (wave(solution, position + dx, side) - wave(solution, position - dx, side)) / (2 * dx)
                residual = -1j * ALPHA @ derivative + BETA @ psi + (potential - 2) * psi
                error = float(np.linalg.norm(residual) / np.linalg.norm(psi))
                order = None if previous is None else math.log2(previous / error)
                derivative_rows.append({"step_height": step, "region_side": side, "position": position, "dx": dx,
                                        "stationary_equation_residual_per_spinor_norm": error, "observed_order": order})
                previous = error
            checks.condition(f"stationary derivative second order V={step} side={side}", abs(order - 2) < .01, observed_order=order)
    return threshold_rows, limit_rows, derivative_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 parser():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument("--output-dir", type=Path, help="new output files; default fresh temporary directory")
    p.add_argument("--json-output", type=Path, help="override JSON destination")
    p.add_argument("--sweep-csv", type=Path, help="override step-sweep CSV destination")
    p.add_argument("--mass", type=float, default=1.0, help="m=0 or 1e-6..10 (default 1)")
    p.add_argument("--energy", type=float, default=5 / 3, help="incident E>m and <=100 (default 5/3)")
    p.add_argument("--step-max", type=float, default=6.0, help="maximum potential-energy step, 0<Vmax<=1e4 (default 6)")
    p.add_argument("--samples", type=int, default=121, help="uniform step-height samples, 9..2001; thresholds also recorded")
    p.add_argument("--x-span", type=float, default=12.0, help="profile sampling half-span, .1..100 (default 12)")
    p.add_argument("--x-samples", type=int, default=241, help="odd analytic-profile samples, 9..2001 (default 241)")
    return p


def main(argv=None):
    p = parser()
    args = p.parse_args(argv)
    if not math.isfinite(args.mass) or not (args.mass == 0 or 1e-6 <= args.mass <= 10):
        p.error("mass must be zero or finite in 1e-6..10")
    if not math.isfinite(args.energy) or args.energy > 100 or args.energy - args.mass <= 1e-8 * max(1, args.energy, args.mass):
        p.error("energy must exceed mass by a resolved positive gap and be <=100")
    if not math.isfinite(args.step_max) or not 0 < args.step_max <= 1e4:
        p.error("step-max must be finite, positive and <=1e4")
    if not 9 <= args.samples <= 2001 or not 9 <= args.x_samples <= 2001 or args.x_samples % 2 != 1:
        p.error("samples must be 9..2001; x-samples must additionally be odd")
    if not math.isfinite(args.x_span) or not .1 <= args.x_span <= 100:
        p.error("x-span must be finite and in .1..100")
    output = args.output_dir or Path(tempfile.mkdtemp(prefix="klein-step-"))
    paths = {"sweep": args.sweep_csv or output / "sweep.csv", "profiles": output / "profiles.csv",
             "thresholds": output / "thresholds.csv", "limits": output / "limits.csv", "derivatives": output / "derivatives.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()
    exact_count = exact_fixtures(checks)
    thresholds, limits, derivatives = fixed_studies(checks)
    heights = sorted(set(np.linspace(0, args.step_max, args.samples).tolist() +
                         [v for v in (args.energy - args.mass, args.energy + args.mass) if 0 <= v <= args.step_max]))
    sweep = []
    for number, height in enumerate(heights):
        try:
            solution = match_step(args.energy, args.mass, height)
        except ValueError as error:
            sweep.append({"case": "configured_sweep", "energy": args.energy, "mass": args.mass, "step_height": height,
                          "epsilon_right": args.energy - height, "excluded": True, "region": "threshold_excluded",
                          "reason": str(error), "R": None, "T": None})
            continue
        check_solution(checks, solution, f"sweep {number}")
        sweep.append(solution_row(solution, "configured_sweep"))
    profile_steps = [(args.energy - args.mass) / 2, args.energy, args.energy + 1.25 * args.mass] if args.mass else [args.energy / 2, 1.5 * args.energy]
    profiles = []
    for height in profile_steps:
        solution = match_step(args.energy, args.mass, height)
        for x in np.linspace(-args.x_span, args.x_span, args.x_samples):
            side = "left" if x < 0 else "right"
            psi = wave(solution, float(x), side)
            profiles.append({"step_height": height, "region_right": solution["transmitted"]["region"], "x": float(x), "side": side,
                             "upper_real": float(psi[0].real), "upper_imag": float(psi[0].imag),
                             "lower_real": float(psi[1].real), "lower_imag": float(psi[1].imag),
                             "density": density(psi), "probability_current": current(psi)})
    for name, rows in (("sweep", sweep), ("profiles", profiles), ("thresholds", thresholds), ("limits", limits), ("derivatives", derivatives)):
        write_csv(paths[name], rows)
    report = {"program": Path(__file__).name, "passed": checks.passed, "check_count": len(checks.rows),
              "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
              "environment": {"python": platform.python_version(), "numpy": np.__version__, "platform": platform.platform()},
              "randomness": "none", "parameters": {key: getattr(args, key) for key in ("mass", "energy", "step_max", "samples", "x_span", "x_samples")},
              "exact_fraction_checks": exact_count, "excluded_threshold_rows": sum(row["excluded"] for row in sweep), "checks": checks.rows,
              "interpretation": ["Natural units hbar=c=1; V is potential energy; fixed spin at normal incidence.",
                  "Real asymptotic modes have unit absolute Dirac probability current; gap mode has unit interface density.",
                  "The right outgoing negative-energy mode has negative momentum and positive current/group velocity.",
                  "R+T=1 is a positive Dirac probability-current balance, not the scalar signed KG-current convention.",
                  "These stationary one-body coefficients do not count produced pairs or specify a vacuum.",
                  "Threshold rows are explicitly excluded because unit-flux normalization is singular there.",
                  "Profiles are exact stationary solutions sampled on x; no mesh or propagation-step error is implied.",
                  "Finite-difference dx refinement diagnoses the stationary derivative away from the discontinuous interface.",
                  "The separate NR study holds physical mass, kinetic energy and potential fixed as c increases.",
                  "An abrupt semi-infinite step does not model a smooth field, finite barrier or dynamical source."],
              "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())