#!/usr/bin/env python3
"""
WHAT ACTUALLY SETS THE EROSION'S DEPTH?
map_dc_diffusion.gd:602 derives strength from a depth in metres by
strength = DEPTH_METRES / (2 * EROSION_H), EROSION_H = 2625 m
i.e. it believes depth = 2 * strength * H. At the map's settings that
predicts 30 m. sweep_assumed_slope.py measured 0.23 m rms / 0.39 m peak on the
same field with the same parameters -- about 100x short. So the formula is
missing something.
THE SUSPECT is en_erosion_scale, the gully BAND. It shrank by 105x when the
band moved from 6300 m to 60 m (#55 section 3), and 105 is close to the 130x
shortfall. If depth tracks the band as well as the strength, the map's algebra
is wrong by exactly that factor and nobody would see erosion at a 60 m band.
THIS DOES NOT ASSUME THAT. It varies each one alone, over a wide range, and
reports the measured slope in log-log -- an exponent near 1 means proportional,
near 0 means the amplitude does not care about that parameter at all. The two
sweeps cross at the map's own operating point, which is the control: both must
report the same amplitude there or the harness is inconsistent.
python amplitude_law.py
"""
import argparse, os, subprocess
import numpy as np
from sweep_assumed_slope import (render, MAP_PARAMS, EROSION_H, GULLY_METRES,
DEPTH_METRES, TERRAIN_HEIGHT, CELL_SCALE)
MAP = dict(MAP_PARAMS)
MAP_STRENGTH = MAP["en_strength"]
MAP_ESCALE = MAP["en_erosion_scale"]
FACTORS = [1, 2, 4, 8, 16, 32, 64, 128]
def args_with(**over):
d = dict(MAP_PARAMS); d.update(over)
out = []
for k, v in d.items():
out += ["--set", k, str(v)]
return out
def amplitude(binpath, ox, oz, n, mpp, base, **over):
h = render(binpath, ox, oz, n, mpp, "/tmp/amp.f32", args_with(**over))
e = h - base
return float(np.sqrt((e**2).mean())), float(np.abs(e).max())
def loglog_slope(xs, ys):
xs = np.asarray(xs, float); ys = np.asarray(ys, float)
ok = (xs > 0) & (ys > 0)
if ok.sum() < 2: return float("nan")
return float(np.polyfit(np.log(xs[ok]), np.log(ys[ok]), 1)[0])
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--bin", default="/tmp/erosion_render")
ap.add_argument("--ox", type=float, default=50000.0)
ap.add_argument("--oz", type=float, default=34000.0)
ap.add_argument("--n", type=int, default=512)
ap.add_argument("--mpp", type=float, default=15.0)
a = ap.parse_args()
print(f"window {a.n*a.mpp/1000:.2f} km at ({a.ox:.0f}, {a.oz:.0f}) @ {a.mpp} m/px")
print(f"map operating point: strength {MAP_STRENGTH:.6g}, "
f"erosion_scale {MAP_ESCALE:.6g} (gully {GULLY_METRES:g} m)")
print(f"map PREDICTS depth = 2*strength*H = {2*MAP_STRENGTH*EROSION_H:.2f} m\n")
base = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/amp_off.f32",
args_with(en_strength=0))
print("A. STRENGTH alone (erosion_scale pinned at the map's)")
print(f"{'xmap':>6} {'strength':>12} {'rms m':>10} {'peak m':>10} {'2*s*H':>10}")
xs, ys = [], []
for f in FACTORS:
s = MAP_STRENGTH * f
rms, pk = amplitude(a.bin, a.ox, a.oz, a.n, a.mpp, base, en_strength=s)
print(f"{f:>6} {s:>12.6g} {rms:>10.3f} {pk:>10.3f} {2*s*EROSION_H:>10.2f}")
xs.append(f); ys.append(rms)
slope_s = loglog_slope(xs, ys)
print("\nB. EROSION_SCALE alone (strength pinned at the map's)")
print(f"{'xmap':>6} {'gully m':>10} {'esc':>12} {'rms m':>10} {'peak m':>10}")
xs2, ys2 = [], []
for f in FACTORS:
esc = MAP_ESCALE * f
rms, pk = amplitude(a.bin, a.ox, a.oz, a.n, a.mpp, base, en_erosion_scale=esc)
print(f"{f:>6} {GULLY_METRES*f:>10.0f} {esc:>12.6g} {rms:>10.3f} {pk:>10.3f}")
xs2.append(f); ys2.append(rms)
slope_e = loglog_slope(xs2, ys2)
print(f"\nlog-log slope vs strength : {slope_s:+.3f}"
" (1 = proportional, 0 = ignores it)")
print(f"log-log slope vs erosion_scale : {slope_e:+.3f}")
print(f"CONTROL, both sweeps at x1 : {ys[0]:.3f} vs {ys2[0]:.3f} m"
f" {'consistent' if abs(ys[0]-ys2[0]) < 1e-6 else 'INCONSISTENT — harness bug'}")
if ys[0] > 0:
need = DEPTH_METRES / ys[0]
print(f"\nTo reach the map's requested {DEPTH_METRES:g} m at the CURRENT band,"
f" amplitude must rise {need:.1f}x")
if abs(slope_s) > 0.5:
print(f" -> via strength: {MAP_STRENGTH:.6g} * {need**(1/slope_s):.1f}"
f" = {MAP_STRENGTH*need**(1/slope_s):.6g}")
if __name__ == "__main__":
main()