1 files changed, 179 insertions, 0 deletions
diff --git a/mapgen/plateaus.py b/mapgen/plateaus.py
new file mode 100644
index 0000000..9647d04
--- /dev/null
+++ b/mapgen/plateaus.py
@@ -0,0 +1,179 @@
+"""Sunken plateaus: Kerguelen-type continental crust that never rose. Outline,
+surface and volcanic field are functions of position and the plateau's config (deterministic per name), so the world
+build, the seabed pass and the viewer agree at any resolution."""
+from __future__ import annotations
+
+import numpy as np
+
+from .graph import distance_to
+from .noise import fbm, name_seed
+from .sphere import east_north, gc_dist_km, great_circle_point, latlon_to_xyz, tangent_dir, xyz_to_latlon
+from .zones import smootherstep
+
+EDGE_WARP = 0.15        # outline radius ± 15 % (fractal noise): bays and lobes, like a small continent
+MARGIN_KM = 150.0       # the surface blends down to the surrounding sea floor over this distance inside the outline
+RELIEF_M = 400.0        # internal relief (ridges, basins)
+HIDDEN_MAX_M = -550.0   # hidden plateaus: every cell at least this deep (checked at −500 m after erosion's re-solve)
+SITE_RULES = (("land", 400.0), ("a plate boundary", 150.0), ("the trench", 200.0))   # km from the outline
+
+
+def semi_axes(p: dict):
+    """(long, short) semi-axes (km) of the plateau's ellipse: π·a·b = area_km2, a / b = elongation."""
+    e = float(p.get("elongation", 1.0))
+    b = float(np.sqrt(p["area_km2"] / (np.pi * e)))
+    return e * b, b
+
+
+def _near(xyz, p, radius_km, pad_km=0.0):
+    a, _ = semi_axes(p)
+    lim = min(np.pi, (a * (1.0 + EDGE_WARP) / (1.0 - EDGE_WARP) + pad_km) / radius_km)
+    return xyz @ latlon_to_xyz(*p["center"]) > np.cos(lim)
+
+
+def rho(xyz, p: dict, seed: int, radius_km: float):
+    """Normalised radius of points: 0 at the centre, < 1 inside the noise-warped elliptical outline."""
+    xyz = np.asarray(xyz, dtype=np.float64)
+    c = latlon_to_xyz(*p["center"])
+    e, n = east_north(c[None])
+    ang = np.arccos(np.clip(xyz @ c, -1.0, 1.0)) * radius_km            # km from the centre along the surface
+    t = tangent_dir(np.broadcast_to(c, xyz.shape), xyz)
+    x, y = ang * (t @ e[0]), ang * (t @ n[0])                            # east, north (azimuthal equidistant)
+    az = np.radians(p.get("azimuth_deg", 0.0))
+    along, across = x * np.sin(az) + y * np.cos(az), x * np.cos(az) - y * np.sin(az)
+    a, b = semi_axes(p)
+    r = np.sqrt((along / a) ** 2 + (across / b) ** 2)
+    w = np.clip(2.0 * fbm(xyz, name_seed(seed, p["name"]), 4, 25.0), -1.0, 1.0)
+    return r / (1.0 + EDGE_WARP * w)
+
+
+def cell_ids(xyz, plateaus: list, seed: int, radius_km: float):
+    """Plateau index per point (−1 outside every plateau)."""
+    xyz = np.asarray(xyz, dtype=np.float64)
+    ids = np.full(len(xyz), -1, np.int16)
+    for k, p in enumerate(plateaus):
+        idx = np.flatnonzero(_near(xyz, p, radius_km))
+        if len(idx):
+            ids[idx[rho(xyz[idx], p, seed, radius_km) < 1.0]] = k
+    return ids
+
+
+def _place(rng, p, seed, radius_km, origin, max_km, inside=0.85, tries=30):
+    """A random point within max_km of origin inside the outline (rho < inside); None if none is found."""
+    e, n = east_north(origin[None])
+    for _ in range(tries):
+        az, dist = rng.uniform(0.0, 2.0 * np.pi), rng.uniform(0.0, max_km)
+        q = great_circle_point(origin, np.sin(az) * e[0] + np.cos(az) * n[0], dist, radius_km)
+        if rho(q[None], p, seed, radius_km)[0] < inside:
+            return q
+    return None
+
+
+def _ll(q):
+    lat, lon = xyz_to_latlon(np.asarray(q, dtype=np.float64))
+    return round(float(lat), 4), round(float(lon), 4)
+
+
+def features(p: dict, seed: int, radius_km: float) -> dict:
+    """The plateau's volcanic field at world scale (deterministic per name): its own hotspot point, 6–20 cones and
+    1–3 calderas (none when vent = 0); island plateaus lift 3–6 of the cones to +0.6…+1.5 km."""
+    rng = np.random.default_rng(name_seed(seed, p["name"]) + 7)
+    c = latlon_to_xyz(*p["center"])
+    _, b = semi_axes(p)
+    vent = float(p.get("vent", 1.0))
+    hot = _place(rng, p, seed, radius_km, c, 0.3 * b)
+    hot = c if hot is None else hot
+    cones, calderas = [], []
+    if vent > 0:
+        for _ in range(int(rng.integers(6, 21))):
+            q = _place(rng, p, seed, radius_km, hot, 0.6 * b)
+            r_km, h = float(rng.uniform(15.0, 40.0)), float(rng.uniform(500.0, 2000.0)) * vent
+            if q is not None:
+                lat, lon = _ll(q)
+                cones.append({"lat": lat, "lon": lon, "radius_km": r_km, "height_m": h, "island": False})
+        for _ in range(int(rng.integers(1, 4))):
+            q = _place(rng, p, seed, radius_km, hot, 0.5 * b)
+            r_km = float(rng.uniform(20.0, 50.0))
+            if q is not None:
+                lat, lon = _ll(q)
+                calderas.append({"lat": lat, "lon": lon, "radius_km": r_km, "rim_m": 300.0 * vent,
+                                 "floor_m": -400.0 * vent})
+    if p.get("islands", False) and cones:
+        k = min(len(cones), int(rng.integers(3, 7)))
+        for i in rng.choice(len(cones), size=k, replace=False):
+            cones[int(i)].update(island=True, peak_m=float(rng.uniform(600.0, 1500.0)),
+                                 radius_km=float(rng.uniform(40.0, 70.0)))
+    return {"name": p["name"], "center": [float(v) for v in p["center"]], "hotspot": list(_ll(hot)),
+            "cones": cones, "calderas": calderas}
+
+
+def surface(xyz, p: dict, feat: dict, seed: int, radius_km: float, z_floor):
+    """Heights of points inside the outline: the plateau top (depth across top_m by low-frequency noise, ± RELIEF_M)
+    plus its cones and calderas, blended down to z_floor over MARGIN_KM inside the edge."""
+    xyz = np.asarray(xyz, dtype=np.float64)
+    s = name_seed(seed, p["name"])
+    t = np.clip(0.5 + 1.5 * fbm(xyz, s + 1, 3, 8.0), 0.0, 1.0)
+    top0, top1 = p["top_m"]
+    z = -(top0 + (top1 - top0) * t) + RELIEF_M * np.clip(2.5 * fbm(xyz, s + 2, 5, 40.0), -1.0, 1.0)
+    for c in feat["cones"]:
+        if not c["island"]:
+            d = gc_dist_km(xyz, latlon_to_xyz(c["lat"], c["lon"]), radius_km)
+            z = z + c["height_m"] * np.clip(1.0 - d / c["radius_km"], 0.0, 1.0) ** 1.5
+    for c in feat["calderas"]:
+        d = gc_dist_km(xyz, latlon_to_xyz(c["lat"], c["lon"]), radius_km)
+        r = c["radius_km"]
+        z = z + c["rim_m"] * np.exp(-((d - r) / (0.25 * r)) ** 2) + c["floor_m"] * (1.0 - smootherstep(d / r))
+    _, b = semi_axes(p)
+    w = smootherstep((1.0 - rho(xyz, p, seed, radius_km)) * b / MARGIN_KM)
+    return np.asarray(z_floor, dtype=np.float64) + (z - z_floor) * w
+
+
+def islands(xyz, feat: dict, radius_km: float, z):
+    """Island cones lift the ground to their peak_m (small volcanic islands); only raises."""
+    z = np.asarray(z, dtype=np.float64)
+    for c in feat["cones"]:
+        if c["island"]:
+            d = gc_dist_km(np.asarray(xyz, dtype=np.float64), latlon_to_xyz(c["lat"], c["lon"]), radius_km)
+            f = np.clip(d / c["radius_km"], 0.0, 1.0)
+            z = np.where(d < c["radius_km"], np.maximum(z, c["peak_m"] - (c["peak_m"] - z) * f ** 1.2), z)
+    return z
+
+
+def apply(g, z, plateaus: list, plateau_id, seed: int):
+    """Grid heights with the plateaus set (after the world's sea-level solve): surfaces and volcanic fields; hidden
+    plateaus clamped to HIDDEN_MAX_M; islands, each island's nearest cell raised to its peak (so every resolution
+    keeps it)."""
+    z = np.asarray(z, dtype=np.float64).copy()
+    R = g.radius_km
+    for k, p in enumerate(plateaus):
+        idx = np.flatnonzero(np.asarray(plateau_id) == k)
+        if len(idx) == 0:
+            continue
+        feat = features(p, seed, R)
+        zk = surface(g.xyz[idx], p, feat, seed, R, z[idx])
+        if p.get("islands", False):
+            z[idx] = islands(g.xyz[idx], feat, R, zk)
+            for c in feat["cones"]:
+                if c["island"]:
+                    i = g.cell_index(c["lat"], c["lon"])
+                    z[i] = max(z[i], c["peak_m"])
+        else:
+            z[idx] = np.minimum(zk, HIDDEN_MAX_M)
+    return z
+
+
+def site_report(g, data: dict, plateaus: list) -> list:
+    """Plateaus that no longer fit their site on this world, as warnings: outline ≥ 400 km from land,
+    ≥ 150 km from plate boundaries, ≥ 200 km from the sketch's trench. Reported, never moved."""
+    ids = np.asarray(data["plateau_id"])
+    masks = {"land": ~np.asarray(data["ocean"]) & (ids < 0), "a plate boundary": np.asarray(data["bnd_type"]) > 0,
+             "the trench": np.asarray(data.get("sk_trench", np.zeros(g.n))) > 0.5}
+    out = []
+    for what, lim in SITE_RULES:
+        if not masks[what].any():
+            continue
+        d = distance_to(g, masks[what])
+        for k, p in enumerate(plateaus):
+            m = ids == k
+            if m.any() and float(d[m].min()) < lim:
+                out.append(f"plateau {p['name']}: {float(d[m].min()):.0f} km from {what} (rule ≥ {lim:.0f} km)")
+    return out