""" 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 # --------------------------------------------------------------------------- # 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 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) 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}/")