"""Trapped note in a pipe: lowest antisymmetric eigenvalue for a thin plate on the centreline.

Channel -1<y<1 with hard walls; plate on y=0 for |x|<a. Antisymmetric in y, symmetric in x.
Cell-centred finite differences on the quarter domain 0<x<L, 0<y<1:
  hard wall at y=1 and on the plate (no link), u odd across y=0 off the plate,
  even across x=0, u=0 far away at x=L. Below the cutoff k=pi/2 nothing can leak, so an
  eigenvalue k<pi/2 is a trapped mode.
"""
import numpy as np, scipy.sparse as sp, scipy.sparse.linalg as sl, json, sys

def trapped_k(a, h, L=12.0):
    nx, ny = int(round(L / h)), int(round(1 / h))
    N = nx * ny; idx = lambda i, j: j * nx + i
    rows, cols, vals = [], [], []
    for j in range(ny):
        for i in range(nx):
            p = idx(i, j); d = 0.0; x = (i + 0.5) * h
            for di, dj in ((1, 0), (-1, 0), (0, 1), (0, -1)):
                ii, jj = i + di, j + dj
                if ii < 0: continue                     # even across x=0
                if ii >= nx: d += 2; continue           # u=0 at x=L (ghost -u)
                if jj >= ny: continue                   # hard wall y=1
                if jj < 0:
                    if x < a: continue                  # plate: hard
                    d += 2; continue                    # odd across y=0
                rows.append(p); cols.append(idx(ii, jj)); vals.append(-1.0); d += 1
            rows.append(p); cols.append(p); vals.append(d)
    A = sp.csr_matrix((vals, (rows, cols)), shape=(N, N)) / h**2
    w = sl.eigsh(A, k=1, sigma=0, which='LM', return_eigenvectors=False)
    return float(np.sqrt(w[0]))

if __name__ == '__main__':
    # writes data/trapped.json: python3 code/trapped.py data/trapped.json
    out = {}
    for a in (0.5, 1.0, 2.0):
        k1, k2, k3 = (trapped_k(a, h) for h in (0.05, 0.025, 0.0125))
        out[str(a)] = {'h0.05': round(k1, 5), 'h0.025': round(k2, 5), 'h0.0125': round(k3, 5), 'extrapolated': round(2 * k3 - k2, 4)}
    json.dump({'note': 'k for half-width 1; cutoff pi/2. From code/trapped.py (cell-centred FD, eigsh). extrapolated = 2*k(h/2)-k(h), assuming error ~h.',
               'cutoff': 1.5708, 'plates': out}, open(sys.argv[1] if len(sys.argv) > 1 else 'trapped.json', 'w'), indent=1)