#!/usr/bin/env python3
"""Relativistic Landau levels: independent spectrum and oscillator-cutoff tests.

Run with --help. Standard library + NumPy only; no network or plotting.
Uniform prescribed B=B*z_hat, B>0, signed q!=0, Phi=0, pi=p-q*A.
H=c*alpha.pi+beta*m*c^2 uses the Dirac basis and minimal coupling only.
There is no anomalous moment, photon emission or radiative energy shift.
One guiding-center orbital is represented; this code does NOT compute the
macroscopic orbital degeneracy |q|B/(2*pi*hbar) per unit transverse area.

Build a truncated oscillator a (states n=0,...,K-1), then form pi_x,pi_y
and diagonalize the complete 4K-dimensional first-order Hamiltonian with
numpy.linalg.eigh. Separately derive its physical invariant Landau blocks:
N=n+(1-sign(q)*s)/2, E_N^2=m^2*c^4+c^2*(pz^2+2*|q|*hbar*B*N).
N=0 has one internal state per NONZERO branch; N>=1 has two.
At m=pz=N=0 the two branches coalesce into a two-dimensional zero space.
No division by E, positive/negative assignment, or duplicate counting is used
there. Analytic normalized spinors check the first-order equation, including
spin mixing; the seed spin label is generally not the full-spinor spin.

Crucial cutoff diagnostic: [a,a_dagger]=I-K*|K-1><K-1|. The full truncated
matrix contains a spurious opposite-spin top-oscillator state with the SAME
energy as the physical lowest level. Blind eigenvalue counting doubles its
multiplicity even as K grows. A second boundary-adjacent physical N=K-1
block is exact in this model but is excluded from reported low-sector
convergence. Cutoff tests compare fixed N<=K-2 across independent matrices.

Degenerate eigensolvers can mix physical and spurious vectors. Consequently
classification uses traces of sector projectors on each complete degenerate
eigenspace, not arbitrary individual eigenvector labels. Raw eigenvalues and
residuals remain exported. This is an oscillator-basis cutoff study, not a
finite spatial box or gauge-coordinate wavefunction simulation.

Exact Fraction fixtures use an unnormalized ladder basis with its factorial
Gram matrix, independently checking ordered blocks and weighted Hermiticity.
Decimal square-root remainders check the massive nonrelativistic limit.
CSV/JSON exports include versions, script hash, CSV hashes and tolerances.
Only new output files are written; defaults use a fresh temporary directory.
Owner: /relativistic-qm/relativistic-landau-levels/.
References: Thaller, The Dirac Equation (1992); Bjorken & Drell,
Relativistic Quantum Mechanics (1964). Tested: Python 3.12.14, NumPy 2.3.5.
"""
from __future__ import annotations

import argparse
import csv
from decimal import Decimal, localcontext
from fractions import Fraction
import hashlib
import json
import math
from pathlib import Path
import platform
import sys
import tempfile

import numpy as np

PAULI = (np.array([[0, 1], [1, 0]], complex), np.array([[0, -1j], [1j, 0]], complex),
         np.diag([1, -1]).astype(complex))
ZERO2 = np.zeros((2, 2), complex)
ALPHA = tuple(np.block([[ZERO2, s], [s, ZERO2]]) for s in PAULI)
BETA = np.diag([1, 1, -1, -1]).astype(complex)
SIGMA_Z = np.diag([1, -1, 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=3e-12):
        absolute = float(np.linalg.norm(np.asarray(actual) - np.asarray(expected)))
        reference = max(1.0, float(np.linalg.norm(expected)))
        self.condition(name, absolute <= tolerance * reference, absolute_error=absolute,
                       scaled_error=absolute / reference, tolerance=tolerance)

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


def oscillator_hamiltonian(cutoff, mass, charge, field, pz, hbar=1.0, speed=1.0):
    """Full finite ladder matrix, independent of the analytic energy formula."""
    a = np.diag(np.sqrt(np.arange(1, cutoff)), 1).astype(complex)
    identity = np.eye(cutoff)
    magnitude = abs(charge) * hbar * field
    sign = 1 if charge > 0 else -1
    px = math.sqrt(magnitude / 2) * (a + a.conj().T)
    py = -1j * sign * math.sqrt(magnitude / 2) * (a - a.conj().T)
    momenta = (px, py, pz * identity)
    h = mass * speed ** 2 * np.kron(BETA, identity)
    for alpha, momentum in zip(ALPHA, momenta):
        h += speed * np.kron(alpha, momentum)
    return h, a, momenta


def sector_indices(cutoff, charge, level):
    aligned = 0 if charge > 0 else 1
    seeds = [(aligned, 0)] if level == 0 else [(aligned, level), (1 - aligned, level - 1)]
    return [(spin + lower) * cutoff + n for lower in (0, 2) for spin, n in seeds], seeds


def analytic_block(level, mass, charge, field, pz, hbar, speed):
    sign = 1 if charge > 0 else -1
    transverse = math.sqrt(2 * abs(charge) * hbar * field * level)
    t = np.array([[sign * pz]], complex) if level == 0 else np.array([[sign * pz, transverse], [transverse, -sign * pz]], complex)
    identity = np.eye(len(t))
    rest = mass * speed ** 2
    block = np.block([[rest * identity, speed * t], [speed * t, -rest * identity]])
    # Preserve tiny nonzero gaps so the CLI can reject unresolved branches;
    # squaring first could underflow and falsely label a nonzero gap as apex.
    magnitude = math.hypot(rest, speed * pz, speed * transverse)
    return block, t, magnitude


def analytic_spinors(t, magnitude, mass, speed):
    if magnitude == 0:
        # One coalesced zero eigenspace, not two independently labelled branches.
        return [("zero", index, np.eye(2, dtype=complex)[:, index]) for index in range(2)]
    rest = mass * speed ** 2
    normalization = math.sqrt((magnitude + rest) / (2 * magnitude))
    columns = []
    for index, seed in enumerate(np.eye(len(t), dtype=complex)):
        lower = speed * (t @ seed) / (magnitude + rest)
        columns.append(("positive", index, normalization * np.concatenate((seed, lower))))
        columns.append(("negative", index, normalization * np.concatenate((-lower, seed))))
    return columns


def eigen_clusters(values, tolerance):
    clusters = []
    for index, value in enumerate(values):
        # Compare with the first member, avoiding tolerance-chain merging.
        if not clusters or abs(value - values[clusters[-1][0]]) > tolerance:
            clusters.append([index])
        else:
            clusters[-1].append(index)
    return clusters


def full_spectrum_study(checks, label, cutoff, mass, charge, field, pz, hbar=1.0, speed=1.0, max_level=3):
    h, a, momenta = oscillator_hamiltonian(cutoff, mass, charge, field, pz, hbar, speed)
    dimension = len(h)
    identity = np.eye(cutoff)
    top = np.zeros((cutoff, cutoff))
    top[-1, -1] = 1
    comm = a @ a.conj().T - a.conj().T @ a
    checks.close(label + " ladder commutator with boundary", comm, identity - cutoff * top)
    checks.close(label + " trace ladder commutator zero", np.trace(comm), 0)
    checks.condition(label + " finite ladder cannot obey CCR globally", np.linalg.norm(comm - identity) > cutoff - 0.1)
    checks.close(label + " H Hermitian", h.conj().T, h)
    pi2 = sum(p @ p for p in momenta)
    magnetic = charge * hbar * field * speed ** 2
    naive_square = mass ** 2 * speed ** 4 * np.eye(dimension) + speed ** 2 * np.kron(np.eye(4), pi2) - magnetic * np.kron(SIGMA_Z, identity)
    boundary_correction = magnetic * cutoff * np.kron(SIGMA_Z, top)
    checks.close(label + " H squared with cutoff correction", h @ h, naive_square + boundary_correction)
    checks.close(label + " naive square error only at top", h @ h - naive_square, boundary_correction)
    values, vectors = np.linalg.eigh(h)
    checks.close(label + " full independent eigen residual", h @ vectors, vectors * values)
    checks.close(label + " full eigenbasis orthonormal", vectors.conj().T @ vectors, np.eye(dimension))
    # Form a complete invariant partition, including the two artifact states.
    physical_ids = [sector_indices(cutoff, charge, n)[0] for n in range(cutoff)]
    opposite = 1 if charge > 0 else 0
    artifact_ids = [opposite * cutoff + cutoff - 1, (opposite + 2) * cutoff + cutoff - 1]
    all_ids = [i for ids in physical_ids for i in ids] + artifact_ids
    checks.condition(label + " sector partition complete and disjoint", sorted(all_ids) == list(range(dimension)))
    spectral_scale = max(1.0, float(np.max(np.abs(values))))
    tolerance = 2e-11 * spectral_scale
    groups = eigen_clusters(values, tolerance)
    cluster_rows, raw_rows, level_rows, spinor_rows = [], [], [], []
    for group_number, group in enumerate(groups):
        columns = vectors[:, group]
        magnitude = float(np.mean(values[group]))
        weights = [float(np.sum(np.abs(columns[ids, :]) ** 2)) for ids in physical_ids]
        artifact_weight = float(np.sum(np.abs(columns[artifact_ids, :]) ** 2))
        expected_weights = []
        for n in range(cutoff):
            _, _, e = analytic_block(n, mass, charge, field, pz, hbar, speed)
            multiplicity = 1 if n == 0 else 2
            expected = 2 * multiplicity if e == 0 and abs(magnitude) < tolerance else multiplicity if min(abs(magnitude - e), abs(magnitude + e)) < tolerance else 0
            expected_weights.append(expected)
        _, _, e0 = analytic_block(0, mass, charge, field, pz, hbar, speed)
        expected_artifact = 2 if e0 == 0 and abs(magnitude) < tolerance else 1 if min(abs(magnitude - e0), abs(magnitude + e0)) < tolerance else 0
        checks.close(f"{label} cluster {group_number} invariant physical multiplicities", weights, expected_weights, tolerance=2e-10)
        checks.close(f"{label} cluster {group_number} artifact multiplicity", artifact_weight, expected_artifact, tolerance=2e-10)
        checks.close(f"{label} cluster {group_number} partition trace", sum(weights) + artifact_weight, len(group))
        # A complex unitary rotation inside the cluster changes eigenvectors,
        # but cannot change its physical/artifact projector trace.
        indices = np.arange(len(group))
        rotation = np.exp(2j * np.pi * np.outer(indices, indices) / len(group)) / math.sqrt(len(group))
        rotated = columns @ rotation
        rotated_weights = [float(np.sum(np.abs(rotated[ids, :]) ** 2)) for ids in physical_ids]
        rotated_artifact = float(np.sum(np.abs(rotated[artifact_ids, :]) ** 2))
        checks.close(f"{label} cluster {group_number} basis-independent classification",
                     rotated_weights + [rotated_artifact], weights + [artifact_weight])
        cluster_rows.append({"case": label, "cutoff": cutoff, "cluster": group_number, "energy_mean": magnitude,
                             "energy_min": float(values[group[0]]), "energy_max": float(values[group[-1]]),
                             "raw_multiplicity": len(group), "physical_N0_weight": weights[0],
                             "physical_interior_weight": sum(weights[:-1]), "boundary_adjacent_NKminus1_weight": weights[-1],
                             "spurious_top_weight": artifact_weight,
                             "physical_level_weights": json.dumps(weights, separators=(",", ":"))})
        for index in group:
            residual = float(np.linalg.norm(h @ vectors[:, index] - values[index] * vectors[:, index]))
            raw_rows.append({"case": label, "cutoff": cutoff, "eigenvalue_index": index, "cluster": group_number,
                             "energy": float(values[index]), "absolute_eigen_residual": residual,
                             "classification": "use complete cluster projector traces; individual eigenvectors can mix"})
    checks.close(label + " total artifact dimensions", sum(row["spurious_top_weight"] for row in cluster_rows), 2)
    for level in range(max_level + 1):
        ids, seeds = sector_indices(cutoff, charge, level)
        actual = h[np.ix_(ids, ids)]
        block, t, magnitude = analytic_block(level, mass, charge, field, pz, hbar, speed)
        complement = sorted(set(range(dimension)) - set(ids))
        checks.close(f"{label} N={level} sector invariant", h[np.ix_(complement, ids)], np.zeros((len(complement), len(ids))))
        checks.close(f"{label} N={level} first-order block", actual, block)
        checks.close(f"{label} N={level} squared energy", actual @ actual, magnitude ** 2 * np.eye(len(ids)))
        block_eigenvalues = np.linalg.eigvalsh(actual)
        multiplicity = len(seeds)
        expected_values = np.array([-magnitude] * multiplicity + [magnitude] * multiplicity)
        checks.close(f"{label} N={level} independent block spectrum", block_eigenvalues, expected_values)
        sigmaz = np.diag([1 if (index // cutoff) % 2 == 0 else -1 for index in ids])
        spinors = analytic_spinors(t, magnitude, mass, speed)
        columns = np.column_stack([v for _, _, v in spinors])
        checks.close(f"{label} N={level} analytic modes orthonormal and complete", columns.conj().T @ columns, np.eye(len(ids)))
        if magnitude:
            for branch in (-1, 1):
                projector = (np.eye(len(ids)) + branch * actual / magnitude) / 2
                checks.close(f"{label} N={level} projector {branch}", projector @ projector, projector)
                checks.close(f"{label} N={level} projector trace {branch}", np.trace(projector), multiplicity)
        else:
            checks.condition(label + " zero apex has one coalesced 2D space", len(spinors) == 2 and all(branch == "zero" for branch, _, _ in spinors))
        for branch, seed, vector in spinors:
            signed_energy = magnitude if branch == "positive" else -magnitude if branch == "negative" else 0.0
            residual = float(np.linalg.norm(h[:, ids] @ vector - signed_energy * np.eye(dimension)[:, ids] @ vector))
            checks.close(f"{label} N={level} {branch} seed={seed} full first-order spinor", residual, 0)
            spin_mean = float(np.vdot(vector, sigmaz @ vector).real)
            variance = max(0.0, 1 - spin_mean ** 2)
            if level == 0:
                checks.close(f"{label} LLL full-spin label {branch} {seed}", spin_mean, 1 if charge > 0 else -1)
            level_rows.append({"case": label, "cutoff": cutoff, "N": level, "branch": branch, "seed": seed,
                               "energy": signed_energy, "internal_multiplicity_per_nonzero_branch": multiplicity if magnitude else None,
                               "zero_space_dimension": 2 if not magnitude else None,
                               "spin_z_mean_in_hbar_over2_units": spin_mean, "spin_z_variance_same_units_squared": variance,
                               "first_order_residual": residual, "below_boundary_adjacent_sector": level <= cutoff - 2})
            for index, value in zip(ids, vector):
                spinor_rows.append({"case": label, "N": level, "branch": branch, "seed": seed,
                                    "dirac_component_1_based": index // cutoff + 1, "oscillator_n": index % cutoff,
                                    "coefficient_real": float(value.real), "coefficient_imag": float(value.imag)})
        for spin, n in seeds:
            sign_spin = 1 if spin == 0 else -1
            mapped = n + (1 - (1 if charge > 0 else -1) * sign_spin) // 2
            checks.condition(f"{label} seed n={n},s={sign_spin} maps to N={level}", mapped == level)
        if level > 0:
            checks.condition(f"{label} N={level} seed spin is not conserved full spin", np.linalg.norm(actual @ sigmaz - sigmaz @ actual) > 1e-6)
    return {"clusters": cluster_rows, "raw": raw_rows, "levels": level_rows, "spinors": spinor_rows}


def exact_fixtures(checks):
    """Unnormalized basis has a|n>=n|n-1>, a_dagger|n>=|n+1>.

    Divide the factorial Gram matrix in N>=1 by (N-1)!: G_seed=diag(N,1).
    This exact rational fixture is not tested with the wrong Euclidean norm.
    """
    f = Fraction
    def eye(n):
        return [[f(int(i == j)) for j in range(n)] for i in range(n)]
    def mul(a, b):
        return [[sum((a[i][k] * b[k][j] for k in range(len(b))), f()) for j in range(len(b[0]))] for i in range(len(a))]
    def scale(a, x):
        return [[x * v for v in row] for row in a]
    def transpose(a):
        return [list(row) for row in zip(*a)]
    def block(a, b, c, d):
        return [ar + br for ar, br in zip(a, b)] + [cr + dr for cr, dr in zip(c, d)]
    count = 0
    for charge_sign in (-1, 1):
        for kappa in (f(1), f(3, 2)):
            for n in range(5):
                mass, speed, pz = f(3, 4), f(7, 5), f(5, 3)
                t = [[charge_sign * pz]] if n == 0 else [[charge_sign * pz, kappa], [n * kappa, -charge_sign * pz]]
                identity = eye(len(t))
                zero = scale(identity, 0)
                gram = [[f(1)]] if n == 0 else [[f(n), f()], [f(), f(1)]]
                h = block(scale(identity, mass * speed ** 2), scale(t, speed), scale(t, speed), scale(identity, -mass * speed ** 2))
                full_gram = block(gram, zero, zero, gram)
                e2 = mass ** 2 * speed ** 4 + speed ** 2 * (pz ** 2 + n * kappa ** 2)
                fixtures = [("T squared", mul(t, t), scale(identity, pz ** 2 + n * kappa ** 2)),
                            ("H squared", mul(h, h), scale(eye(len(h)), e2)),
                            ("weighted Hermiticity", mul(transpose(h), full_gram), mul(full_gram, h))]
                for name, actual, expected in fixtures:
                    checks.condition(f"exact sign={charge_sign} kappa={kappa} N={n} {name}", actual == expected, arithmetic="Fraction")
                    count += 1
    return count


def cutoff_study(checks):
    rows = []
    spectra = []
    for cutoff in (6, 10, 14):
        result = full_spectrum_study(checks, f"cutoff_{cutoff}", cutoff, 1.0, -1.0, 0.7, 0.4, max_level=3)
        for level in range(4):
            states = [r for r in result["levels"] if r["N"] == level]
            h, _, _ = oscillator_hamiltonian(cutoff, 1.0, -1.0, 0.7, 0.4)
            ids, _ = sector_indices(cutoff, -1.0, level)
            numerical = np.linalg.eigvalsh(h[np.ix_(ids, ids)])
            exact_energy = math.sqrt(1 + 0.4 ** 2 + 1.4 * level)
            error = float(np.max(np.abs(np.abs(numerical) - exact_energy)))
            rows.append({"cutoff": cutoff, "N": level, "maximum_energy_error": error,
                         "low_sector_dimension": len(states), "spurious_total_dimension": sum(r["spurious_top_weight"] for r in result["clusters"]),
                         "boundary_adjacent_N": cutoff - 1, "comparison": "fixed physical N across independently built matrices"})
            spectra.append((cutoff, level, numerical))
    for level in range(4):
        selected = [v for _, n, v in spectra if n == level]
        for index in (1, 2):
            checks.close(f"cutoff-independent physical N={level} comparison={index}", selected[index], selected[0])
    return rows


def nonrelativistic_study(checks):
    rows = []
    with localcontext() as context:
        context.prec = 90
        mass, charge_abs, field, pz, hbar = map(Decimal, ("1", "1", ".7", ".4", "1"))
        for n in range(5):
            value = pz ** 2 + 2 * charge_abs * field * hbar * n
            previous = None
            for speed in map(Decimal, ("2", "4", "8", "16", "32")):
                energy = (mass ** 2 * speed ** 4 + speed ** 2 * value).sqrt()
                rest_removed = energy - mass * speed ** 2
                pauli = value / (2 * mass)
                corrected = pauli - value ** 2 / (8 * mass ** 3 * speed ** 2)
                error = rest_removed - corrected
                r = value / (mass ** 2 * speed ** 2)
                upper = value ** 3 / (16 * mass ** 5 * speed ** 4)
                lower = upper / ((1 + r) ** 2 * (1 + r).sqrt())
                checks.condition(f"NR N={n} c={speed} exact Taylor bounds", lower <= error <= upper, arithmetic="Decimal 90 digits")
                pauli_error = pauli - rest_removed
                order = None if previous is None else math.log2(float(previous / error))
                rows.append({"study": "fixed_massive_NR_fixture", "mass": str(mass), "N": n, "c": str(speed),
                             "F_eigenvalue": str(value), "r_F_over_m2c2": str(r), "exact_rest_removed_energy": str(rest_removed),
                             "Pauli_energy": str(pauli), "FW_corrected_energy": str(corrected),
                             "Pauli_error": str(pauli_error), "FW_remainder": str(error), "lower_bound": str(lower),
                             "upper_bound": str(upper), "FW_observed_order": order})
                previous = error
            checks.condition(f"NR N={n} approaches fourth order", abs(rows[-1]["FW_observed_order"] - 4) < .04,
                             observed_order=rows[-1]["FW_observed_order"])
    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 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 summary JSON destination")
    p.add_argument("--spectrum-csv", type=Path, help="override full raw eigenvalue CSV destination")
    p.add_argument("--cutoff", type=int, default=12, help="oscillator dimension K, 4..64 (default 12)")
    p.add_argument("--max-level", type=int, default=4, help="highest exported low physical N, 0..K-2 (default 4)")
    p.add_argument("--mass", type=float, default=1.0, help="m, 0..5 (default 1); m=0 supported")
    p.add_argument("--charge", type=float, default=-1.0, help="signed q, 0.1<=|q|<=3 (default -1)")
    p.add_argument("--field", type=float, default=.7, help="positive B, 0.05..5 (default .7)")
    p.add_argument("--pz", type=float, default=.4, help="longitudinal momentum, -5..5 (default .4)")
    p.add_argument("--hbar", type=float, default=1.0, help="hbar, 0.5..2 (default 1)")
    p.add_argument("--c", type=float, default=1.0, help="c, 0.5..5 (default 1)")
    return p


def main(argv=None):
    p = parser()
    args = p.parse_args(argv)
    bounds = {"mass": (0, 5), "field": (.05, 5), "pz": (-5, 5), "hbar": (.5, 2), "c": (.5, 5)}
    for name, (low, high) in bounds.items():
        value = getattr(args, name)
        if not math.isfinite(value) or not low <= value <= high:
            p.error(f"{name} must be finite and between {low} and {high}")
    if not math.isfinite(args.charge) or not .1 <= abs(args.charge) <= 3:
        p.error("charge must be finite and satisfy 0.1<=abs(charge)<=3")
    if not 4 <= args.cutoff <= 64 or not 0 <= args.max_level <= args.cutoff - 2:
        p.error("cutoff must be 4..64 and max-level must be 0..cutoff-2")
    magnitudes = [analytic_block(n, args.mass, args.charge, args.field, args.pz, args.hbar, args.c)[2]
                  for n in range(args.cutoff)]
    if (args.mass != 0 or args.pz != 0) and magnitudes[0] == 0:
        p.error("nonzero gap underflowed during scaling; only exact mass=pz=0 is the zero apex")
    cluster_tolerance = 2e-11 * max(1.0, magnitudes[-1])
    distinct_gaps = [b - a for a, b in zip(magnitudes, magnitudes[1:])]
    if magnitudes[0] > 0:
        distinct_gaps.append(2 * magnitudes[0])
    if min(distinct_gaps) <= 4 * cluster_tolerance:
        p.error("distinct energy gaps are unresolved at the declared clustering tolerance; use exact m=pz=0 for the apex")
    output = args.output_dir or Path(tempfile.mkdtemp(prefix="landau-levels-"))
    paths = {"raw": args.spectrum_csv or output / "spectrum.csv", "clusters": output / "clusters.csv",
             "levels": output / "levels.csv", "spinors": output / "spinors.csv", "cutoff": output / "cutoff.csv",
             "nonrelativistic": output / "nonrelativistic.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)
    result = full_spectrum_study(checks, "configured", args.cutoff, args.mass, args.charge, args.field, args.pz,
                                 args.hbar, args.c, args.max_level)
    for charge in (-1.0, 1.0):
        for mass in (0.0, 1.0):
            for pz in (0.0, .4):
                full_spectrum_study(checks, f"fixed_q={charge}_m={mass}_pz={pz}", 7, mass, charge, .7, pz, max_level=3)
    cutoff = cutoff_study(checks)
    nonrelativistic = nonrelativistic_study(checks)
    for name in ("raw", "clusters", "levels", "spinors"):
        write_csv(paths[name], result[name])
    write_csv(paths["cutoff"], cutoff)
    write_csv(paths["nonrelativistic"], nonrelativistic)
    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 ("cutoff", "max_level", "mass", "charge", "field", "pz", "hbar", "c")},
              "exact_fraction_checks": exact_count, "checks": checks.rows,
              "interpretation": ["One guiding-center label only; internal multiplicities do not include orbital degeneracy.",
                  "Finite ladder CCR defect generates two spurious top-boundary Dirac states, degenerate with the true LLL.",
                  "Individual degenerate eigenvectors are not classified; complete eigenspace projector traces are basis independent.",
                  "The N=K-1 block is boundary-adjacent and excluded from the low-sector cutoff comparison.",
                  "Massless pz=0 LLL is a coalesced two-dimensional zero space, not distinct nonzero branches.",
                  "A Pauli seed spin label is generally not a full-Dirac-spinor Sigma_z eigenvalue.",
                  "Numerical parameter ranges keep tested gaps resolved; no arbitrary-parameter precision guarantee is made.",
                  "Exact rational tests use the factorial Gram metric in the unnormalized oscillator basis.",
                  "The NR table is a separate fixed massive fixture even when the configured output has mass zero.",
                  "No coordinate wavefunctions, finite spatial-box convergence, anomalous moment or dynamical photons are simulated."],
              "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())