jpmschweitzerandClaude Opus 4.6
f84275edf3
Full terrain-to-render pipeline at tooling/planet-gen/: - body_definition_parser: reads system index.md → body definitions - planet_simulation: FBM + Voronoi ridges, Whittaker biome classification, D8 river routing, physics-driven craters (atmo × tectonics scaling) - render_heightmap: 4096×2048 cartographic maps with hillshade - planet_renderer: 512×512 globe with terrain UV mapping, clouds, rings - biomes.toml: externalized color/classification tables (single source) - scaffold_bodies: creates per-body index.md with YAML frontmatter - batch: unattended processing with error handling, resume, determinism check Supports all body types: temperate, arid, frozen, volcanic, barren, oceanic, gas giant (banded + ringed), and moons (half-size, grey, cratered). Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
957 lines
40 KiB
Python
957 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
|
|
"river_grid": bool (H, W) river cell mask
|
|
"rivers": list of [(row,col), ...] polylines in grid coords
|
|
"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__)
|
|
GRID_W = 512
|
|
GRID_H = 256
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# 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()
|
|
for _ in range(passes):
|
|
gy, gx = np.gradient(result)
|
|
slope = np.sqrt(gx**2 + gy**2)
|
|
smooth = gaussian_filter(result, sigma=1.2)
|
|
weight = np.clip(slope * 6.0, 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)
|
|
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 (2-5 cells), a few are large (15-30)
|
|
raw_sizes = rng.power(0.4, n_craters) # skewed toward 0
|
|
sizes = (2 + raw_sizes * 28).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).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)
|
|
# 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)
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# 5. Rivers
|
|
# ---------------------------------------------------------------------------
|
|
|
|
def compute_rivers(body_def, elevation, sea_level, moisture,
|
|
max_rivers=12):
|
|
atmo = body_def["physical"]["atmosphere"]
|
|
planet_class = body_def["planet_class"].replace("_ringed", "")
|
|
hydro = body_def.get("environment", {}).get("hydrosphere", "none")
|
|
|
|
if atmo == "none" or hydro in ("none", "subsurface", "ice"):
|
|
return []
|
|
|
|
river_cap = {"barren": 2, "volcanic": 3, "arid": 3, "frozen": 2}
|
|
max_rivers = river_cap.get(planet_class, max_rivers)
|
|
|
|
H, W = elevation.shape
|
|
land_mask = elevation >= sea_level
|
|
seed = body_def["seed"]
|
|
rng = _rng(seed, 500)
|
|
|
|
from scipy.ndimage import maximum_filter
|
|
local_max = (elevation == maximum_filter(elevation, size=8)) & land_mask
|
|
moist_ok = moisture > 0.35
|
|
candidates = np.argwhere(local_max & moist_ok)
|
|
if len(candidates) == 0:
|
|
candidates = np.argwhere(land_mask)
|
|
|
|
np.random.default_rng(seed).shuffle(candidates)
|
|
sources = candidates[:min(max_rivers, len(candidates))]
|
|
|
|
D8 = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]
|
|
rivers = []
|
|
|
|
for src in sources:
|
|
r, c = int(src[0]), int(src[1])
|
|
path = [(r, c)]
|
|
visited = {(r, c)}
|
|
|
|
for _ in range(GRID_W * 2):
|
|
if elevation[r, c] < sea_level:
|
|
break
|
|
best_drop = 0.0; best_nr = -1; best_nc = -1
|
|
for dr, dc in D8:
|
|
nr = r + dr; nc = (c + dc) % W
|
|
if nr < 0 or nr >= H or (nr, nc) in visited:
|
|
continue
|
|
drop = elevation[r, c] - elevation[nr, nc]
|
|
drop += float(rng.uniform(-0.005, 0.005))
|
|
if drop > best_drop:
|
|
best_drop = drop; best_nr = nr; best_nc = nc
|
|
if best_nr < 0:
|
|
break
|
|
r, c = best_nr, best_nc
|
|
visited.add((r, c))
|
|
path.append((r, c))
|
|
|
|
if len(path) > 5:
|
|
rivers.append(path)
|
|
|
|
return rivers
|
|
|
|
|
|
def _rivers_to_grid(rivers, H, W):
|
|
grid = np.zeros((H, W), dtype=bool)
|
|
for path in rivers:
|
|
for r, c in path:
|
|
if 0 <= r < H and 0 <= c < W:
|
|
grid[r, c] = True
|
|
return grid
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# 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
|
|
ring1 = binary_dilation(oasis_water, iterations=2) & ~oasis_water & land
|
|
ring2 = binary_dilation(oasis_water, iterations=4) & ~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)
|
|
|
|
rivers = compute_rivers(body_def, elevation, sea_level, moisture)
|
|
river_grid = _rivers_to_grid(rivers, GRID_H, GRID_W)
|
|
|
|
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,
|
|
"river_grid": river_grid,
|
|
"rivers": rivers,
|
|
"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, json, time, 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()}")
|
|
print(f" rivers: {len(terrain['rivers'])} polylines")
|
|
|
|
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}/")
|