// =========================================================================
// hypso_curve.inc — the analytic land/sea transfer curve.
//
// ONE TEXT, compiled as GLSL (injected into the density shader) and as C++
// (through erosion_cpp_shim.h, for the oracle). Same rule as
// erosion_filter.inc: only the subset both languages agree on.
//
// WHAT IT IS FOR. A plain fBm has a UNIMODAL height distribution, so
// intersecting it with a water plane at any level gives a coastline of
// uniform fractal character everywhere -- no shelves, no basins, and no
// reason for land to clump. Earth's hypsography is BIMODAL: ~85% of the
// surface sits in two narrow bands, because continental crust and oceanic
// crust have different densities. This curve puts that shape into the
// height field directly, with no real-world data.
//
// WHY A BLEND AND NOT A PIECEWISE JOIN. The obvious construction --
// `h < sea ? ocean(h) : land(h)` -- is C0 but not C1: the slope jumps at
// the join, which is EXACTLY the waterline, and a slope discontinuity along
// every coastline renders as a crease ring. Blending two regimes with a
// smoothstep is C1 by construction because smoothstep is.
//
// WHERE THE BIMODALITY COMES FROM. Not from a threshold. Output values
// cluster where the curve is FLAT (low |dhv/dh|), so a flat abyssal
// plateau and a land ramp separated by a narrow STEEP section produce two
// modes with few samples between them -- the continental slope. `shelf` is
// the half-width of that steep section: smaller = sharper coast = more
// strongly bimodal.
//
// THE WATERLINE IS DERIVED. The curve is monotonic, so sea level in OUTPUT
// space is simply hypso_curve(sea).x. Do not also hand-set a water height;
// read it from here, or the two disagree and the shoreline renders where
// the terrain is not.
//
// PEAK LIFT is additive on the upper tail only, so it raises summits
// without steepening mid-slopes or deepening valleys. That matters beyond
// looks: total relief is bounded by the physics band (character movement
// degrades with |y|), and lifting only the tail is how peak height stops
// competing with valley depth for the same budget.
//
// h input height, nominally [-1, 1]
// sea input value at the middle of the continental slope
// shelf half-width of that slope, in input units (>0)
// abyss OUTPUT value of the abyssal plateau (negative)
// land_slope output units per input unit above sea
// peak_onset input value where the summit lift begins
// peak_gain output units added at h = 1
//
// Returns (hv, dhv/dh). The derivative is not decoration: the erosion
// filter is slope-aware and the wrapper hands it a gradient, so remapping
// the height without remapping the gradient by the chain rule would carve
// gullies aligned to a surface that is no longer there.
// =========================================================================
// WHY THREE REGIMES AND NOT TWO. A first version had an abyssal plateau and
// a linear land ramp, and the gate measured its output histogram as ONE mode
// plus a smear -- because a ramp has no mode. Earth's bimodality is two
// PLATEAUS: the continental platform (~0..+1 km) and the abyssal plain
// (~-4..-5 km), with the continental slope between them. So the land side
// needs a low-gradient shelf of its own, and the mountains rise ABOVE that.
//
// abyss --[continental slope]-- platform --[range front]-- mountains
// flat steep gentle steeper + peak lift
//
// Two smoothsteps composing three lines. Both are C1, so the whole is.
// smoothstep's derivative. GLSL has smoothstep but not its derivative, and
// writing 6t(1-t)/(hi-lo) at three call sites is how one of them ends up
// missing the 1/(hi-lo).
// ⚠ DO NOT NAME A LOCAL `out`. It is a GLSL storage qualifier and a perfectly
// ordinary identifier in C++, so a C++-only test compiles it and the SHADER
// dies with "syntax error, unexpected OUT" -- which is what happened: the GPU
// fell back to the CPU mesher and the only symptom was that terrain generation
// got slow. Same class as erosion_filter.inc's PN_TAU note. `in`, `inout`,
// `filter`, `sample` and `active` are reserved the same way.
float hypso_dsmooth(float lo, float hi, float x) {
float d = hi - lo;
if (d <= 0.0) return 0.0;
float t = clamp((x - lo) / d, 0.0, 1.0);
return 6.0 * t * (1.0 - t) / d;
}
// h input height, nominally [-1, 1]
// sea input value that IS sea level; the curve returns 0 there
// shelf_width how far below sea the gentle shelf runs before dropping
// slope_width half-width of the steep drop (smaller = sharper break)
// abyss OUTPUT value of the abyssal plain (negative)
// plat_slope output per input on the continental platform (small)
// mtn_onset input value where the mountain front begins
// mtn_width half-width of that front
// mtn_slope output per input in the mountains (larger)
// peak_onset input value where the summit lift begins
// peak_gain output added at h = 1
//
// Returns (hv, dhv/dh). The derivative is not decoration: the erosion
// filter is slope-aware and the wrapper hands it a gradient, so remapping
// the height without remapping the gradient by the chain rule would carve
// gullies aligned to a surface that is no longer there.
vec2 hypso_curve(float h, float sea, float shelf_width, float slope_width,
float abyss,
float plat_slope, float mtn_onset, float mtn_width,
float mtn_slope, float peak_onset, float peak_gain) {
// --- regime 1: the continental platform, a gentle line through (sea, 0)
float plat = plat_slope * (h - sea);
// --- the continental slope, CENTRED BELOW SEA LEVEL.
//
// ⚠ THIS IS WHY THE FIRST VERSION WAS A CLIFF. The blend was centred ON
// `sea`, so the steep section straddled the waterline and land met deep
// water with nothing in between -- a wall, roughly 50x steeper than the
// platform behind it. That steepness IS the bimodality, so it cannot be
// removed; it has to be MOVED.
//
// On Earth the steep part is submerged: land -> beach -> a gentle shelf
// ~100 m deep and tens of km wide -> shelf break -> continental slope ->
// abyssal plain. Putting the blend `shelf_width` BELOW sea leaves the
// gentle platform line running down through the waterline, and that
// stretch is the continental shelf.
//
// Requires shelf_width > slope_width, or the drop reaches back up through
// sea level and the shelf does not exist.
float slope_mid = sea - shelf_width;
float c0 = slope_mid - slope_width;
float c1 = slope_mid + slope_width;
float w = smoothstep(c0, c1, h);
float dw = hypso_dsmooth(c0, c1, h);
float base = abyss + (plat - abyss) * w;
float dbase = (plat - abyss) * dw + w * plat_slope;
// --- the mountain front, ADDED not blended to. See the note above:
// blending toward a steeper LINE folds the terrain below its onset,
// measured at dhv/dh = -0.112. An additive ramp cannot fold.
float a = mtn_onset - mtn_width;
float b = mtn_onset + mtn_width;
float sm = smoothstep(a, b, h);
float dsm = hypso_dsmooth(a, b, h);
// (h-a)*smoothstep: exactly 0 below a (smoothstep clamps), asymptotically
// linear with unit slope above b, C1 throughout, never decreasing.
float ramp = (h - a) * sm;
float dramp = sm + (h - a) * dsm;
float extra = mtn_slope - plat_slope; // the slope the range ADDS
float hv = base + extra * ramp;
float dhv = dbase + extra * dramp;
// --- summit lift, upper tail only. Same additive shape, same reason.
if (peak_gain != 0.0) {
hv += peak_gain * smoothstep(peak_onset, 1.0, h);
dhv += peak_gain * hypso_dsmooth(peak_onset, 1.0, h);
}
return vec2(hv, dhv);
}