"""R3 — Sensor elevation, terrain occlusion and placement.
Synthetic terrain: 30 fractional-Brownian surfaces (spectral synthesis, beta=3.2), 3 km x 3 km at 15 m cells,
scaled to 12 m std relief, plus 0-60 random building blocks (8-14 m tall, 30-60 m footprint) in the inner area.
Viewshed by radial ray-march (1440 azimuths, 15 m steps), observer at height h above ground, target 1.7 m person.
Experiment A: visible fraction of ground within 1.2 km vs mast height (centre sensor).
Experiment B: k sensors at h=2 m (ground post), 9 m (mast) and perched sentries (roof/high-point, +2 m) chosen
greedily from 80 candidates vs randomly; coverage of the 2.4 km square inner area.
Seeded (seed=11). Outputs research/out/elevation_coverage.json."""
import json
import numpy as np

rng = np.random.default_rng(11)
N, CELL = 200, 15.0
yy, xx = np.mgrid[0:N, 0:N]

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)
    obst = (surface[ix, iy] - z0) / steps                      # slope to obstacle top
    horizon = np.maximum.accumulate(obst, axis=1)
    horizon_prev = np.concatenate([np.full((len(AZ), 1), -np.inf), horizon[:, :-1]], axis=1)
    is_ground = surface[ix, iy] - ground[ix, iy] < 0.5       # building roofs are not "ground"
    tslope = (ground[ix, iy] + target - z0) / steps
    seen = (tslope >= horizon_prev) & inside & is_ground
    vis[ix[seen], iy[seen]] = True
    return vis

dist = np.hypot(xx - 100, yy - 100) * CELL
heights = [2, 5, 9, 15, 25, 40]
A = {h: [] for h in heights}
B = {"post_2m": {"greedy": [], "random": []}, "mast_9m": {"greedy": [], "random": []}, "perch_high": {"greedy": [], "random": []}}
K = 6
inner = (xx >= 20) & (xx < 180) & (yy >= 20) & (yy < 180)
for trial in range(30):
    g, s = terrain(int(rng.integers(0, 60)))
    groundmask = (s - g < 0.5) & (dist <= 1200)
    for h in heights:
        v = viewshed(s, g, 100, 100, h)
        A[h].append(float((v & groundmask).sum() / groundmask.sum()))
    # candidates: 40 highest local cells (incl. roofs) + 40 random, inside inner area
    cand_idx = np.argwhere(inner)
    hi = cand_idx[np.argsort(-s[inner])[:400:10]]
    rnd = cand_idx[rng.choice(len(cand_idx), 40, replace=False)]
    cands = np.vstack([hi, rnd])
    gmask = inner & (s - g < 0.5)
    for mode, h in [("post_2m", 2), ("mast_9m", 9), ("perch_high", 2)]:
        pool = cands if mode == "perch_high" else rnd           # posts/masts go where ground access is (random cells)
        views = [viewshed(s, g, int(a), int(b), h) & gmask for a, b in pool]
        cov = np.zeros((N, N), bool); greedy = []
        chosen = set()
        for k in range(K):
            gains = [(-1 if i in chosen else (v | cov).sum()) for i, v in enumerate(views)]
            i = int(np.argmax(gains)); chosen.add(i); cov |= views[i]
            greedy.append(float(cov.sum() / gmask.sum()))
        B[mode]["greedy"].append(greedy)
        order = rng.permutation(len(views))[:K]; cov = np.zeros((N, N), bool); rr = []
        for i in order:
            cov |= views[i]; rr.append(float(cov.sum() / gmask.sum()))
        B[mode]["random"].append(rr)

def ms(a): a = np.array(a); return [round(float(x), 4) for x in a.mean(0)], [round(float(x), 4) for x in a.std(0)]
out = {"A": [{"h": h, "mean": round(float(np.mean(A[h])), 4), "sd": round(float(np.std(A[h])), 4)} for h in heights], "B": {}}
for mode in B:
    out["B"][mode] = {}
    for strat in B[mode]:
        m, sd = ms(B[mode][strat]); out["B"][mode][strat] = {"mean": m, "sd": sd}
out["horizon_km"] = [{"h": h, "geometric": round(3.57 * h ** 0.5, 2), "refracted": round(4.12 * h ** 0.5, 2)} for h in heights]
json.dump(out, open("out/elevation_coverage.json", "w"), indent=1)
print(json.dumps(out["A"])); [print(m, out["B"][m]["greedy"]["mean"], out["B"][m]["random"]["mean"]) for m in out["B"]]