# Computes V(ell)/ell of the crossover theorem (t:crossover) by quadrature, for the figure f:intro-crossover.
# Run with: python3 crossover_ratio.py   (writes crossover_ratio.tsv; the figure inlines the same table)
import numpy as np
from scipy.special import i0e
from scipy.integrate import quad
def V(l):
    f = lambda u: i0e(2.0 / u) * np.sin(np.pi * l * u)**2 / (np.pi * u)**2
    tot, _ = quad(f, 1e-12, 1.0 / l, limit=200)
    b = 1.0 / l
    k = 1
    while True:
        c = (k + 1) / l
        v, _ = quad(f, b, c, limit=100)
        tot += v
        b = c
        k += 1
        if c > 2e3 and k > 50:
            break
    # tail with S(u) = i0e(2/u) ~ 1 - 2/u ; average of sin^2 = 1/2
    tot += 1.0 / (2 * np.pi**2 * b) - 2.0 / (4 * np.pi**2 * b**2)
    return 2 * tot
ls = 10.0 ** np.linspace(-3, 3, 37)
out = []
prev = 2.0
for l in ls:
    v = V(l)
    r = v / l
    assert r < prev and 0 < r < 1, (l, r)
    prev = r
    out.append((l, r))
with open('crossover_ratio.tsv', 'w') as fh:
    fh.write('ell\tratio\n')
    for l, r in out:
        fh.write(f'{l:.6e}\t{r:.6e}\n')
print(out[0], out[18], out[-1])