"""Tacoma Narrows toy model: a deck cross-section on one-sided hangers (Lazer-McKenna / McKenna).
State s = (y, y', theta, theta'); y is the downward heave of the centre, theta > 0 when the right edge is down.
Run:  python3 tacoma_bridge.py      (needs only numpy; takes a few seconds)"""
import numpy as np

m, K, ell, delta, G = 2500.0, 8000.0, 6.0, 0.15, 9.8
y0 = m * G / (2 * K)                       # static stretch of each hanger
wy, wt = np.sqrt(2 * K / m), np.sqrt(6 * K / m)   # heave and torsion frequencies when the hangers stay taut


def rhs(s, t, eps, mu, d=delta, slack=True):
    y, v, th, w = s
    fL, fR = y0 + y - ell * np.sin(th), y0 + y + ell * np.sin(th)
    if slack:
        fL, fR = max(fL, 0.0), max(fR, 0.0)           # a hanger pulls but never pushes
    return np.array([v,
                     -(K / m) * (fL + fR - 2 * y0) - d * v + eps * np.sin(mu * t),
                     w,
                     3 * K * np.cos(th) * (fL - fR) / (m * ell) - d * w])


def step(method, s, t, h, *a, **k):
    f = lambda u, tt: rhs(u, tt, *a, **k)
    if method == "euler":
        return s + h * f(s, t)
    if method == "heun":
        k1 = f(s, t); return s + 0.5 * h * (k1 + f(s + h * k1, t + h))
    k1 = f(s, t); k2 = f(s + 0.5 * h * k1, t + 0.5 * h); k3 = f(s + 0.5 * h * k2, t + 0.5 * h); k4 = f(s + h * k3, t + h)
    return s + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)


def energy(s, slack=True):
    y, v, th, w = s
    fL, fR = y0 + y - ell * np.sin(th), y0 + y + ell * np.sin(th)
    if slack:
        fL, fR = max(fL, 0.0), max(fR, 0.0)
    return 0.5 * m * v**2 + 0.5 * (m * ell**2 / 3) * w**2 + 0.5 * K * (fL**2 + fR**2 - 2 * y0**2) - m * G * y


def run(method, s0, T, h, **k):
    s, t = np.array(s0, float), 0.0
    for i in range(int(round(T / h))):
        s = step(method, s, t, h, **k); t = (i + 1) * h
    return s


def max_theta(eps, mu, T=120.0, h=0.02, stop=0.1, method="rk4", slack=True):
    s, t, big = np.array([0, 0, 0.001, 0.0]), 0.0, 0.0
    for i in range(int(round(T / h))):
        s = step(method, s, t, h, eps, mu, slack=slack); t = (i + 1) * h
        big = max(big, abs(s[2]))
        if not np.isfinite(big) or big > stop:
            break
    return big


print(f"y0 = {y0:.3f} m,  omega_y = {wy:.4f} rad/s (period {2*np.pi/wy:.3f} s),  omega_theta = {wt:.4f} rad/s (period {2*np.pi/wt:.3f} s)")

# 1. energy test: free fall, no wind, no damping
s0 = [-2.5 * y0, 0, 0.3, 0]
print("\nenergy drift over 30 s (max |E - E0| / E0), h = 0.02")
for meth in ("euler", "heun", "rk4"):
    s, t, worst, E0 = np.array(s0, float), 0.0, 0.0, energy(s0)
    for i in range(1500):
        s = step(meth, s, t, 0.02, 0.0, 1.0, d=0.0); t = (i + 1) * 0.02
        worst = max(worst, abs(energy(s) - E0) / E0)
    print(f"  {meth:6s} {worst:.2e}")

# 2. orders of convergence against RK4 with h = 1e-4 (T = 2 s)
print("\nconvergence orders (T = 2 s, error = max(|dy|, ell |dtheta|)), taut vs slack-crossing")
for slack in (False, True):
    ref = run("rk4", s0, 2.0, 1e-4, eps=0.0, mu=1.0, d=0.0, slack=slack)
    for meth in ("euler", "heun", "rk4"):
        hs = np.array([2.0 / n for n in (56, 112, 224, 448, 896)])
        er = [max(abs(a[0] - ref[0]), ell * abs(a[2] - ref[2]))
              for a in (run(meth, s0, 2.0, h, eps=0.0, mu=1.0, d=0.0, slack=slack) for h in hs)]
        p = np.polyfit(np.log(hs[-4:]), np.log(er[-4:]), 1)[0]
        print(f"  {'slack' if slack else 'taut ':5s} {meth:6s} order {p:5.2f}")

# 3. critical wind at one frequency: scan eps upward on a geometric grid, then bisect the first crossing
def critical(ratio, **kw):
    grid = 0.02 * (8 / 0.02) ** (np.arange(14) / 13)
    for i, e in enumerate(grid):
        if max_theta(e, ratio * wy, **kw) > 0.1:
            lo, hi = (grid[i - 1] if i else 0.0), e
            for _ in range(9):
                mid = 0.5 * (lo + hi)
                lo, hi = (lo, mid) if max_theta(mid, ratio * wy, **kw) > 0.1 else (mid, hi)
            return hi
    return None

print("\ncritical gust strength eps (m/s^2) that twists the deck beyond 0.1 rad within 120 s")
for r in (0.8, 1.0, 1.3):
    print(f"  mu/omega_y = {r}:  RK4 {critical(r)}")
print("  mu/omega_y = 1.0, hangers always taut:", critical(1.0, slack=False))