"""Bounded finite-emitter frequency-domain qualification, not a PF-001 result.

Natural units set wave speed and angular frequency to one. The prescribed
source is two normalized Gaussian scalar emitters. Source work and reaction
come from the outgoing Green function; far-field flux comes from a separate
angular Fourier sum over the discretized emitters.
"""

from __future__ import annotations

import argparse
import json
import math
from itertools import product

WAVE_NUMBER = 1.0
SEPARATION = math.pi / 2
WIDTH = 0.25
QUALIFIED_WIDTHS = (0.2, 0.25, 0.3)
FOUR_PI = 4 * math.pi


def cloud(
    nodes_per_axis: int, center_z: float, width: float = WIDTH, edge_rule: str = "full",
) -> list[tuple[float, float, float, float]]:
    if nodes_per_axis not in (3, 5, 7, 9):
        raise ValueError("nodes_per_axis must be one of 3, 5, 7, 9")
    if width not in QUALIFIED_WIDTHS:
        raise ValueError("source width is not a qualified setting")
    if edge_rule not in ("full", "half"):
        raise ValueError("unregistered source endpoint rule")
    offsets = [width * 3 * (2 * i / (nodes_per_axis - 1) - 1) for i in range(nodes_per_axis)]
    weights = [math.exp(-0.5 * (offset / width) ** 2) for offset in offsets]
    if edge_rule == "half":
        weights[0] *= 0.5
        weights[-1] *= 0.5
    total = sum(weights) ** 3
    return [
        (x, y, center_z + z, wx * wy * wz / total)
        for (x, wx), (y, wy), (z, wz) in product(zip(offsets, weights), repeat=3)
    ]


def source_kernels(
    left: list[tuple[float, float, float, float]],
    right: list[tuple[float, float, float, float]],
) -> tuple[float, float, float]:
    """Return imaginary self Green kernel, cross kernel, and z-gradient.

    The real coincident Green singularity does not enter cycle-mean source
    work. Its imaginary limit is k/(4 pi). Self-force pairs cancel exactly.
    """
    self_imag = 0.0
    for i, (x, y, z, weight) in enumerate(left):
        self_imag += weight * weight * WAVE_NUMBER / FOUR_PI
        for xx, yy, zz, other_weight in left[i + 1 :]:
            radius = math.dist((x, y, z), (xx, yy, zz))
            self_imag += 2 * weight * other_weight * math.sin(radius) / (FOUR_PI * radius)

    cross_imag = 0.0
    gradient_imag = 0.0
    for x, y, z, weight in left:
        for xx, yy, zz, other_weight in right:
            radius = math.dist((x, y, z), (xx, yy, zz))
            pair_weight = weight * other_weight
            cross_imag += pair_weight * math.sin(radius) / (FOUR_PI * radius)
            gradient_imag += (
                pair_weight
                * (z - zz)
                * (radius * math.cos(radius) - math.sin(radius))
                / (FOUR_PI * radius**3)
            )
    return self_imag, cross_imag, gradient_imag


def angular_flux(
    left: list[tuple[float, float, float, float]],
    right: list[tuple[float, float, float, float]],
    phase: float,
    channel_overlap: float,
    polar_samples: int,
    azimuth_samples: int,
) -> tuple[float, float]:
    if not (24 <= polar_samples <= 96 and 8 <= azimuth_samples <= 48):
        raise ValueError("angular sample counts exceed the qualified bounds")
    if not -1 <= channel_overlap <= 1:
        raise ValueError("channel overlap must be within [-1, 1]")
    power = 0.0
    momentum_z = 0.0
    for polar in range(polar_samples):
        u = -1 + (polar + 0.5) * 2 / polar_samples
        transverse = math.sqrt(1 - u * u)
        for azimuth in range(azimuth_samples):
            angle = 2 * math.pi * (azimuth + 0.5) / azimuth_samples
            nx = transverse * math.cos(angle)
            ny = transverse * math.sin(angle)
            left_amplitude = sum(
                weight * complex(math.cos(theta), -math.sin(theta))
                for x, y, z, weight in left
                for theta in [nx * x + ny * y + u * z]
            )
            right_amplitude = sum(
                weight * complex(math.cos(theta), -math.sin(theta))
                for x, y, z, weight in right
                for theta in [nx * x + ny * y + u * z]
            ) * complex(math.cos(phase), math.sin(phase))
            intensity = (
                abs(left_amplitude) ** 2
                + abs(right_amplitude) ** 2
                + 2 * channel_overlap * (left_amplitude.conjugate() * right_amplitude).real
            )
            power += intensity
            momentum_z += u * intensity
    solid_angle = 4 * math.pi / (polar_samples * azimuth_samples)
    coefficient = solid_angle / (32 * math.pi**2)
    return coefficient * power, -coefficient * momentum_z


def analytic(phase: float, channel_overlap: float, width: float = WIDTH) -> tuple[float, float]:
    """Continuum Gaussian far-field result at fixed source amplitude."""
    if width not in QUALIFIED_WIDTHS:
        raise ValueError("source width is not a qualified setting")
    form_factor_squared = math.exp(-width**2)
    sinc = math.sin(SEPARATION) / SEPARATION
    shape = (math.sin(SEPARATION) - SEPARATION * math.cos(SEPARATION)) / SEPARATION**2
    return (
        form_factor_squared / FOUR_PI * (1 + channel_overlap * math.cos(phase) * sinc),
        -form_factor_squared / FOUR_PI * channel_overlap * math.sin(phase) * shape,
    )


def evaluate(nodes_per_axis: int, polar_samples: int, azimuth_samples: int) -> dict:
    left = cloud(nodes_per_axis, -SEPARATION / 2)
    right = cloud(nodes_per_axis, SEPARATION / 2)
    self_imag, cross_imag, gradient_imag = source_kernels(left, right)
    cases = (
        ("forward", -math.pi / 2, 1.0),
        ("reverse", math.pi / 2, 1.0),
        ("in_phase", 0.0, 1.0),
        ("opposed", math.pi, 1.0),
        ("orthogonal", -math.pi / 2, 0.0),
    )
    findings = {}
    for name, phase, overlap in cases:
        flux_power, flux_force = angular_flux(left, right, phase, overlap, polar_samples, azimuth_samples)
        source_power = self_imag + overlap * math.cos(phase) * cross_imag
        source_force = -overlap * math.sin(phase) * gradient_imag
        expected_power, expected_force = analytic(phase, overlap)
        findings[name] = {
            "source_work": source_power,
            "radiated_power": flux_power,
            "source_reaction_z": source_force,
            "radiation_recoil_z": flux_force,
            "continuum_power": expected_power,
            "continuum_recoil_z": expected_force,
            "power_residual": source_power - flux_power,
            "momentum_residual": source_force - flux_force,
        }
    return {
        "status": "qualification only; not a canonical result or physical validation",
        "model": "massless linear scalar field with two prescribed finite Gaussian sources",
        "nodes_per_axis": nodes_per_axis,
        "polar_samples": polar_samples,
        "azimuth_samples": azimuth_samples,
        "cases": findings,
    }


def assess(run: dict) -> list[str]:
    """Return failed qualification checks; never assigns a canonical state."""
    cases = run.get("cases", {})
    required = {"forward", "reverse", "in_phase", "opposed", "orthogonal"}
    if set(cases) != required:
        return ["incomplete qualification cohort"]
    failures = []
    required_fields = {
        "source_work", "radiated_power", "source_reaction_z", "radiation_recoil_z",
        "continuum_power", "continuum_recoil_z", "power_residual", "momentum_residual",
    }
    for name, case in cases.items():
        if not isinstance(case, dict) or set(case) != required_fields:
            failures.append(f"{name}: incomplete or extra measurements")
            continue
        if any(not isinstance(value, (int, float)) or not math.isfinite(value) for value in case.values()):
            failures.append(f"{name}: non-finite measurement")
            continue
        if abs(case["power_residual"] - (case["source_work"] - case["radiated_power"])) > 1e-12:
            failures.append(f"{name}: power residual is inconsistent")
        if abs(case["momentum_residual"] - (case["source_reaction_z"] - case["radiation_recoil_z"])) > 1e-12:
            failures.append(f"{name}: momentum residual is inconsistent")
        if abs(case["source_work"] - case["radiated_power"]) > 1e-4:
            failures.append(f"{name}: source work and radiation disagree")
        if abs(case["source_reaction_z"] - case["radiation_recoil_z"]) > 1e-4:
            failures.append(f"{name}: source reaction and recoil disagree")
        if abs(case["radiated_power"] - case["continuum_power"]) > 1e-3:
            failures.append(f"{name}: finite-source power differs from continuum control")
        if abs(case["radiation_recoil_z"] - case["continuum_recoil_z"]) > 1e-3:
            failures.append(f"{name}: finite-source recoil differs from continuum control")
    if all(isinstance(case, dict) and isinstance(case.get("source_reaction_z"), (int, float)) for case in cases.values()):
        if cases["forward"]["source_reaction_z"] <= 0 or cases["reverse"]["source_reaction_z"] >= 0:
            failures.append("phase reversal did not reverse source reaction")
        for name in ("in_phase", "opposed", "orthogonal"):
            if abs(cases[name]["source_reaction_z"]) > 1e-8:
                failures.append(f"{name}: null control produced source reaction")
    return failures


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--nodes", type=int, choices=(3, 5, 7, 9), default=5)
    parser.add_argument("--polar", type=int, default=48)
    parser.add_argument("--azimuth", type=int, default=16)
    args = parser.parse_args()
    run = evaluate(args.nodes, args.polar, args.azimuth)
    run["qualification_failures"] = assess(run)
    print(json.dumps(run, indent=2, allow_nan=False))
    if run["qualification_failures"]:
        raise SystemExit(1)


if __name__ == "__main__":
    main()