jpmschweitzerandClaude Fable 5
b165c8038d
Guard becomes land_districts <= 1 (both reviewers converged — a lone island definitionally cannot show two distinct directions; same nothing-to-vary condition one value short), with a lone-island vacuous- pass fixture; golden confirmed untouched. Oasis scaling adjudicated as LIVE, not future — GRID_W is already 1024 on main, so ring iterations change 2/4 -> 4/8 today: extracted a pure oasis_ring_iterations() helper pinned by tests at both 512 and 1024, and traced exactly why the determinism hash stayed green (it reads only elevation; the rings touch only biome — a genuinely different array, not a coincidence). The drainage merge-logic question answered byte-precisely: zero logic changed vs main (comment-only diff) — and the deeper dig PROVED the 'isolated basin with another basin to escape to' branch is mathematically unreachable for any connected grid (contracting vertex groups of a connected graph cannot disconnect it), so the comment now states that instead of narrating a divergence that never fires; two direct merge-target tests added regardless. Wrap test renamed to what it actually pins (non-wrap-awareness). D-010 docstring softened to same-process purity, naming the cascade golden as the cross-run layer. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
937 lines
40 KiB
Python
937 lines
40 KiB
Python
"""
|
||
planet_simulation.py
|
||
--------------------
|
||
Terrain simulation stack for the Settled Reach planet generator.
|
||
|
||
Consumes a body_definition dict (output of body_definition_parser.py)
|
||
and produces a terrain dict consumed by planet_renderer.render_globe().
|
||
|
||
Output terrain dict:
|
||
{
|
||
"elevation": float32 (H, W) [0, 1] normalised elevation
|
||
"temperature": float32 (H, W) [0, 1] 0=coldest, 1=hottest
|
||
"moisture": float32 (H, W) [0, 1] 0=driest, 1=wettest
|
||
"biome": int8 (H, W) biome class index
|
||
"surface_water": bool (H, W) ocean/lake mask
|
||
"hillshade": float32 (H, W) [0, 1] lighting from slope+aspect
|
||
"sea_level": float elevation threshold
|
||
}
|
||
|
||
Pipeline:
|
||
1. Elevation - continent mask + domain-warped FBM + tectonic ridges + erosion
|
||
2. Temperature - analytical formula: star + latitude + altitude
|
||
3. Moisture - Hadley cells + ocean proximity + rain shadow
|
||
4. Hillshade - surface normals from elevation gradient
|
||
5. Rivers - downhill carving from moisture-seeded sources
|
||
6. Biome - extended Whittaker lookup + modifier stack
|
||
|
||
Grid: 512 x 256 (longitude x latitude), equirectangular.
|
||
Row 0 = north pole, row 255 = south pole.
|
||
Col 0 = 180W, col 511 = 180E.
|
||
"""
|
||
|
||
import logging
|
||
import math
|
||
import numpy as np
|
||
from scipy.ndimage import gaussian_filter
|
||
|
||
from biome_config import (
|
||
WHITTAKER_TABLE, CLASS_T_BAND, EXOTIC_CLASSES, CRATER_SCALING,
|
||
)
|
||
|
||
log = logging.getLogger(__name__)
|
||
# Canonical heightmap grid (D-202 amended, #963): bumped to 1024×512 so the
|
||
# stored elevation has real mid-scale detail for the lower cascade layers. The
|
||
# elevation noise is resolution-independent (normalized coords + absolute
|
||
# frequencies), so a higher grid samples the SAME continuous terrain at finer
|
||
# density — features keep their physical size and generation stays deterministic.
|
||
# Pixel-unit operations (gaussian sigma, crater radii, filter windows) scale by
|
||
# `GRID_W / 512` so smoothing/morphology behave identically at any resolution.
|
||
GRID_W = 1024
|
||
GRID_H = 512
|
||
|
||
|
||
def oasis_ring_iterations(grid_w: int) -> tuple[int, int]:
|
||
"""Pixel-unit `binary_dilation` iteration counts for the oasis vegetation
|
||
ring (T-964, PR #210 review — hoshe), scaled by `grid_w / 512` per the
|
||
module's GRID_W pixel-unit convention (see the comment above). Base counts
|
||
(2, 4) are the values authored at the original 512-wide grid; `round(...)`
|
||
is used (not floor/ceil) so the scaling is symmetric around exact
|
||
doubling/halving points, and `max(1, ...)` guards against a degenerate
|
||
zero-iteration dilation at very small `grid_w`.
|
||
|
||
Parameterized on `grid_w` (not reading the module-global `GRID_W`
|
||
directly) so it is independently unit-testable at both the historical
|
||
512 baseline and the current production 1024 value without needing to
|
||
monkeypatch the module global.
|
||
"""
|
||
scale = grid_w / 512.0
|
||
ring1_iters = max(1, round(2 * scale))
|
||
ring2_iters = max(1, round(4 * scale))
|
||
return ring1_iters, ring2_iters
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Seeded RNG
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def _rng(seed: int, salt: int = 0) -> np.random.Generator:
|
||
return np.random.default_rng(seed ^ (salt * 2654435761))
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Noise primitives
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def _hash2(x: np.ndarray, y: np.ndarray, seed: int) -> np.ndarray:
|
||
s = np.int64(seed & 0xFFFF)
|
||
h = (x.astype(np.int64) * np.int64(1619) +
|
||
y.astype(np.int64) * np.int64(31337) +
|
||
s * np.int64(6971)) & np.int64(0xFFFFFFFF)
|
||
h = ((h >> 16) ^ h) * np.int64(0x45d9f3b) & np.int64(0xFFFFFFFF)
|
||
h = ((h >> 16) ^ h) * np.int64(0x45d9f3b) & np.int64(0xFFFFFFFF)
|
||
return (h & np.int64(0xFFFF)).astype(np.float32) / 65535.0
|
||
|
||
|
||
def _vnoise(u, v, freq, seed):
|
||
"""Standard 2D value noise — NOT seamless. Use _vnoise_s for longitude axis."""
|
||
uf = u * freq; vf = v * freq
|
||
x0 = np.floor(uf).astype(np.int32); y0 = np.floor(vf).astype(np.int32)
|
||
x1 = x0 + 1; y1 = y0 + 1
|
||
tx = uf - x0; ty = vf - y0
|
||
tx = tx * tx * (3.0 - 2.0 * tx)
|
||
ty = ty * ty * (3.0 - 2.0 * ty)
|
||
v00 = _hash2(x0, y0, seed); v10 = _hash2(x1, y0, seed)
|
||
v01 = _hash2(x0, y1, seed); v11 = _hash2(x1, y1, seed)
|
||
return (v00*(1-tx)*(1-ty) + v10*tx*(1-ty) +
|
||
v01*(1-tx)*ty + v11*tx*ty).astype(np.float32)
|
||
|
||
|
||
def _hash3(x, y, z, seed):
|
||
"""Hash for 3D integer coords."""
|
||
s = np.int64(seed & 0xFFFF)
|
||
h = (x.astype(np.int64) * np.int64(1619) +
|
||
y.astype(np.int64) * np.int64(31337) +
|
||
z.astype(np.int64) * np.int64(49979) +
|
||
s * np.int64(6971)) & np.int64(0xFFFFFFFF)
|
||
h = ((h >> 16) ^ h) * np.int64(0x45d9f3b) & np.int64(0xFFFFFFFF)
|
||
h = ((h >> 16) ^ h) * np.int64(0x45d9f3b) & np.int64(0xFFFFFFFF)
|
||
return (h & np.int64(0xFFFF)).astype(np.float32) / 65535.0
|
||
|
||
|
||
def _vnoise_seamless(u, v, freq, seed):
|
||
"""
|
||
Seamless value noise in the U (longitude) axis only.
|
||
Maps u -> (cos(u*2π), sin(u*2π)) before hashing, so the noise
|
||
field is periodic in U with period 1 — no seam at the date line.
|
||
V (latitude) is not periodic — poles are endpoints, not a loop.
|
||
"""
|
||
# Project U onto a circle: (cx, cy)
|
||
# Divide circle radius by 2π so one full revolution spans the same
|
||
# distance as freq units on the flat V axis — corrects aspect ratio.
|
||
angle = u * (2.0 * math.pi)
|
||
r = freq / (2.0 * math.pi)
|
||
cx = np.cos(angle) * r
|
||
cy = np.sin(angle) * r
|
||
vf = v * freq
|
||
|
||
# Integer lattice in 3D (cx, cy, vf)
|
||
x0 = np.floor(cx).astype(np.int32); x1 = x0 + 1
|
||
y0 = np.floor(cy).astype(np.int32); y1 = y0 + 1
|
||
z0 = np.floor(vf).astype(np.int32); z1 = z0 + 1
|
||
|
||
# Smoothstep weights
|
||
tx = cx - x0; tx = tx * tx * (3.0 - 2.0 * tx)
|
||
ty = cy - y0; ty = ty * ty * (3.0 - 2.0 * ty)
|
||
tz = vf - z0; tz = tz * tz * (3.0 - 2.0 * tz)
|
||
|
||
# Trilinear interpolation over 8 corners
|
||
v000 = _hash3(x0, y0, z0, seed); v100 = _hash3(x1, y0, z0, seed)
|
||
v010 = _hash3(x0, y1, z0, seed); v110 = _hash3(x1, y1, z0, seed)
|
||
v001 = _hash3(x0, y0, z1, seed); v101 = _hash3(x1, y0, z1, seed)
|
||
v011 = _hash3(x0, y1, z1, seed); v111 = _hash3(x1, y1, z1, seed)
|
||
|
||
return (v000*(1-tx)*(1-ty)*(1-tz) + v100*tx*(1-ty)*(1-tz) +
|
||
v010*(1-tx)*ty*(1-tz) + v110*tx*ty*(1-tz) +
|
||
v001*(1-tx)*(1-ty)*tz + v101*tx*(1-ty)*tz +
|
||
v011*(1-tx)*ty*tz + v111*tx*ty*tz).astype(np.float32)
|
||
|
||
|
||
def _fbm(u, v, seed, octaves=6, lacunarity=2.0, gain=0.50, base_freq=2.0):
|
||
"""FBM using seamless noise in U — no longitude seam."""
|
||
result = np.zeros_like(u, dtype=np.float32)
|
||
amp = 1.0; freq = base_freq; total = 0.0
|
||
rng = np.random.default_rng(seed)
|
||
for _ in range(octaves):
|
||
oct_seed = int(rng.integers(0, 0x7FFFFFFF))
|
||
result += amp * _vnoise_seamless(u, v, freq, oct_seed)
|
||
total += amp
|
||
amp *= gain; freq *= lacunarity
|
||
return result / (total + 1e-9)
|
||
|
||
|
||
def _domain_warp(u, v, seed, strength=0.35):
|
||
"""Domain warp using seamless FBM — preserves no-seam property."""
|
||
wu = _fbm(u + 1.7, v + 9.2, seed + 1, octaves=4) * 2.0 - 1.0
|
||
wv = _fbm(u + 8.3, v + 2.8, seed + 2, octaves=4) * 2.0 - 1.0
|
||
# Only warp u periodically — keep v warp non-periodic (poles stay poles)
|
||
return (u + wu * strength) % 1.0, np.clip(v + wv * strength * 0.5, 0.0, 1.0)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Coordinate grids
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def _make_grids():
|
||
u_1d = np.linspace(0, 1, GRID_W, dtype=np.float32)
|
||
v_1d = np.linspace(0, 1, GRID_H, dtype=np.float32)
|
||
u, v = np.meshgrid(u_1d, v_1d)
|
||
lat_frac = -(v - 0.5) * 2.0 # +1 = north, -1 = south
|
||
lon_frac = (u - 0.5) * 2.0
|
||
lat_rad = lat_frac * (math.pi / 2.0)
|
||
return u, v, lat_frac, lon_frac, lat_rad
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 1. Elevation
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def _continent_mask(u, v, seed, land_fraction):
|
||
def _norm(a):
|
||
lo, hi = a.min(), a.max()
|
||
return (a - lo) / (hi - lo + 1e-9)
|
||
|
||
def _contrast(a, strength=3.0):
|
||
"""
|
||
S-curve contrast: pushes highs toward 1 and lows toward 0
|
||
regardless of the field mean. More reliable than power curves
|
||
which behave differently depending on the field's distribution.
|
||
strength controls steepness — higher = sharper separation.
|
||
"""
|
||
# Sigmoid centred at 0.5: f(x) = 1/(1+exp(-k*(x-0.5)))
|
||
k = strength * 8.0
|
||
return 1.0 / (1.0 + np.exp(-k * (a - 0.5)))
|
||
|
||
# Primary: large continental plates
|
||
wu1, wv1 = _domain_warp(u, v, seed, strength=0.45)
|
||
primary = _norm(_fbm(wu1, wv1, seed + 10, octaves=5, gain=0.58, base_freq=1.2))
|
||
|
||
# Secondary: independent medium-scale field.
|
||
# S-curve contrast gives reliable highs and lows regardless of seed.
|
||
wu2, wv2 = _domain_warp(u, v, seed + 11, strength=0.40)
|
||
sec_raw = _norm(_fbm(wu2, wv2, seed + 20, octaves=5, gain=0.55, base_freq=1.8))
|
||
secondary = _contrast(sec_raw, strength=2.5)
|
||
|
||
# Rift: anisotropic thin elongated features
|
||
wu3, wv3 = _domain_warp(u, v, seed + 17, strength=0.30)
|
||
rift = _norm(_fbm(wu3, wv3 * 0.35, seed + 30, octaves=4, gain=0.52, base_freq=3.5))
|
||
|
||
# Multiplicative gate: secondary zeroes kill primary → ocean channels
|
||
separated = primary * (0.4 + secondary * 0.6)
|
||
combined = separated * 0.82 + (rift - 0.5) * 0.18
|
||
|
||
return _norm(combined).astype(np.float32)
|
||
|
||
|
||
def _tectonic_ridges(u, v, seed, n_plates=8):
|
||
rng = _rng(seed, 99)
|
||
px = rng.uniform(0, 1, n_plates).astype(np.float32)
|
||
py = rng.uniform(0, 1, n_plates).astype(np.float32)
|
||
H, W = u.shape
|
||
|
||
# Domain-warp coords before Voronoi — bends ridge positions into curves
|
||
wu1 = _fbm(u * 1.5 + 3.1, v * 1.5 + 7.4, seed + 201, octaves=3,
|
||
gain=0.55, base_freq=1.8) * 2.0 - 1.0
|
||
wv1 = _fbm(u * 1.5 + 8.6, v * 1.5 + 2.2, seed + 202, octaves=3,
|
||
gain=0.55, base_freq=1.8) * 2.0 - 1.0
|
||
wu2 = _fbm(u * 4.0 + 1.3, v * 4.0 + 5.7, seed + 203, octaves=2,
|
||
gain=0.50, base_freq=3.5) * 2.0 - 1.0
|
||
wv2 = _fbm(u * 4.0 + 6.1, v * 4.0 + 0.9, seed + 204, octaves=2,
|
||
gain=0.50, base_freq=3.5) * 2.0 - 1.0
|
||
|
||
uw = (u + wu1 * 0.22 + wu2 * 0.08) % 1.0
|
||
vw = np.clip(v + wv1 * 0.18 + wv2 * 0.06, 0.0, 1.0)
|
||
|
||
dist1 = np.full((H, W), np.inf, dtype=np.float32)
|
||
dist2 = np.full((H, W), np.inf, dtype=np.float32)
|
||
for i in range(n_plates):
|
||
du = np.minimum(np.abs(uw - px[i]), 1.0 - np.abs(uw - px[i]))
|
||
dv = np.abs(vw - py[i])
|
||
d = np.sqrt(du**2 + dv**2)
|
||
mask = d < dist1
|
||
dist2 = np.where(mask, dist1, np.minimum(dist2, d))
|
||
dist1 = np.where(mask, d, dist1)
|
||
|
||
# Two ridge widths: broad ranges + sharp collision zones
|
||
broad = np.exp(-((dist2 - dist1) / 0.06) ** 2) * 0.5
|
||
sharp = np.exp(-((dist2 - dist1) / 0.025) ** 2) * 1.0
|
||
ridge_raw = np.clip(broad + sharp, 0, 1)
|
||
|
||
# Amplitude variation along ridge
|
||
ridge_noise = _fbm(u, v, seed + 50, octaves=4, gain=0.55, base_freq=4.0)
|
||
|
||
# Fracture zones — cross-cutting features (transform faults, rift valleys)
|
||
# Anisotropic: stretch u relative to v for elongated cross features
|
||
fracture = _fbm(u * 0.4, v, seed + 77, octaves=3, gain=0.6, base_freq=6.0)
|
||
fracture = np.clip(fracture - 0.55, 0, 1) * 2.0
|
||
|
||
return np.clip(ridge_raw * (0.35 + 0.65 * ridge_noise)
|
||
+ fracture * 0.20, 0, 1).astype(np.float32)
|
||
|
||
|
||
def _erode(terrain, passes, seed):
|
||
result = terrain.copy()
|
||
# Resolution scale: smoothing radii and the per-pixel slope (which halves as
|
||
# the grid doubles, since np.gradient is in pixel units) scale with width so
|
||
# erosion behaves identically at any GRID size.
|
||
scale = result.shape[1] / 512.0
|
||
for _ in range(passes):
|
||
gy, gx = np.gradient(result)
|
||
slope = np.sqrt(gx**2 + gy**2)
|
||
smooth = gaussian_filter(result, sigma=1.2 * scale)
|
||
weight = np.clip(slope * 6.0 * scale, 0.0, 1.0)
|
||
result = result * (1.0 - weight * 0.35) + smooth * (weight * 0.35)
|
||
gy, gx = np.gradient(result)
|
||
slope = np.sqrt(gx**2 + gy**2)
|
||
flow = gaussian_filter(slope, sigma=3.0 * scale)
|
||
flow = (flow - flow.min()) / (flow.max() - flow.min() + 1e-9)
|
||
result = result - flow * 0.06
|
||
return np.clip(result, 0.0, 1.0)
|
||
|
||
|
||
def compute_elevation(body_def, u, v, lat_frac):
|
||
seed = body_def["seed"]
|
||
planet_class = body_def["planet_class"].replace("_ringed", "")
|
||
land_frac = body_def["terrain"]["land_fraction"]
|
||
tectonics = body_def["terrain"].get("tectonics", "active")
|
||
|
||
plate_map = {"extreme": 12, "active": 8, "low": 5, "none": 3}
|
||
erosion_map = {"extreme": 1, "active": 3, "low": 4, "none": 2}
|
||
n_plates = plate_map.get(tectonics, 8)
|
||
erosion_p = erosion_map.get(tectonics, 3)
|
||
|
||
ocean_pct = (1.0 - land_frac) * 100.0
|
||
detail = _fbm(u, v, seed + 300, octaves=5, gain=0.45, base_freq=4.0)
|
||
|
||
if tectonics == "none":
|
||
# No tectonic activity: gentle base terrain, no ridges, no continents.
|
||
# Craters dominate on these worlds.
|
||
base = _fbm(u, v, seed + 100, octaves=4, gain=0.50, base_freq=1.5)
|
||
elev = base * 0.60 + detail * 0.40
|
||
else:
|
||
# Tectonic worlds: continent mask + ridges scaled by activity level.
|
||
cont = _continent_mask(u, v, seed, land_frac)
|
||
ridges = _tectonic_ridges(u, v, seed, n_plates=n_plates)
|
||
|
||
sea_level_est = float(np.percentile(cont, ocean_pct))
|
||
land_mask = cont >= sea_level_est
|
||
|
||
# Ridge prominence scales with tectonic activity
|
||
ridge_weight = {"low": 0.12, "active": 0.25, "extreme": 0.38}
|
||
rw = ridge_weight.get(tectonics, 0.25)
|
||
|
||
elev = (cont * (0.80 - rw)
|
||
+ ridges * rw * land_mask
|
||
+ detail * 0.20)
|
||
|
||
# Craters happen everywhere. Atmosphere controls how many impactors
|
||
# survive entry; tectonics controls how many craters get resurfaced.
|
||
# Both reduce density, neither toggles craters off entirely.
|
||
crater_factor = (CRATER_SCALING["atmosphere"].get(body_def["physical"]["atmosphere"], 0.25)
|
||
* CRATER_SCALING["tectonics"].get(tectonics, 0.3))
|
||
|
||
if crater_factor > 0.02:
|
||
rng = _rng(seed, 77)
|
||
base_count = CRATER_SCALING["base_count"]
|
||
n_craters = max(5, int(base_count * crater_factor))
|
||
cy_c = rng.uniform(0, GRID_H, n_craters).astype(np.float32)
|
||
cx_c = rng.uniform(0, GRID_W, n_craters).astype(np.float32)
|
||
# Power-law: most craters are small, a few large. Radii are in pixels,
|
||
# so scale with resolution to keep craters the same physical size.
|
||
raw_sizes = rng.power(0.4, n_craters) # skewed toward 0
|
||
sizes = ((2 + raw_sizes * 28) * (GRID_W / 512.0)).astype(np.float32)
|
||
# Depth scales with crater factor — eroded worlds have shallower craters
|
||
depth_scale = 0.5 + 0.5 * crater_factor
|
||
depths = ((0.05 + raw_sizes * 0.20) * depth_scale).astype(np.float32)
|
||
|
||
rows = np.arange(GRID_H, dtype=np.float32)
|
||
cols = np.arange(GRID_W, dtype=np.float32)
|
||
rr, cc = np.meshgrid(rows, cols, indexing='ij')
|
||
# Latitude correction: scale longitude distance by cos(lat) so
|
||
# craters are circular on the sphere, not stretched at the poles.
|
||
lat_rad = (0.5 - rr / GRID_H) * math.pi # +pi/2 at north, -pi/2 at south
|
||
cos_lat = np.cos(lat_rad)
|
||
cos_lat = np.clip(cos_lat, 0.1, 1.0) # avoid division issues at poles
|
||
craters = np.zeros_like(elev)
|
||
for i in range(n_craters):
|
||
dy = rr - cy_c[i]
|
||
dx = cc - cx_c[i]
|
||
# Wrap longitude for craters near the date line
|
||
dx = np.minimum(np.abs(dx), GRID_W - np.abs(dx))
|
||
# Scale dx by cos(lat) at the crater center
|
||
center_lat = (0.5 - cy_c[i] / GRID_H) * math.pi
|
||
dx_scaled = dx / max(math.cos(center_lat), 0.1)
|
||
d = np.sqrt(dy**2 + dx_scaled**2)
|
||
r = sizes[i]
|
||
dep = depths[i]
|
||
# Crater profile: flat floor inside 0.6r, raised rim at 0.9-1.1r,
|
||
# smooth falloff outside. More realistic than gaussian dimple.
|
||
floor = np.clip(1.0 - d / (r * 0.6), 0, 1)
|
||
rim = np.exp(-((d - r) / (r * 0.25))**2)
|
||
craters -= dep * floor * 0.8 # excavate floor
|
||
craters += dep * rim * 0.3 # raise rim
|
||
elev = elev + craters
|
||
elev = np.clip(elev, 0.0, None) # floor at 0
|
||
elif planet_class == "frozen":
|
||
elev = gaussian_filter(elev, sigma=1.5 * (GRID_W / 512.0)).astype(np.float32)
|
||
elif planet_class == "volcanic":
|
||
erosion_p = max(1, erosion_p - 1)
|
||
|
||
elev = _erode(elev, passes=erosion_p, seed=seed)
|
||
|
||
lo, hi = elev.min(), elev.max()
|
||
elev = (elev - lo) / (hi - lo + 1e-9)
|
||
sea_level = float(np.percentile(elev, ocean_pct))
|
||
|
||
# Polar ice — smooth land elevation toward a low plateau at high latitudes.
|
||
# Only applies when there's an atmosphere to deliver precipitation/ice.
|
||
# Airless bodies have no polar caps — cold rock stays rock.
|
||
hydro = body_def.get("environment", {}).get("hydrosphere", "ocean")
|
||
atmo = body_def["physical"]["atmosphere"]
|
||
has_polar_ice = atmo not in ("none",) and hydro not in ("none", "subsurface")
|
||
|
||
if has_polar_ice:
|
||
ice_lat = body_def["terrain"].get("polar_ice_lat", 0.80)
|
||
lat_abs = np.abs(lat_frac)
|
||
ice_blend = np.clip((lat_abs - ice_lat) / (1.0 - ice_lat + 0.01), 0, 1)
|
||
if planet_class == "frozen":
|
||
ice_blend = np.clip(ice_blend * 2.0, 0, 1)
|
||
land_mask = elev >= sea_level
|
||
ice_target = sea_level + 0.05
|
||
elev = np.where(
|
||
land_mask,
|
||
elev * (1.0 - ice_blend * 0.6) + ice_target * (ice_blend * 0.6),
|
||
elev)
|
||
elev = np.clip(elev, 0.0, 1.0).astype(np.float32)
|
||
|
||
sea_level = float(np.percentile(elev, ocean_pct))
|
||
# Dry worlds (no/subsurface hydrosphere): low elevation is dry basin, not ocean.
|
||
hydro = body_def.get("environment", {}).get("hydrosphere", "ocean")
|
||
if hydro in ("none", "subsurface"):
|
||
surf_water = np.zeros_like(elev, dtype=bool)
|
||
else:
|
||
surf_water = elev < sea_level
|
||
return elev, sea_level, surf_water
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 2. Temperature
|
||
# ---------------------------------------------------------------------------
|
||
|
||
# Stellar luminosity relative to Sol (approximate midpoint per spectral type)
|
||
STAR_LUMINOSITY = {
|
||
"O": 100000.0, "B": 1000.0, "A": 10.0,
|
||
"F": 2.5, "G": 1.0, "K": 0.4, "M": 0.04,
|
||
}
|
||
|
||
|
||
# CLASS_T_BAND loaded from biomes.toml via biome_config
|
||
|
||
def compute_temperature(body_def, elevation, sea_level, lat_frac):
|
||
star_type = body_def["star"]["type"]
|
||
distance_au = body_def["orbit"]["distance_au"]
|
||
axial_tilt = body_def["orbit"]["axial_tilt_deg"]
|
||
atmo = body_def["physical"]["atmosphere"]
|
||
planet_class = body_def["planet_class"].replace("_ringed", "")
|
||
geothermal = body_def.get("environment", {}).get("geothermal_flux", "low")
|
||
|
||
# Equilibrium temperature — descriptor-anchored.
|
||
#
|
||
# We compute the raw stellar physics (Stefan-Boltzmann) to get a
|
||
# physically grounded value, then clamp it to the temperature band
|
||
# appropriate for the planet_class. This ensures the wiki's descriptors
|
||
# (temperate, frozen, arid…) are always honoured even when orbital
|
||
# parameters were set with "close enough" precision.
|
||
#
|
||
# Within the clamped band, the raw value still drives relative warmth:
|
||
# a close-in temperate world sits at the warm end of the temperate band,
|
||
# a far-out one at the cool end. The fiction wins; physics sets the gradient.
|
||
lum = body_def.get("star", {}).get("luminosity_solar",
|
||
STAR_LUMINOSITY.get(star_type, 1.0))
|
||
t_raw = 278.5 * (lum ** 0.25) / math.sqrt(max(distance_au, 0.01))
|
||
|
||
greenhouse = {"none": 0, "thin": 8, "standard": 33, "thick": 80}
|
||
t_raw += greenhouse.get(atmo, 0)
|
||
|
||
|
||
temperature_clamped = False
|
||
temperature_raw_K = float(t_raw)
|
||
|
||
if planet_class in CLASS_T_BAND:
|
||
t_lo, t_hi = CLASS_T_BAND[planet_class]
|
||
t_base = float(np.clip(t_raw, t_lo, t_hi))
|
||
if t_raw < t_lo or t_raw > t_hi:
|
||
temperature_clamped = True
|
||
log.debug(f" T_raw={t_raw:.0f}K clamped to [{t_lo},{t_hi}] "
|
||
f"for {planet_class} ({body_def.get('id','')})")
|
||
else:
|
||
t_base = t_raw
|
||
|
||
tilt_factor = 1.0 - (axial_tilt / 90.0) * 0.5
|
||
# Atmosphere controls heat redistribution — thicker atmo = smaller
|
||
# equator-pole gradient. Thin/no atmo = extreme day/night but we
|
||
# still want the planet class to read correctly at the poles.
|
||
atmo_gradient_scale = {"none": 0.6, "thin": 0.7, "standard": 1.0, "thick": 1.2}
|
||
lat_gradient = 60.0 * tilt_factor * atmo_gradient_scale.get(atmo, 1.0)
|
||
t_lat = t_base - lat_gradient * np.abs(lat_frac)
|
||
|
||
max_relief_km = body_def.get("terrain", {}).get("max_elevation_km", 10.0)
|
||
elev_land = np.where(elevation >= sea_level,
|
||
(elevation - sea_level) / (1.0 - sea_level + 1e-9), 0.0)
|
||
elev_km = elev_land * max_relief_km
|
||
lapse = 6.5 if atmo != "none" else 2.0
|
||
t_final = t_lat - lapse * elev_km
|
||
|
||
class_offset = {"frozen": -30, "volcanic": 20, "arid": 10}
|
||
t_final += class_offset.get(planet_class, 0)
|
||
|
||
geo_boost = {"low": 0, "moderate": 5, "high": 15, "extreme": 35}
|
||
t_final += geo_boost.get(geothermal, 0)
|
||
|
||
# Soft floor: prevent planet class from being contradicted at the poles.
|
||
# An arid world shouldn't have ice caps; a volcanic world shouldn't freeze.
|
||
# Clamp the minimum temperature to the class band's lower bound.
|
||
if planet_class in CLASS_T_BAND:
|
||
t_floor = CLASS_T_BAND[planet_class][0]
|
||
t_final = np.maximum(t_final, t_floor)
|
||
|
||
# Return absolute Kelvin grid plus audit metadata.
|
||
# Biome lookup needs absolute values; renderer normalises for display.
|
||
return t_final.astype(np.float32), temperature_clamped, temperature_raw_K
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 3. Moisture
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def compute_moisture(body_def, elevation, sea_level, temperature,
|
||
lat_frac, lon_frac):
|
||
# Normalise temperature locally for moisture computation
|
||
t_norm = np.clip((temperature - temperature.min()) /
|
||
(temperature.max() - temperature.min() + 1e-9), 0, 1)
|
||
atmo = body_def["physical"]["atmosphere"]
|
||
planet_class = body_def["planet_class"].replace("_ringed", "")
|
||
hydro = body_def.get("environment", {}).get("hydrosphere", "none")
|
||
|
||
if atmo == "none":
|
||
return np.zeros((GRID_H, GRID_W), dtype=np.float32)
|
||
|
||
lat_abs = np.abs(lat_frac)
|
||
|
||
# Hadley cell bands
|
||
itcz = np.clip(1.0 - (lat_abs / 0.33), 0, 1)
|
||
subtr = np.clip(1.0 - np.abs(lat_abs - 0.50) / 0.17, 0, 1)
|
||
polar = np.clip((lat_abs - 0.67) / 0.33, 0, 1)
|
||
hadley = np.clip(itcz * 0.85 + subtr * 0.10 + polar * 0.40, 0, 1)
|
||
|
||
# Ocean proximity
|
||
surf_water = elevation < sea_level
|
||
if surf_water.any():
|
||
from scipy.ndimage import distance_transform_edt
|
||
dist = distance_transform_edt(~surf_water).astype(np.float32)
|
||
ocean_prox = 1.0 - np.clip(dist / (dist.max() * 0.5 + 1e-9), 0, 1)
|
||
else:
|
||
ocean_prox = np.zeros((GRID_H, GRID_W), dtype=np.float32)
|
||
|
||
# Rain shadow — westerly winds: windward (west face) is wet
|
||
shift = max(1, GRID_W // 80)
|
||
elev_above = np.clip(elevation - sea_level, 0, None)
|
||
elev_sh = np.clip(np.roll(elevation, shift, axis=1) - sea_level, 0, None)
|
||
shadow_raw = np.clip(elev_sh - elev_above * 0.5, 0, None)
|
||
shadow_raw = shadow_raw / (shadow_raw.max() + 1e-9)
|
||
rain_shadow = 1.0 - shadow_raw * 0.70
|
||
|
||
moisture = (hadley * 0.40
|
||
+ ocean_prox * 0.45
|
||
+ t_norm * 0.15) * rain_shadow
|
||
|
||
class_scale = {
|
||
"arid": 0.25, "oceanic": 1.30, "forest": 1.30,
|
||
"frozen": 0.55, "volcanic": 0.40, "barren": 0.05,
|
||
}
|
||
moisture *= class_scale.get(planet_class, 1.0)
|
||
|
||
hydro_scale = {
|
||
"ocean": 1.2, "liquid_water": 1.2,
|
||
"subsurface": 0.1, "none": 0.05,
|
||
}
|
||
moisture *= hydro_scale.get(hydro, 1.0)
|
||
|
||
moisture = gaussian_filter(moisture.astype(np.float32), sigma=2.0 * (GRID_W / 512.0))
|
||
# Only normalize if the raw range is substantial — otherwise the
|
||
# normalization re-inflates near-zero moisture on dry worlds back to [0,1].
|
||
m_min, m_max = moisture.min(), moisture.max()
|
||
if m_max > 0.05:
|
||
moisture = ((moisture - m_min) / (m_max - m_min + 1e-9)).astype(np.float32)
|
||
else:
|
||
# Effectively dry — clamp to near-zero
|
||
moisture = np.clip(moisture / 0.05, 0, 1).astype(np.float32)
|
||
return moisture
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 4. Hillshade
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def compute_hillshade(elevation,
|
||
sun_azimuth_deg=315.0,
|
||
sun_altitude_deg=45.0):
|
||
scale = GRID_W / 8.0
|
||
gy, gx = np.gradient(elevation * scale)
|
||
mag = np.sqrt(gx**2 + gy**2 + 1.0)
|
||
nx = -gx / mag; ny = -gy / mag; nz = 1.0 / mag
|
||
|
||
az = math.radians(sun_azimuth_deg)
|
||
alt = math.radians(sun_altitude_deg)
|
||
lx = math.cos(alt) * math.cos(az)
|
||
ly = math.cos(alt) * math.sin(az)
|
||
lz = math.sin(alt)
|
||
|
||
diffuse = np.clip(nx * lx + ny * ly + nz * lz, 0.0, 1.0)
|
||
return (0.25 + 0.75 * diffuse).astype(np.float32)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Rivers are NOT computed here (D-208, #963): river networks are derived by the
|
||
# Rust cascade's D8 drainage from the heightmap — the single source of river
|
||
# truth, with mouths that reach the sea by construction. The old heuristic
|
||
# `compute_rivers` was removed to avoid implying the Python sim owns rivers.
|
||
# ---------------------------------------------------------------------------
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 6. Biome
|
||
# ---------------------------------------------------------------------------
|
||
|
||
# WHITTAKER_TABLE, EXOTIC_CLASSES loaded from biomes.toml via biome_config
|
||
|
||
|
||
def compute_biome(body_def, elevation, sea_level, surface_water,
|
||
temperature, moisture):
|
||
H, W = elevation.shape
|
||
biome = np.zeros((H, W), dtype=np.int8)
|
||
land = ~surface_water
|
||
|
||
atmo = body_def["physical"]["atmosphere"]
|
||
|
||
# --- Atmosphere gate ---
|
||
# Worlds with no or thin atmosphere can't support vegetation.
|
||
# Skip the Whittaker table entirely — classify by elevation and
|
||
# temperature only, using rock/dust/ice classes.
|
||
if atmo in ("none", "thin"):
|
||
# Dry terrain classes: 27=dust plain, 28=rocky highland,
|
||
# 29=warm dust, 30=cold rock. No vegetation possible.
|
||
# No ice on airless worlds — cold rock stays rock.
|
||
hydro = body_def.get("environment", {}).get("hydrosphere", "none")
|
||
has_ice_source = hydro not in ("none", "subsurface") or atmo == "thin"
|
||
|
||
elev_norm = np.where(land,
|
||
(elevation - sea_level) / (1.0 - sea_level + 1e-9),
|
||
0.0)
|
||
cf = np.full(land.sum(), 28, dtype=np.int8) # default: rocky highland
|
||
tf = temperature[land].ravel()
|
||
en = elev_norm[land].ravel()
|
||
|
||
# Moon vs planet: moons use grey lunar palette, planets use warm rock
|
||
is_lunar = body_def.get("body_type") == "moon"
|
||
|
||
if is_lunar:
|
||
# Lunar classes: 31=highland, 32=mare (dark basin), 33=midland
|
||
cf[:] = 33 # default: midland grey
|
||
cf[en > 0.50] = 31 # highland
|
||
cf[en < 0.20] = 32 # mare (dark basin floor)
|
||
if has_ice_source:
|
||
cf[tf < 200] = 17 # ice (only if water source)
|
||
else:
|
||
# Temperature-based classification using dry terrain classes
|
||
if has_ice_source:
|
||
cf[tf < 200] = 17 # ice/snow (only if water source)
|
||
else:
|
||
cf[tf < 200] = 30 # cold rock (no water = no ice)
|
||
cf[(tf >= 200) & (tf < 260)] = 30 # cold rock
|
||
cf[(tf >= 260) & (tf < 310)] = 28 # rocky highland
|
||
cf[(tf >= 310) & (tf < 340)] = 29 # warm dust
|
||
cf[tf >= 340] = 15 # hot desert (scorched)
|
||
|
||
# Elevation variation
|
||
if has_ice_source:
|
||
cf[(en > 0.70) & (tf < 273)] = 17 # high + cold = ice cap
|
||
cf[(en < 0.20) & (tf >= 260)] = 27 # low elevation = dust plain
|
||
|
||
biome[land] = cf
|
||
else:
|
||
# --- Standard Whittaker lookup for breathable/toxic atmospheres ---
|
||
tf = temperature[land].ravel()
|
||
mf = moisture[land].ravel()
|
||
cf = np.full(tf.shape, 17, dtype=np.int8) # default: ice
|
||
|
||
# Temperature fed to biome is absolute Kelvin — compare directly
|
||
for (tlo, thi, mlo, mhi, cls) in WHITTAKER_TABLE:
|
||
mask = (tf >= tlo) & (tf <= thi) & (mf >= mlo) & (mf <= mhi)
|
||
cf[mask] = cls
|
||
|
||
biome[land] = cf
|
||
|
||
# Ocean depth bands
|
||
if surface_water.any():
|
||
depth = np.clip((sea_level - elevation) / (sea_level + 1e-9), 0, 1)
|
||
biome[surface_water & (depth < 0.15)] = 2
|
||
biome[surface_water & (depth >= 0.15) & (depth < 0.50)] = 1
|
||
biome[surface_water & (depth >= 0.50)] = 0
|
||
|
||
# Frozen ocean — override ocean biome with ice shelf (class 26).
|
||
# Distinct from land ice (17) — slightly different appearance,
|
||
# blue tint suggests ocean beneath.
|
||
# Add noise to the freeze threshold so the boundary isn't a straight
|
||
# latitude line — ice edges are irregular in reality.
|
||
seed = body_def["seed"]
|
||
u_grid, v_grid, _, _, _ = _make_grids()
|
||
ice_noise = _fbm(u_grid, v_grid, seed + 900, octaves=4,
|
||
gain=0.5, base_freq=3.0) * 2.0 - 1.0
|
||
freeze_threshold = 271.0 + ice_noise * 8.0 # ±8K variation
|
||
frozen_ocean = surface_water & (temperature < freeze_threshold)
|
||
biome[frozen_ocean] = 26
|
||
|
||
# Very cold override — only on worlds with atmosphere (ice needs deposition)
|
||
if atmo not in ("none",):
|
||
biome[(temperature < 243.0) & land] = 17 # below -30C → ice
|
||
|
||
# Elevation overrides — mountain rock and permanent snow.
|
||
# Only apply snow on worlds with atmosphere (ice needs deposition).
|
||
elev_norm = np.where(land,
|
||
(elevation - sea_level) / (1.0 - sea_level + 1e-9),
|
||
0.0)
|
||
hydro_here = body_def.get("environment", {}).get("hydrosphere", "none")
|
||
has_ice_deposition = atmo not in ("none",) and hydro_here not in ("none", "subsurface")
|
||
if has_ice_deposition:
|
||
biome[land & (elev_norm > 0.85)] = 17
|
||
biome[land & (elev_norm > 0.65) & (temperature < 0.35)] = 18
|
||
|
||
# ── Modifier stack ─────────────────────────────────────────────────────
|
||
env = body_def.get("environment", {})
|
||
geothermal = env.get("geothermal_flux", "low")
|
||
chemosyn = env.get("chemosynthetic", False)
|
||
uv_index = env.get("uv_index", "moderate")
|
||
substrate = env.get("substrate", "silicate")
|
||
atmo = body_def["physical"]["atmosphere"]
|
||
planet_class = body_def["planet_class"].replace("_ringed", "")
|
||
|
||
# Geothermal: volcanic worlds get lava/ash at high elevations
|
||
if geothermal in ("extreme", "high") and planet_class == "volcanic":
|
||
biome[land & (elev_norm > 0.75)] = EXOTIC_CLASSES["lava_field"]
|
||
biome[land & (elev_norm > 0.45) & (elev_norm <= 0.75)] = EXOTIC_CLASSES["ash_field"]
|
||
|
||
# Thermophilic fields near heat vents on any high-geothermal world
|
||
if geothermal in ("extreme", "high") and not chemosyn:
|
||
hot = (temperature > 303.0) & land & (elev_norm < 0.45)
|
||
biome[hot] = EXOTIC_CLASSES["thermophilic_field"]
|
||
|
||
# Chemosynthetic worlds (Europa-type): cold surface, geothermal warmth
|
||
if chemosyn:
|
||
geo_warm = (temperature > 263.0) & (temperature < 293.0) & land
|
||
biome[geo_warm] = EXOTIC_CLASSES["chemosynthetic_mat"]
|
||
|
||
# UV radiation: cryptobiotic crust on exposed terrain with thin/no atmo
|
||
if uv_index in ("extreme", "high") and atmo in ("none", "thin"):
|
||
exposed = (land & (elev_norm > 0.15) & (elev_norm < 0.65)
|
||
& (moisture < 0.30)
|
||
& (biome != 17) & (biome != 18) & (biome != 19))
|
||
biome[exposed] = EXOTIC_CLASSES["cryptobiotic_crust"]
|
||
|
||
# Sulfuric substrate: scrub on volcanic mid-elevations
|
||
if substrate == "sulfuric":
|
||
scrub = land & (elev_norm > 0.25) & (elev_norm < 0.65) & (temperature > 0.35)
|
||
biome[scrub & (biome == 18)] = EXOTIC_CLASSES["sulfuric_scrub"]
|
||
|
||
# ── Anomaly scatter ─────────────────────────────────────────────────
|
||
# Sparse micro-features that break biome uniformity and tell stories.
|
||
# A high-frequency noise field selects ~2-5% of cells for anomaly
|
||
# replacement. The anomaly type depends on the surrounding biome context.
|
||
if atmo not in ("none",):
|
||
seed = body_def["seed"]
|
||
u_grid, v_grid, _, _, _ = _make_grids()
|
||
scatter_noise = _fbm(u_grid, v_grid, seed + 800, octaves=3,
|
||
gain=0.6, base_freq=12.0)
|
||
# High threshold = sparse features (~3% of land)
|
||
scatter_mask = (scatter_noise > 0.72) & land
|
||
|
||
if scatter_mask.any():
|
||
b_local = biome[scatter_mask]
|
||
t_local = temperature[scatter_mask]
|
||
m_local = moisture[scatter_mask]
|
||
e_local = elev_norm[scatter_mask]
|
||
new_b = b_local.copy()
|
||
|
||
# Temperate/forest → volcanic vent (lava at high elevation)
|
||
veg_mask = np.isin(b_local, [5, 6, 7, 8, 9, 10, 11])
|
||
new_b[veg_mask & (e_local > 0.50)] = 19 # lava field
|
||
new_b[veg_mask & (e_local > 0.35) & (e_local <= 0.50)] = 25 # ash
|
||
|
||
# Desert/dry → oasis with vegetation ring (only if water exists)
|
||
hydro = body_def.get("environment", {}).get("hydrosphere", "none")
|
||
has_water = hydro not in ("none", "subsurface")
|
||
dry_mask = np.isin(b_local, [13, 14, 15, 27, 28, 29])
|
||
if has_water:
|
||
# Very rare lake in desert lowlands
|
||
new_b[dry_mask & (e_local < 0.10) & (m_local > 0.20)] = 2 # shallow water
|
||
# Vegetation around moisture (oasis fringe — works even without
|
||
# standing water, represents subsurface moisture reaching roots)
|
||
new_b[dry_mask & (m_local > 0.15) & (e_local >= 0.10)] = 7 # savanna
|
||
|
||
# Frozen → geothermal hotspot with pioneer vegetation
|
||
cold_mask = np.isin(b_local, [16, 17])
|
||
new_b[cold_mask & (t_local > 260)] = 12 # shrubland (hardy plants)
|
||
|
||
# Volcanic → cooling zone with pioneer life
|
||
lava_mask = np.isin(b_local, [19, 25])
|
||
new_b[lava_mask & (t_local < 310) & (m_local > 0.30)] = 24 # lithic pioneer
|
||
|
||
biome[scatter_mask] = new_b
|
||
|
||
# Vegetation ring around oasis lakes: dilate water cells from the
|
||
# scatter pass and assign graduated vegetation to the ring.
|
||
# water → coast vegetation → savanna/shrub → original biome
|
||
oasis_water = (biome == 2) & land # scattered lake cells on land
|
||
if oasis_water.any():
|
||
from scipy.ndimage import binary_dilation
|
||
# Pixel-unit iteration counts scale by GRID_W/512 (module
|
||
# convention, see the GRID_W comment above) so the vegetation
|
||
# ring covers the same physical distance at any resolution — a
|
||
# fixed pixel count would shrink the ring (relative to the
|
||
# planet's real size) as GRID_W grows past the 512 baseline.
|
||
# LIVE BEHAVIOR CHANGE (T-964, PR #210 review — hoshe): GRID_W is
|
||
# ALREADY 1024 in production (not a future value), so this scales
|
||
# ring1/ring2 from the historical 2/4 iterations to 4/8 on EVERY
|
||
# planet generated today — see `oasis_ring_iterations`'s test
|
||
# coverage (`test_planet_simulation.py`) pinning both the 512
|
||
# (historical) and 1024 (current production) values explicitly.
|
||
ring1_iters, ring2_iters = oasis_ring_iterations(GRID_W)
|
||
ring1 = binary_dilation(oasis_water, iterations=ring1_iters) & ~oasis_water & land
|
||
ring2 = (
|
||
binary_dilation(oasis_water, iterations=ring2_iters)
|
||
& ~oasis_water
|
||
& ~ring1
|
||
& land
|
||
)
|
||
# Inner ring: lush vegetation (coast/lowland green)
|
||
biome[ring1] = 4 # lowland
|
||
# Outer ring: transitional (savanna/shrub)
|
||
biome[ring2] = 12 # shrubland
|
||
|
||
return biome
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# Top-level simulate()
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def simulate(body_def: dict) -> dict:
|
||
"""
|
||
Run the full simulation stack for one body.
|
||
|
||
Parameters
|
||
----------
|
||
body_def : dict — from body_definition_parser.parse_system()
|
||
|
||
Returns
|
||
-------
|
||
dict terrain dict consumed by planet_renderer.render_globe()
|
||
Empty dict for gas giants (renderer handles those procedurally).
|
||
"""
|
||
planet_class = body_def.get("planet_class", "barren").replace("_ringed", "")
|
||
if planet_class == "gas_giant":
|
||
return {}
|
||
|
||
u, v, lat_frac, lon_frac, lat_rad = _make_grids()
|
||
|
||
elevation, sea_level, surface_water = compute_elevation(
|
||
body_def, u, v, lat_frac)
|
||
|
||
temperature, temp_clamped, temp_raw_K = compute_temperature(
|
||
body_def, elevation, sea_level, lat_frac)
|
||
|
||
moisture = compute_moisture(
|
||
body_def, elevation, sea_level, temperature, lat_frac, lon_frac)
|
||
|
||
hillshade = compute_hillshade(elevation)
|
||
|
||
biome = compute_biome(
|
||
body_def, elevation, sea_level, surface_water, temperature, moisture)
|
||
|
||
# Normalise temperature to [0,1] for renderer display — biome already computed
|
||
t_min, t_max = temperature.min(), temperature.max()
|
||
temperature_norm = ((temperature - t_min) / (t_max - t_min + 1e-9)).astype(np.float32)
|
||
|
||
return {
|
||
"elevation": elevation,
|
||
"temperature": temperature_norm, # normalised [0,1] for renderer
|
||
"moisture": moisture,
|
||
"biome": biome,
|
||
"surface_water": surface_water,
|
||
"hillshade": hillshade,
|
||
"sea_level": sea_level,
|
||
"_grid_w": GRID_W,
|
||
"_grid_h": GRID_H,
|
||
# Audit trail
|
||
"temperature_clamped": temp_clamped,
|
||
"temperature_raw_K": round(temp_raw_K, 1),
|
||
"temperature_band_K": list(CLASS_T_BAND.get(
|
||
body_def.get("planet_class","").replace("_ringed",""), [None,None])),
|
||
}
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# CLI
|
||
# ---------------------------------------------------------------------------
|
||
|
||
if __name__ == "__main__":
|
||
import sys
|
||
import json
|
||
import time
|
||
import os
|
||
from PIL import Image
|
||
|
||
if len(sys.argv) < 2:
|
||
print("Usage: python3 planet_simulation.py body_def.json [--save-grids]")
|
||
sys.exit(1)
|
||
|
||
with open(sys.argv[1]) as f:
|
||
bd = json.load(f)
|
||
|
||
save_grids = "--save-grids" in sys.argv
|
||
|
||
print(f"Simulating: {bd['id']} ({bd['planet_class']})")
|
||
t0 = time.time()
|
||
terrain = simulate(bd)
|
||
|
||
if not terrain:
|
||
print("Gas giant — no terrain simulation.")
|
||
sys.exit(0)
|
||
|
||
dt = time.time() - t0
|
||
print(f"Done in {dt:.1f}s")
|
||
print(f" sea_level: {terrain['sea_level']:.3f}")
|
||
print(f" land cells: {(~terrain['surface_water']).sum()}")
|
||
|
||
ids, counts = np.unique(terrain['biome'], return_counts=True)
|
||
print(f" biomes: {list(zip(ids.tolist(), counts.tolist()))}")
|
||
|
||
if save_grids:
|
||
out = f"/tmp/{bd['id']}_grids"
|
||
os.makedirs(out, exist_ok=True)
|
||
for name in ("elevation", "temperature", "moisture", "hillshade"):
|
||
arr = terrain[name]
|
||
Image.fromarray((arr * 255).astype("uint8"), "L").save(
|
||
f"{out}/{name}.png")
|
||
print(f"Grids saved → {out}/")
|