"""Export one synthetic site from the R3 viewshed model for the interactive 3D viewer.
Same terrain generator and ray-march viewshed as elevation_coverage.py (functions copied so the
experiment isn't re-run on import). Seed 23. Grid 200x200 at 15 m, exported at 100x100 (every 2nd cell).
Outputs research/out/viewshed_site.json: surface heights + visibility masks for sensor heights 2, 9, 25 m."""
import json, numpy as np
rng = np.random.default_rng(23)
N, CELL = 200, 15.0
def terrain(buildings):
    kx = np.fft.fftfreq(N)[:, None]; ky = np.fft.fftfreq(N)[None, :]
    k = np.sqrt(kx ** 2 + ky ** 2); k[0, 0] = 1
    spec = (rng.standard_normal((N, N)) + 1j * rng.standard_normal((N, N))) / k ** (3.2 / 2); spec[0, 0] = 0
    z = np.real(np.fft.ifft2(spec)); z = (z - z.mean()) / z.std() * 12.0
    roof = z.copy()
    for _ in range(buildings):
        cx, cy = rng.integers(40, 160, 2); w, l = rng.integers(2, 5, 2); hgt = rng.uniform(8, 14)
        roof[cx:cx + w, cy:cy + l] = z[cx:cx + w, cy:cy + l].max() + hgt
    return z, roof
AZ = np.linspace(0, 2 * np.pi, 1440, endpoint=False)
def viewshed(surface, ground, ox, oy, h, rmax=1200.0, target=1.7):
    vis = np.zeros((N, N), bool); z0 = surface[ox, oy] + h
    steps = np.arange(1, int(rmax / CELL) + 1) * CELL
    px = ox + np.outer(np.cos(AZ), steps) / CELL; py = oy + np.outer(np.sin(AZ), steps) / CELL
    ix = np.clip(np.rint(px).astype(int), 0, N - 1); iy = np.clip(np.rint(py).astype(int), 0, N - 1)
    inside = (px >= 0) & (px < N) & (py >= 0) & (py < N)
    horizon = np.maximum.accumulate((surface[ix, iy] - z0) / steps, axis=1)
    prev = np.concatenate([np.full((len(AZ), 1), -np.inf), horizon[:, :-1]], axis=1)
    is_ground = surface[ix, iy] - ground[ix, iy] < 0.5
    seen = ((ground[ix, iy] + target - z0) / steps >= prev) & inside & is_ground
    vis[ix[seen], iy[seen]] = True; return vis
g, s = terrain(30)
out = {"n": 100, "cell_m": 30, "z": np.round(s[::2, ::2], 1).flatten().tolist(), "building": ((s - g)[::2, ::2] > 0.5).astype(int).flatten().tolist(), "views": {}}
mask = np.hypot(*np.mgrid[0:N, 0:N] - 100) * CELL <= 1200
for h in (2, 9, 25):
    v = viewshed(s, g, 100, 100, h)
    out["views"][str(h)] = v[::2, ::2].astype(int).flatten().tolist()
    gm = mask & (s - g < 0.5)
    out.setdefault("frac", {})[str(h)] = round(float((v & gm).sum() / gm.sum()), 4)
json.dump(out, open("out/viewshed_site.json", "w"))
print(out["frac"])