"""Where the squares sit: a check of the wound-chain reading.

The lab's quantum mode uses three hopping rules and never sees the thread's
shape. This script is the one place on the findings page where geometry
enters the equations. It places the 108 beads exactly where the lab draws
them on the torus preset (ring radius 1, tube radius r), then asks two
questions for a range of tube radii:

1. Which beads are nearest each other in space? The squares (k, k+27, k+54,
   k+81) are the four corners of the tube at one station; k+28 is the same
   corner at the next station.
2. If every pair of beads coupled with strength falling off as 1/d^3 (a
   dipolar law) below a cutoff of 2.5 bead spacings, would the square qubit
   and perfect state transfer survive?

This is one simple model, chosen by us. It is a computation, not an
established result. Requires Python 3 and NumPy. Prints JSON.
Written 2026-09-23 for the findings page, section 9. Extended the same day
(patch 03): tube radii 0.08 and 0.02, and the one-corner transfer under a
1/d^6 (Rydberg-type) law with the same cutoff and time window.
"""

import json
import math

import numpy as np

P, Q, NB = 4, 27, 108


def core(u, R):
    return np.array([R * math.cos(u), R * math.sin(u), 0.0])


def surf(u, v, R, r):
    # Same frame as engine.ts surf(): tangent by central difference, normal
    # pointing out from the ring's axis, binormal completing the frame.
    c, a, b = core(u, R), core(u + 1e-3, R), core(u - 1e-3, R)
    t = (a - b) / np.linalg.norm(a - b)
    n = np.array([math.cos(u), math.sin(u), 0.0])
    n = n - n.dot(t) * t
    n /= np.linalg.norm(n)
    bn = np.cross(t, n)
    return c + n * r * math.cos(v) + bn * r * math.sin(v)


def beads(r, R=1.0):
    return np.array([surf(P * 2 * math.pi * k / NB, Q * 2 * math.pi * k / NB, R, r) for k in range(NB)])


def distances(r):
    x = beads(r)
    return np.linalg.norm(x[:, None] - x[None], axis=2)


def step_distances(d, s):
    return np.array([d[k, (k + s) % NB] for k in range(NB)])


def nearest_is_square(d):
    """Share of beads whose nearest other bead is k +/- 27."""
    hits = 0
    for k in range(NB):
        j = int(np.argsort(d[k])[1])
        s = (j - k) % NB
        hits += min(s, NB - s) == 27
    return hits / NB


def single_corner(r, alpha, cutoff=2.5, t_max=40.0, samples=8001):
    """Best P(bead 55) from a bead-1 start, J = (d27/d)^alpha below the cutoff."""
    d = distances(r)
    spacing = step_distances(d, 1).mean()
    safe = np.where(d > 0, d, 1.0)
    j = np.where((d > 0) & (d < cutoff * spacing), (d[0, 27] / safe) ** alpha, 0.0)
    w, v = np.linalg.eigh(-j)
    ts = np.linspace(0, t_max, samples)
    amp55 = (np.exp(-1j * np.outer(ts, w)) * v[0]) @ v[54]
    return float((np.abs(amp55) ** 2).max())


def dipolar_check(r, cutoff=2.5, t_max=40.0, samples=8001):
    d = distances(r)
    spacing = step_distances(d, 1).mean()
    d27 = d[0, 27]
    safe = np.where(d > 0, d, 1.0)
    j = np.where((d > 0) & (d < cutoff * spacing), (d27 / safe) ** 3, 0.0)  # square edge = 1
    w, v = np.linalg.eigh(-j)
    ts = np.linspace(0, t_max, samples)
    phases = np.exp(-1j * np.outer(ts, w))  # samples x NB

    corner = np.zeros(NB)
    corner[0] = 1
    amp55 = (phases * (v.T @ corner)) @ v[54]
    single = np.abs(amp55) ** 2

    qubit = np.zeros(NB)
    qubit[0] = qubit[54] = 1 / math.sqrt(2)
    cq = phases * (v.T @ qubit)
    empty = np.abs(cq @ v[27]) ** 2 + np.abs(cq @ v[81]) ** 2

    return {
        "coupling_k28_over_k27": round(float(j[0, 28]), 3),
        "coupling_k54_over_k27": round(float(j[0, 54]), 3),
        "best_single_corner_to_bead_55": round(float(single.max()), 3),
        "best_square_qubit_P_empty": round(float(empty.max()), 3),
    }


report = {"model": "lab torus preset, ring radius 1; J = (d27/d)^3 for d < 2.5 bead spacings; t from 0 to 40"}
for r in (0.4, 0.25, 0.15, 0.1, 0.08, 0.05, 0.02):
    d = distances(r)
    row = {
        f"distance_step_{s}": round(float(step_distances(d, s).mean()), 3) for s in (1, 26, 27, 28, 54)
    }
    row["distance_step_28_min"] = round(float(step_distances(d, 28).min()), 3)
    row["share_nearest_is_square"] = round(nearest_is_square(d), 3)
    row.update(dipolar_check(r))
    row["best_single_corner_to_bead_55_d6"] = round(single_corner(r, 6), 4)
    report[f"tube_{r}"] = row

# The tube radius below which every bead's nearest neighbour is its square.
lo, hi = 0.05, 0.4
for _ in range(40):
    mid = (lo + hi) / 2
    if nearest_is_square(distances(mid)) == 1.0:
        lo = mid
    else:
        hi = mid
report["tube_radius_below_which_all_nearest_are_square"] = round(lo, 3)

print(json.dumps(report, indent=1))