"""Millennium Bridge toy model (Strogatz, Abrams, McRobie, Eckhardt and Ott, Nature 2005): a bridge mode driven by walkers who adjust their pace to the sway.
Bridge:   M X'' + B X' + K X = G sum sin(Theta_i)       Walker:  Theta_i' = omega_i + C A sin(Psi - Theta_i + alpha),  X = A sin Psi, X' = Omega0 A cos Psi.
Run:  python3 mb_bridge.py      (needs only numpy; takes about ten seconds)"""
import numpy as np

M, f0, zeta, G, C, alpha, gamma = 1.13e5, 0.8, 0.01, 35.0, 16.0, np.pi / 2, 0.8   # illustrative values
Om = 2 * np.pi * f0
B, K = 2 * zeta * M * Om, M * Om**2
T = min(5 * gamma, 0.8 * Om)                      # the Lorentzian of walking frequencies is cut at +-T
Z = 2 / np.pi * np.arctan(T / gamma)              # share of the Lorentzian that is kept
Nc = Z * 2 * gamma * B * Om / (G * C * np.sin(alpha))
print(f"omega0 = {Om:.3f} rad/s, B = {B:.0f} kg/s, Z = {Z:.4f}, critical crowd Nc = {Nc:.1f}")


def walkers(N):
    """natural frequencies at the quantiles of the cut Lorentzian, golden-ratio phases (no random numbers)"""
    u = (np.arange(N) + 0.5) / N
    w = Om + gamma * np.tan(np.arctan(T / gamma) * (2 * u - 1))
    th = 2 * np.pi * ((np.arange(N) * 0.6180339887498949) % 1)
    p = (np.arange(N) * 7919) % N, (np.arange(N) * 4973) % N
    return w[p[1]], th[p[0]]


def rhs(s, w, zeta=zeta):
    X, V, th = s[0], s[1], s[2:]
    A, psi = np.hypot(X, V / Om), np.arctan2(X, V / Om)
    return np.concatenate(([V, (G * np.sin(th).sum() - 2 * zeta * M * Om * V - K * X) / M], w + C * A * np.sin(psi - th + alpha)))


def rk4(s, w, h, **k):
    k1 = rhs(s, w, **k); k2 = rhs(s + h / 2 * k1, w, **k); k3 = rhs(s + h / 2 * k2, w, **k); k4 = rhs(s + h * k3, w, **k)
    return s + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)


def wobble(N, Tend=1500.0, h=0.08, **k):
    """mean amplitude A (m) and order parameter r over the last 150 s"""
    w, th = walkers(N); s = np.concatenate(([1e-3, 0.0], th)); n = round(Tend / h); a = []; r = []
    for i in range(n):
        s = rk4(s, w, h, **k)
        if i * h > Tend - 150 and i % 5 == 0:
            A, psi = np.hypot(s[0], s[1] / Om), np.arctan2(s[0], s[1] / Om); a.append(A); r.append(abs(np.exp(1j * (s[2:] - psi)).mean()))
    return np.mean(a), np.mean(r)


print("\n crowd   A (cm) simulated   A (cm) theory   r simulated   r theory")
for N in (100, 150, 200, 250):
    A, r = wobble(N)
    At = G / (Z * B * Om) * np.sqrt(N * (N - Nc)) if N > Nc else 0.0
    rt = np.sqrt(1 - Nc / N) / Z if N > Nc else 0.0
    print(f" {N:5d}   {100 * A:12.2f}   {100 * At:13.2f}   {r:11.3f}   {rt:8.3f}")

print("\nthe fix: with zeta = 3 % the critical crowd is", round(3 * Nc, 0), "- a crowd of 250 no longer wobbles")
A, r = wobble(250, zeta=0.03)
print(f"  simulated: A = {100 * A:.2f} cm, r = {r:.2f}")

# a bad solver can fake a wobble: the bridge alone (no walkers), explicit Euler
for h in (0.002, 0.004, 0.01):
    s = np.array([0.01, 0.0])
    for _ in range(round(300 / h)):
        d = np.array([s[1], (-2 * zeta * M * Om * s[1] - K * s[0]) / M]); s = s + h * d
    print(f"Euler, bridge alone, h = {h}: amplitude after 300 s = {np.hypot(s[0], s[1] / Om):.2e} m (exact decay: {0.01 * np.exp(-zeta * Om * 300):.1e} m)")