From e71f145b71de168508fb6b4bdda5df7e7f69cad3 Mon Sep 17 00:00:00 2001
From: Joshua Yanovski <pythonesque@gmail.com>
Date: Tue, 3 Dec 2019 02:07:44 +0100
Subject: [PATCH] Sediment transport, plus many other things.

---
 server/src/lib.rs          |   7 +-
 world/examples/view.rs     |   7 +-
 world/examples/water.rs    | 195 +++++++++---
 world/src/block/mod.rs     |  14 +-
 world/src/column/mod.rs    |  98 ++++--
 world/src/config.rs        |   8 +-
 world/src/lib.rs           |   4 +-
 world/src/sim/diffusion.rs |  76 ++++-
 world/src/sim/erosion.rs   | 606 ++++++++++++++++++++++++++-----------
 world/src/sim/mod.rs       | 286 +++++++++++++----
 world/src/sim/util.rs      |  22 +-
 11 files changed, 995 insertions(+), 328 deletions(-)

diff --git a/server/src/lib.rs b/server/src/lib.rs
index d71c67b025..9abc4c90dd 100644
--- a/server/src/lib.rs
+++ b/server/src/lib.rs
@@ -46,7 +46,7 @@ use std::{
 };
 use uvth::{ThreadPool, ThreadPoolBuilder};
 use vek::*;
-use world::{sim::WORLD_SIZE, World};
+use world::{sim::{WORLD_SIZE, WorldOpts}, World};
 const CLIENT_TIMEOUT: f64 = 20.0; // Seconds
 
 pub enum Event {
@@ -104,7 +104,10 @@ impl Server {
         state.ecs_mut().register::<RegionSubscription>();
         state.ecs_mut().register::<Client>();
 
-        let world = World::generate(settings.world_seed);
+        let world = World::generate(settings.world_seed, WorldOpts {
+            seed_elements: true,
+            ..WorldOpts::default()
+        });
 
         let spawn_point = {
             // NOTE: all of these `.map(|e| e as [type])` calls should compile into no-ops,
diff --git a/world/examples/view.rs b/world/examples/view.rs
index c1f52ec8e6..eab8ef1398 100644
--- a/world/examples/view.rs
+++ b/world/examples/view.rs
@@ -1,12 +1,15 @@
 use std::ops::{Add, Mul, Sub};
 use vek::*;
-use veloren_world::{util::Sampler, World};
+use veloren_world::{sim::WorldOpts, util::Sampler, World};
 
 const W: usize = 640;
 const H: usize = 480;
 
 fn main() {
-    let world = World::generate(0);
+    let world = World::generate(0, WorldOpts {
+        seed_elements: true,
+        ..WorldOpts::default()
+    });
 
     let sampler = world.sample_columns();
 
diff --git a/world/examples/water.rs b/world/examples/water.rs
index 72900e489e..818262f183 100644
--- a/world/examples/water.rs
+++ b/world/examples/water.rs
@@ -1,37 +1,70 @@
 use common::{terrain::TerrainChunkSize, vol::RectVolSize};
-use std::f32;
+// use self::Mode::*;
+use std::{f32, f64};
 use vek::*;
 use veloren_world::{
-    sim::{RiverKind, WORLD_SIZE},
+    sim::{RiverKind, WorldOpts, WORLD_SIZE},
+    util::Sampler,
     World, CONFIG,
 };
 
 const W: usize = 1024;
 const H: usize = 1024;
 
+/* enum Mode {
+    /// Directional keys affect position of the camera.
+    ///
+    /// (W A S D move left and right, F B zoom in and out).
+    Alt,
+    /// Directional keys affect angle of the lens
+    ///
+    /// (W
+    Lens,
+    /// Directional keys affect light direction.
+    ///
+    /// (W A S D move left and right, F B move towards and awaay).
+    Light,
+}; */
+
 fn main() {
     pretty_env_logger::init();
 
-    let world = World::generate(1337);
+    let world = World::generate(1337, WorldOpts {
+        seed_elements: false,
+        ..WorldOpts::default()
+    });
 
     let sampler = world.sim();
 
     let mut win =
         minifb::Window::new("World Viewer", W, H, minifb::WindowOptions::default()).unwrap();
 
-    let mut focus = Vec2::zero();
+    let mut focus = Vec3::new(0.0, 0.0, CONFIG.sea_level as f64);
+    // Altitude is divided by gain and clamped to [0, 1]; thus, decreasing gain makes
+    // smaller differences in altitude appear larger.
     let mut gain = CONFIG.mountain_scale;
-    let mut scale = (WORLD_SIZE.x / W) as i32;
+    // The Z component during normal calculations is multiplied by gain; thus,
+    let mut lgain = 1.0;
+    let mut scale = (WORLD_SIZE.x / W) as f64;
 
-    let light_direction = Vec3::new(-0.8, -1.0, 0.3).normalized();
+    // Right-handed coordinate system: light is going left, down, and "backwards" (i.e. on the
+    // map, where we translate the y coordinate on the world map to z in the coordinate system,
+    // the light comes from -y on the map and points towards +y on the map).  In a right
+    // handed coordinate system, the "camera" points towards -z, so positive z is backwards
+    // "into" the camera.
+    //
+    // "In world space the x-axis will be pointing east, the y-axis up and the z-axis will be pointing south"
+    let mut light_direction = Vec3::new(-0.8, -1.0, 0.3);
     let light_res = 3;
 
     let mut is_basement = false;
+    let mut is_water = true;
     let mut is_shaded = true;
     let mut is_temperature = true;
     let mut is_humidity = true;
 
     while win.is_open() {
+        let light = light_direction.normalized();
         let mut buf = vec![0; W * H];
         const QUADRANTS: usize = 4;
         let mut quads = [[0u32; QUADRANTS]; QUADRANTS];
@@ -40,10 +73,12 @@ fn main() {
         let mut oceans = 0u32;
 
         // let water_light = (light_direction.z + 1.0) / 2.0 * 0.8 + 0.2;
+        let focus_rect = Vec2::from(focus);
+        let true_sea_level = (CONFIG.sea_level as f64 - focus.z) / gain as f64;
 
         for i in 0..W {
             for j in 0..H {
-                let pos = focus + Vec2::new(i as i32, j as i32) * scale;
+                let pos = (focus_rect + Vec2::new(i as f64, j as f64) * scale).map(|e: f64| e as i32);
                 /* let top_left = pos;
                 let top_right = focus + Vec2::new(i as i32 + light_res, j as i32) * scale;
                 let bottom_left = focus + Vec2::new(i as i32, j as i32 + light_res) * scale; */
@@ -90,26 +125,53 @@ fn main() {
                     .get(cross_pos)
                     .map(|s| if is_basement { s.basement } else { s.alt })
                     .unwrap_or(CONFIG.sea_level);
+                // Pointing downhill, forward
+                // (index--note that (0,0,1) is backward right-handed)
                 let forward_vec = Vec3::new(
                     (downhill_pos.x - pos.x) as f64,
+                    (downhill_alt - alt) as f64 * lgain,
                     (downhill_pos.y - pos.y) as f64,
-                    (downhill_alt - alt) as f64 * (CONFIG.mountain_scale as f64 / gain as f64),
                 );
+                // Pointing 90 degrees left (in horizontal xy) of downhill, up
+                // (middle--note that (1,0,0), 90 degrees CCW backward, is right right-handed)
                 let up_vec = Vec3::new(
                     (cross_pos.x - pos.x) as f64,
+                    (cross_alt - alt) as f64 * lgain,
                     (cross_pos.y - pos.y) as f64,
-                    (cross_alt - alt) as f64 * (CONFIG.mountain_scale as f64 / gain as f64),
                 );
+                // Then cross points "to the right" (upwards) on a right-handed coordinate system.
+                // (right-handed coordinate system means (0, 0, 1.0) is "forward" into the screen).
                 let surface_normal = forward_vec.cross(up_vec).normalized();
+                // f = (0, alt_bl - alt_tl, 1) [backward right-handed = (0,0,1)]
+                // u = (1, alt_tr - alt_tl, 0) [right (90 degrees CCW backward) = (1,0,0)]
+                // (f × u in right-handed coordinate system: pointing up)
+                //
+                // f × u =
+                //   (a.y*b.z - a.z*b.y,
+                //    a.z*b.x - a.x*b.z,
+                //    a.x*b.y - a.y*b.x,
+                //   )
+                // =
+                //   (-(alt_tr - alt_tl),
+                //    1,
+                //    -(alt_bl - alt_tl),
+                //   )
+                // =
+                //   (alt_tl - alt_tr,
+                //    1,
+                //    alt_tl - alt_bl,
+                //   )
+                //
                 // let surface_normal = Vec3::new((alt_tl - alt_tr) as f64, 1.0, (alt_tl - alt_bl) as f64).normalized();
-                let light = (surface_normal.dot(light_direction) + 1.0) / 2.0;
-                let light = (light * 0.8) + 0.2;
+                let light = (surface_normal.dot(light) + 1.0) / 2.0;
+                let light = (light * 0.9) + 0.1;
 
-                let water_alt = ((alt.max(water_alt) - CONFIG.sea_level) as f64 / gain as f64)
+                let true_water_alt = (alt.max(water_alt) as f64 - focus.z) / gain as f64;
+                let true_alt = (alt as f64 - focus.z) / gain as f64;
+                let water_depth = (true_water_alt - true_alt)
                     .min(1.0)
                     .max(0.0);
-                let true_alt = (alt - CONFIG.sea_level) as f64 / gain as f64;
-                let water_depth = (water_alt - true_alt)
+                let water_alt = true_water_alt
                     .min(1.0)
                     .max(0.0);
                 let alt = true_alt
@@ -133,20 +195,8 @@ fn main() {
                     None => {}
                 }
 
-                buf[j * W + i] = match river_kind {
-                    Some(RiverKind::Ocean) => u32::from_le_bytes([
-                        ((64.0 - water_depth * 64.0) * 1.0) as u8,
-                        ((32.0 - water_depth * 32.0) * 1.0) as u8,
-                        0,
-                        255,
-                    ]),
-                    Some(RiverKind::River { .. }) => u32::from_le_bytes([
-                        64 + (alt * 191.0) as u8,
-                        32 + (alt * 95.0) as u8,
-                        0,
-                        255,
-                    ]),
-                    None if water_alt >= water_depth => {
+                buf[j * W + i] = match (river_kind, (is_water, true_alt >= true_sea_level)) {
+                    (_, (false, _)) | ( None, (_, true)) => {
                         let (r, g, b) = (
                             (if is_shaded { alt } else { alt } * if is_temperature { temperature as f64 } else if is_shaded { alt } else { 0.0 }).sqrt(),
                             if is_shaded { 0.2 + (alt * 0.8) } else { alt },
@@ -170,7 +220,19 @@ fn main() {
                             255,
                         ]) */
                     },
-                    None | Some(RiverKind::Lake { .. }) => u32::from_le_bytes([
+                    (Some(RiverKind::Ocean), _) => u32::from_le_bytes([
+                        ((64.0 - water_depth * 64.0) * 1.0) as u8,
+                        ((32.0 - water_depth * 32.0) * 1.0) as u8,
+                        0,
+                        255,
+                    ]),
+                    (Some(RiverKind::River { .. }), _) => u32::from_le_bytes([
+                        64 + (alt * 191.0) as u8,
+                        32 + (alt * 95.0) as u8,
+                        0,
+                        255,
+                    ]),
+                    (None, _) | (Some(RiverKind::Lake { .. }), _) => u32::from_le_bytes([
                         (((64.0 + water_alt * 191.0) + (- water_depth * 64.0)) * 1.0) as u8,
                         (((32.0 + water_alt * 95.0) + (- water_depth * 32.0)) * 1.0) as u8,
                         0,
@@ -180,12 +242,14 @@ fn main() {
             }
         }
 
-        let spd = 32;
+        let spd = 32.0;
+        let lspd = 0.1;
         if win.is_key_down(minifb::Key::P) {
             println!(
                 "\
-                 Gain: {:?}\n\
-                 Scale: {:?}\n\
+                 Gain / Shade gain: {:?} / {:?}\n\
+                 Scale / Focus: {:?} / {:?}\n\
+                 Light: {:?}
                  Land(adjacent): (X = temp, Y = humidity): {:?}\n\
                  Rivers: {:?}\n\
                  Lakes: {:?}\n\
@@ -193,7 +257,10 @@ fn main() {
                  Total water: {:?}\n\
                  Total land(adjacent): {:?}",
                 gain,
+                lgain,
                 scale,
+                focus,
+                light_direction,
                 quads,
                 rivers,
                 lakes,
@@ -204,13 +271,14 @@ fn main() {
         }
         if win.get_mouse_down(minifb::MouseButton::Left) {
             if let Some((mx, my)) = win.get_mouse_pos(minifb::MouseMode::Clamp) {
-                let pos = focus + Vec2::new(mx as i32, my as i32) * scale;
+                let pos = (focus_rect + (Vec2::new(mx as f64, my as f64) * scale)).map(|e| e as i32);
                 println!(
                     "Chunk position: {:?}",
                     pos.map2(TerrainChunkSize::RECT_SIZE, |e, f| e * f as i32)
                 );
             }
         }
+        let is_camera = win.is_key_down(minifb::Key::C);
         if win.is_key_down(minifb::Key::B) {
             is_basement ^= true;
         }
@@ -220,32 +288,77 @@ fn main() {
         if win.is_key_down(minifb::Key::T) {
             is_temperature ^= true;
         }
+        if win.is_key_down(minifb::Key::O) {
+            is_water ^= true;
+        }
         if win.is_key_down(minifb::Key::L) {
             is_shaded ^= true;
         }
         if win.is_key_down(minifb::Key::W) {
-            focus.y -= spd * scale;
+            if is_camera {
+                light_direction.z -= lspd;
+            } else {
+                focus.y -= spd * scale;
+            }
         }
         if win.is_key_down(minifb::Key::A) {
-            focus.x -= spd * scale;
+            if is_camera {
+                light_direction.x -= lspd;
+            } else {
+                focus.x -= spd * scale;
+            }
         }
         if win.is_key_down(minifb::Key::S) {
-            focus.y += spd * scale;
+            if is_camera {
+                light_direction.z += lspd;
+            } else {
+                focus.y += spd * scale;
+            }
         }
         if win.is_key_down(minifb::Key::D) {
-            focus.x += spd * scale;
+            if is_camera {
+                light_direction.x += lspd;
+            } else {
+                focus.x += spd * scale;
+            }
         }
         if win.is_key_down(minifb::Key::Q) {
-            gain += 64.0;
+            if is_camera {
+                if (lgain * 2.0).is_normal() {
+                    lgain *= 2.0;
+                }
+            } else {
+                gain += 64.0;
+            }
         }
         if win.is_key_down(minifb::Key::E) {
-            gain = (gain - 64.0).max(64.0);
+            if is_camera {
+                if (lgain / 2.0).is_normal() {
+                    lgain /= 2.0;
+                }
+            } else {
+                gain = (gain - 64.0).max(64.0);
+            }
         }
         if win.is_key_down(minifb::Key::R) {
-            scale += 1;
+            if is_camera {
+                focus.z += spd * scale;
+            } else {
+                if (scale * 2.0).is_normal() {
+                    scale *= 2.0;
+                }
+                // scale += 1;
+            }
         }
         if win.is_key_down(minifb::Key::F) {
-            scale = (scale - 1).max(0);
+            if is_camera {
+                focus.z -= spd * scale;
+            } else {
+                if (scale / 2.0).is_normal() {
+                    scale /= 2.0;
+                }
+                // scale = (scale - 1).max(0);
+            }
         }
 
         win.update_with_buffer(&buf).unwrap();
diff --git a/world/src/block/mod.rs b/world/src/block/mod.rs
index 17d22bf217..de87f2d28f 100644
--- a/world/src/block/mod.rs
+++ b/world/src/block/mod.rs
@@ -149,6 +149,7 @@ impl<'a> BlockGen<'a> {
         let sample = &z_cache?.sample;
         let &ColumnSample {
             alt,
+            basement,
             chaos,
             water_level,
             warp_factor,
@@ -177,10 +178,10 @@ impl<'a> BlockGen<'a> {
         let wposf = wpos.map(|e| e as f64);
 
         let (block, height) = if !only_structures {
-            let (_definitely_underground, height, on_cliff, water_height) =
+            let (_definitely_underground, height, on_cliff, basement_height, water_height) =
                 if (wposf.z as f32) < alt - 64.0 * chaos {
                     // Shortcut warping
-                    (true, alt, false, water_level)
+                    (true, alt, false, basement, water_level)
                 } else {
                     // Apply warping
                     let warp = world
@@ -223,6 +224,7 @@ impl<'a> BlockGen<'a> {
                         false,
                         height,
                         on_cliff,
+                        basement.min(basement + warp),
                         (if water_level <= alt {
                             water_level + warp
                         } else {
@@ -247,10 +249,11 @@ impl<'a> BlockGen<'a> {
 
             let water = Block::new(BlockKind::Water, Rgb::new(60, 90, 190));
 
-            let grass_depth = 1.5 + 2.0 * chaos;
+            let grass_depth = (1.5 + 2.0 * chaos).min(height - basement_height);
             let block = if (wposf.z as f32) < height - grass_depth {
                 let col = Lerp::lerp(
-                    saturate_srgb(sub_surface_color, 0.45).map(|e| (e * 255.0) as u8),
+                    // saturate_srgb(sub_surface_color, 0.45).map(|e| (e * 255.0) as u8),
+                    sub_surface_color.map(|e| (e * 255.0) as u8),
                     stone_col,
                     (height - grass_depth - wposf.z as f32) * 0.15,
                 );
@@ -272,7 +275,8 @@ impl<'a> BlockGen<'a> {
                 // Surface
                 Some(Block::new(
                     BlockKind::Normal,
-                    saturate_srgb(col, 0.45).map(|e| (e * 255.0) as u8),
+                    // saturate_srgb(col, 0.45).map(|e| (e * 255.0) as u8),
+                    col.map(|e| (e * 255.0) as u8),
                 ))
             } else if (wposf.z as f32) < height + 0.9
                 && temp < CONFIG.desert_temp
diff --git a/world/src/column/mod.rs b/world/src/column/mod.rs
index eb40894805..5d6719f938 100644
--- a/world/src/column/mod.rs
+++ b/world/src/column/mod.rs
@@ -443,7 +443,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
             sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
             // sim.get_interpolated(wpos, |chunk| chunk.alt)?
         };
-        let basement = sim./*get_interpolated*/get_interpolated_monotone(wpos, |chunk| chunk.basement)?;
+        let basement = alt + sim./*get_interpolated*/get_interpolated_monotone(wpos, |chunk| chunk.basement.sub(chunk.alt))?;
 
         // Find the average distance to each neighboring body of water.
         let mut river_count = 0.0f64;
@@ -755,7 +755,8 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
                                             lake_dist <= TerrainChunkSize::RECT_SIZE.x as f64 * 0.5;
                                         if gouge_factor == 1.0 {
                                             return Some((
-                                                alt_for_river < lake_water_alt || in_bounds,
+                                                /*alt_for_river < lake_water_alt || in_bounds,*/
+                                                true,
                                                 alt.min(lake_water_alt - 1.0 - river_gouge),
                                                 downhill_water_alt.max(lake_water_alt)
                                                     - river_gouge,
@@ -763,7 +764,8 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
                                             ));
                                         } else {
                                             return Some((
-                                                alt_for_river < lake_water_alt || in_bounds,
+                                                /*alt_for_river < lake_water_alt || in_bounds,*/
+                                                true,
                                                 alt_for_river,
                                                 if in_bounds_ {
                                                     downhill_water_alt.max(lake_water_alt)
@@ -804,6 +806,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
         } else {
             (false, alt_for_river, downhill_water_alt, 1.0)
         };
+        let warp_factor = 0.0;
 
         let riverless_alt_delta = (sim
             .gen_ctx
@@ -867,6 +870,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
         // let wet_grass = Rgb::new(0.1, 0.5, 0.5);
         let cold_stone = Rgb::new(0.57, 0.67, 0.8);
         // let cold_stone = Rgb::new(0.5, 0.5, 0.5);
+        let hot_stone = Rgb::new(0.07, 0.07, 0.06);
         let warm_stone = Rgb::new(0.77, 0.77, 0.64);
         // //let warm_stone = Rgb::new(0.6, 0.6, 0.5);
         // let warm_stone = Rgb::new(0.6, 0.5, 0.1);
@@ -890,7 +894,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
         );
         let tundra = Lerp::lerp(snow, Rgb::new(0.01, 0.3, 0.0), 0.4 + marble * 0.6);
         let dead_tundra = Lerp::lerp(warm_stone, Rgb::new(0.3, 0.12, 0.2), marble);
-        let cliff = Rgb::lerp(cold_stone, warm_stone, marble);
+        let cliff = Rgb::lerp(cold_stone, /*warm_stone*/hot_stone, marble);
 
         let grass = Rgb::lerp(
             cold_grass,
@@ -918,7 +922,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
 
         // For below desert humidity, we are always sand or rock, depending on altitude and
         // temperature.
-        let ground = Rgb::lerp(
+        /* let ground = Rgb::lerp(
             cliff,
             Rgb::lerp(
                 dead_tundra,
@@ -932,7 +936,39 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
             /* alt.sub(CONFIG.sea_level)
                 .sub(CONFIG.mountain_scale * 0.25)
                 .div(CONFIG.mountain_scale * 0.125), */
+        ); */
+        let ground = Lerp::lerp(
+            Lerp::lerp(
+                dead_tundra,
+                sand,
+                temp.sub(CONFIG.snow_temp)
+                    .div(CONFIG.desert_temp.sub(CONFIG.snow_temp))
+                    .mul(/*4.5*/ 0.5),
+            ),
+            dirt,
+            humidity
+                .sub(CONFIG.desert_hum)
+                .div(CONFIG.forest_hum.sub(CONFIG.desert_hum))
+                .mul(1.0),
         );
+
+        let sub_surface_color = Lerp::lerp(
+            cliff,
+            ground,
+            alt.sub(basement).mul(0.25)
+        );
+
+        /*  let ground = Rgb::lerp(
+            dead_tundra,
+            sand,
+            temp.sub(CONFIG.snow_temp)
+                .div(CONFIG.desert_temp.sub(CONFIG.snow_temp))
+                .mul(/*4.5*/ 0.5),
+            /* alt.sub(CONFIG.sea_level)
+                .sub(CONFIG.mountain_scale * 0.25)
+                .div(CONFIG.mountain_scale * 0.125), */
+        ); */
+
         // From desert to forest humidity, we go from tundra to dirt to grass to moss to sand,
         // depending on temperature.
         let ground = Rgb::lerp(
@@ -942,18 +978,20 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
                     Rgb::lerp(
                         Rgb::lerp(
                             tundra,
-                            // snow_temp to 0
+                            // snow_temp to temperate_temp
                             dirt,
                             temp.sub(CONFIG.snow_temp)
-                                .div(CONFIG.snow_temp.neg())
+                                .div(CONFIG.temperate_temp.sub(CONFIG.snow_temp))
                                 /*.sub((marble - 0.5) * 0.05)
                                 .mul(256.0)*/
                                 .mul(1.0),
                             // .mul(2.0),
                         ),
-                        // 0 to tropical_temp
+                        // temperate_temp to tropical_temp
                         grass,
-                        temp.div(CONFIG.tropical_temp).mul(4.0),
+                        temp.sub(CONFIG.temperate_temp)
+                            .div(CONFIG.tropical_temp.sub(CONFIG.temperate_temp))
+                            .mul(4.0),
                     ),
                     // tropical_temp to desert_temp
                     moss,
@@ -983,9 +1021,11 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
                 Rgb::lerp(
                     Rgb::lerp(
                         snow_moss,
-                        // 0 to tropical_temp
+                        // temperate_temp to tropical_temp
                         grass,
-                        temp.div(CONFIG.tropical_temp).mul(4.0),
+                        temp.sub(CONFIG.temperate_temp)
+                            .div(CONFIG.tropical_temp.sub(CONFIG.temperate_temp))
+                            .mul(4.0),
                     ),
                     // tropical_temp to desert_temp
                     tropical,
@@ -1014,9 +1054,11 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
                 Rgb::lerp(
                     Rgb::lerp(
                         snow_moss,
-                        // 0 to tropical_temp
+                        // temperate_temp to tropical_temp
                         rainforest,
-                        temp.div(CONFIG.tropical_temp).mul(4.0),
+                        temp.sub(CONFIG.temperate_temp)
+                            .div(CONFIG.tropical_temp.sub(CONFIG.temperate_temp))
+                            .mul(4.0),
                     ),
                     // tropical_temp to desert_temp
                     tropical,
@@ -1035,15 +1077,35 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
             humidity.sub(CONFIG.jungle_hum).mul(1.0),
         );
 
-        // Snow covering
+        /* // Bedrock
         let ground = Rgb::lerp(
+            cliff,
+            ground,
+            alt.sub(basement)
+                .mul(0.25)
+        ); */
+
+        // Snow covering
+        let snow_cover =
+            temp.sub(CONFIG.snow_temp)
+                .max(-humidity.sub(CONFIG.desert_hum))
+                .mul(16.0)
+                .add((marble_small - 0.5) * 0.5);
+        let (alt, ground, sub_surface_color) = if snow_cover /*< 0.1*/<= 0.5 && alt > water_level {
+            // Allow snow cover.
+            (alt + 1.0 - snow_cover.max(0.0), Rgb::lerp(snow, ground, snow_cover),
+                Lerp::lerp(sub_surface_color, ground, alt.sub(basement).mul(0.15)))
+        } else {
+            (alt, ground, sub_surface_color)
+        };
+        /* let ground = Rgb::lerp(
             snow,
             ground,
             temp.sub(CONFIG.snow_temp)
                 .max(-humidity.sub(CONFIG.desert_hum))
                 .mul(16.0)
                 .add((marble_small - 0.5) * 0.5),
-        );
+        ); */
 
         // Caves
         let cave_at = |wposf: Vec2<f64>| {
@@ -1091,6 +1153,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
 
         Some(ColumnSample {
             alt,
+            basement,
             chaos,
             water_level,
             warp_factor,
@@ -1101,9 +1164,9 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
                 // Beach
                 ((ocean_level - 1.0) / 2.0).max(0.0),
             ),
-            sub_surface_color: /*warm_grass*/Lerp::lerp(cliff, dirt, alt.sub(basement).mul(0.25)),
+            sub_surface_color,// /*warm_grass*/Lerp::lerp(cliff, dirt, alt.sub(basement).mul(0.25)),
             // No growing directly on bedrock.
-            tree_density: if alt - basement < 2.0 { 0.0 } else { tree_density },
+            tree_density: Lerp::lerp(0.0, tree_density, alt.sub(2.0).sub(basement).mul(0.5)),
             forest_kind: sim_chunk.forest_kind,
             close_structures: self.gen_close_structures(wpos),
             cave_xy,
@@ -1139,6 +1202,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
 #[derive(Clone)]
 pub struct ColumnSample<'a> {
     pub alt: f32,
+    pub basement: f32,
     pub chaos: f32,
     pub water_level: f32,
     pub warp_factor: f32,
diff --git a/world/src/config.rs b/world/src/config.rs
index 1e5a372409..d14325521d 100644
--- a/world/src/config.rs
+++ b/world/src/config.rs
@@ -2,6 +2,7 @@ pub struct Config {
     pub sea_level: f32,
     pub mountain_scale: f32,
     pub snow_temp: f32,
+    pub temperate_temp: f32,
     pub tropical_temp: f32,
     pub desert_temp: f32,
     pub desert_hum: f32,
@@ -46,9 +47,10 @@ pub struct Config {
 pub const CONFIG: Config = Config {
     sea_level: 140.0,
     mountain_scale: 2048.0,
-    snow_temp: -0.6,
-    tropical_temp: 0.2,
-    desert_temp: 0.6,
+    snow_temp: -0.8,
+    temperate_temp: -0.4,
+    tropical_temp: 0.4,
+    desert_temp: 0.8,
     desert_hum: 0.15,
     forest_hum: 0.5,
     jungle_hum: 0.85,
diff --git a/world/src/lib.rs b/world/src/lib.rs
index cf68981989..967111076b 100644
--- a/world/src/lib.rs
+++ b/world/src/lib.rs
@@ -36,9 +36,9 @@ pub struct World {
 }
 
 impl World {
-    pub fn generate(seed: u32) -> Self {
+    pub fn generate(seed: u32, opts: sim::WorldOpts) -> Self {
         Self {
-            sim: sim::WorldSim::generate(seed),
+            sim: sim::WorldSim::generate(seed, opts),
         }
     }
 
diff --git a/world/src/sim/diffusion.rs b/world/src/sim/diffusion.rs
index 08d2373556..30a9503f0f 100644
--- a/world/src/sim/diffusion.rs
+++ b/world/src/sim/diffusion.rs
@@ -140,10 +140,10 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
     enddo
 */
         for i in 1..nx-1 {
-            factxp = (kdint[aij(i+1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factxm = (kdint[aij(i-1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factyp = (kdint[aij(i, j+1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
-            factym = (kdint[aij(i, j-1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factxp = (kdint[aij(i+1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxm = (kdint[aij(i-1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factyp = (kdint[aij(i, j+1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factym = (kdint[aij(i, j-1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
             diag[i] = 1.0 + factxp + factxm;
             sup[i] = -factxp;
             inf[i] = -factxm;
@@ -167,12 +167,21 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
     f(1)=zintp(1,j)+factyp*zintp(1,j+1)-(factyp+factym)*zintp(1,j)+factym*zintp(1,j-1)
     endif
 */
-        if true {
+        if false {
             diag[0] = 1.0;
             sup[0] = 0.0;
             f[0] = zintp[aij(0, j)];
         } else {
-            // FIXME: implement reflective boundary
+            // reflective boundary
+            factxp = (kdint[aij(1, j)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factyp = (kdint[aij(0, j+1)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factym = (kdint[aij(0, j-1)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            diag[0] = 1.0 + factxp;
+            sup[0] = -factxp;
+            f[0] =
+                zintp[aij(0, j)] + factyp * zintp[aij(0, j+1)] -
+                (factyp + factym) *
+                zintp[aij(0, j)] + factym * zintp[aij(0, j-1)];
         }
 /*
 ! right bc
@@ -189,12 +198,21 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
     f(nx)=zintp(nx,j)+factyp*zintp(nx,j+1)-(factyp+factym)*zintp(nx,j)+factym*zintp(nx,j-1)
     endif
 */
-        if true {
+        if false {
             diag[nx - 1] = 1.0;
             inf[nx - 1] = 0.0;
             f[nx - 1] = zintp[aij(nx - 1, j)];
         } else {
-            // FIXME: implement reflective boundary
+            // reflective boundary
+            factxm = (kdint[aij(nx-2, j)] + kdint[aij(nx-1, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factyp = (kdint[aij(nx-1, j+1)] + kdint[aij(nx-1, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factym = (kdint[aij(nx-1, j-1)] + kdint[aij(nx-1, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            diag[nx - 1] = 1.0 + factxm;
+            inf[nx - 1] = -factxm;
+            f[nx - 1] =
+                zintp[aij(nx-1, j)] + factyp * zintp[aij(nx-1, j+1)] -
+                (factyp + factym) *
+                zintp[aij(nx-1, j)] + factym * zintp[aij(nx-1, j-1)];
         }
 /*
   call tridag (inf,diag,sup,f,res,nx)
@@ -242,10 +260,10 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
     enddo
 */
         for j in 1..ny-1 {
-            factxp = (kdint[aij(i+1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factxm = (kdint[aij(i-1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
-            factyp = (kdint[aij(i, j+1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
-            factym = (kdint[aij(i, j-1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factxp = (kdint[aij(i+1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxm = (kdint[aij(i-1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factyp = (kdint[aij(i, j+1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factym = (kdint[aij(i, j-1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy);
             diag[j] = 1.0 + factyp + factym;
             sup[j] = -factyp;
             inf[j] = -factym;
@@ -274,7 +292,19 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
             sup[0] = 0.0;
             f[0] = zint[aij(i, 0)];
         } else {
-            // FIXME: implement reflective boundary
+            // reflective boundary
+            // TODO: Check whether this j was actually supposed to be a 0 in the original
+            // (probably).
+            // factxp = (kdint[aij(i+1, 0)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxp = (kdint[aij(i+1, 0)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxm = (kdint[aij(i-1, 0)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factyp = (kdint[aij(i, 1)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            diag[0] = 1.0 + factyp;
+            sup[0] = -factyp;
+            f[0] =
+                zint[aij(i, 0)] + factxp * zint[aij(i+1, 0)] -
+                (factxp + factxm) *
+                zint[aij(i, 0)] + factxm * zint[aij(i-1, 0)];
         }
 /*
 ! top bc
@@ -296,7 +326,16 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
             inf[ny - 1] = 0.0;
             f[ny - 1] = zint[aij(i, ny - 1)];
         } else {
-            // FIXME: implement reflective boundary
+            // reflective boundary
+            factxp = (kdint[aij(i+1, ny-1)] + kdint[aij(i, ny-1)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxm = (kdint[aij(i-1, ny-1)] + kdint[aij(i, ny-1)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factym = (kdint[aij(i, ny-2)] + kdint[aij(i, ny-1)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            diag[ny - 1] = 1.0 + factym;
+            inf[ny - 1] = -factym;
+            f[ny - 1] =
+                zint[aij(i, ny-1)] + factxp * zint[aij(i+1, ny-1)] -
+                (factxp + factxm) *
+                zint[aij(i, ny-1)] + factxm * zint[aij(i-1, ny-1)];
         }
 /*
   call tridag (inf,diag,sup,f,res,ny)
@@ -331,7 +370,7 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
         for i in 0..nx {
             ij = aij(i, j);
             // FIXME: Track total erosion and erosion rate.
-            h[aij(i, j)] = zintp[aij(i, j)] as f32;
+            h[ij] = zintp[ij] as f32;
         }
     }
 
@@ -399,6 +438,13 @@ do 11 j=2,n
         gam[j] = c[j - 1] / bet;
         bet = b[j] - a[j] * gam[j];
         assert!(bet != 0.0);
+        // Round 0: u[0] = r[0] / b[0]
+        //               = r'[0] / b'[0]
+        // Round j: u[j] = (r[j] - a[j] * u'[j - 1]) / b'[j]
+        //               = (r[j] - a[j] * r'[j - 1] / b'[j - 1]) / b'[j]
+        //               = (r[j] - (a[j] / b'[j - 1]) * r'[j - 1]) / b'[j]
+        //               = (r[j] - w[j] * r'[j - 1]) / b'[j]
+        //               = r'[j] / b'[j]
         u[j] = (r[j] - a[j] * u[j - 1]) / bet;
     }
 /*
diff --git a/world/src/sim/erosion.rs b/world/src/sim/erosion.rs
index d9ae16cda9..db9997d45c 100644
--- a/world/src/sim/erosion.rs
+++ b/world/src/sim/erosion.rs
@@ -487,8 +487,8 @@ pub fn get_rivers(
 ///
 /// TODO: See if allocating in advance is worthwhile.
 fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync)) -> Box<[f64]> {
-    const MIN_MAX_ANGLE: f64 = 6.0/*30.0*//*6.0*//*15.0*/ / 360.0 * 2.0 * f64::consts::PI;
-    const MAX_MAX_ANGLE: f64 = 54.0/*50.0*//*54.0*//*45.0*/ / 360.0 * 2.0 * f64::consts::PI;
+    const MIN_MAX_ANGLE: f64 = 15.0/*6.0*//*30.0*//*6.0*//*15.0*/ / 360.0 * 2.0 * f64::consts::PI;
+    const MAX_MAX_ANGLE: f64 = 60.0/*54.0*//*50.0*//*54.0*//*45.0*/ / 360.0 * 2.0 * f64::consts::PI;
     const MAX_ANGLE_RANGE: f64 = MAX_MAX_ANGLE - MIN_MAX_ANGLE;
     let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
     h.par_iter()
@@ -497,7 +497,7 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
             let wposf = uniform_idx_as_vec2(posi).map(|e| e as f64) * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
             let wposz = z as f64 / height_scale;// * CONFIG.mountain_scale as f64;
             // Normalized to be between 6 and and 54 degrees.
-            let div_factor = /*32.0*//*16.0*//*64.0*/256.0/*1.0*//*4.0*/ * /*1.0*/4.0/* TerrainChunkSize::RECT_SIZE.x as f64 / 8.0 */;
+            let div_factor = /*32.0*//*16.0*//*64.0*//*256.0*/8.0/*8.0*//*1.0*//*4.0*/ * /*1.0*/16.0/* TerrainChunkSize::RECT_SIZE.x as f64 / 8.0 */;
             let rock_strength = rock_strength_nz
                 .get([
                     wposf.x, /* / div_factor*/
@@ -531,9 +531,9 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
             // TODO: Make all this stuff configurable... but honestly, it's so complicated that I'm not
             // sure anyone would be able to usefully tweak them on a per-map basis?  Plus it's just a
             // hacky heuristic anyway.
-            let center = 0.25;
+            let center = /*0.25*/0.4;
             let dmin = center - /*0.15;//0.05*/0.05;
-            let dmax = center + 0.05;//0.05;
+            let dmax = center + /*0.05*//*0.10*/0.05;//0.05;
             let log_odds = |x: f64| logit(x) - logit(center);
             let rock_strength = logistic_cdf(
                 1.0 * logit(rock_strength.min(1.0f64 - 1e-7).max(1e-7))
@@ -541,6 +541,7 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
             );
             // let rock_strength = 0.5;
             let max_slope = (rock_strength * MAX_ANGLE_RANGE + MIN_MAX_ANGLE).tan();
+            // let max_slope = /*30.0.to_radians().tan();*/3.0.sqrt() / 3.0;
             max_slope
         })
         .collect::<Vec<_>>()
@@ -554,7 +555,7 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
 ///   dh(p) / dt = uplift(p)−k * A(p)^m * slope(p)^n
 ///
 ///   where A(p) is the drainage area at p, m and n are constants
-///   (we choose m = 0.5 and n = 1), and k is a constant.  We choose
+///   (we choose m = 0.4 and n = 1), and k is a constant.  We choose
 ///
 ///   k = 2.244 * uplift.max() / (desired_max_height)
 ///
@@ -594,6 +595,10 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
 ///
 /// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90
 ///
+/// We also borrow the implementation for sediment transport from
+///
+/// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/StreamPowerLaw.f90
+///
 /// [1] Guillaume Cordonnier, Jean Braun, Marie-Paule Cani, Bedrich Benes, Eric Galin, et al..
 ///     Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion.
 ///     Computer Graphics Forum, Wiley, 2016, Proc. EUROGRAPHICS 2016, 35 (2), pp.165-175.
@@ -608,11 +613,14 @@ fn erode(
     done_val: bool,
     erosion_base: f32,
     max_uplift: f32,
+    max_g: f32,
     kdsed: f64,
     _seed: &RandomField,
     rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync),
-    uplift: impl Fn(usize) -> f32,
-    kd: impl Fn(usize) -> f64,
+    uplift: impl Fn(usize) -> f32 + Sync,
+    kf: impl Fn(usize) -> f32,
+    kd: impl Fn(usize) -> f32,
+    g: impl Fn(usize) -> f32 + Sync,
     is_ocean: impl Fn(usize) -> bool + Sync,
 ) {
     log::debug!("Done draining...");
@@ -636,12 +644,14 @@ fn erode(
     //
     //   max_uplift / height_scale m / dt / (5.010e-4 m / y) = 1
     //   (max_uplift / height_scale / 5.010e-4) y = dt
-    let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / 5.010e-4;
+    // 5e-7
+    let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / /*5.010e-4*/1e-3;
     println!("dt={:?}", dt);
     // Landslide constant: ideally scaled to 10e-2 m / y^-1
     let l = /*200.0 * max_uplift as f64;*/1.0e-2 /*/ CONFIG.mountain_scale as f64*/ * height_scale * dt;
     // Net precipitation rate (m / year)
     let p = 1.0 * height_scale;
+    let m = 0.4;
     // Stream power erosion constant (bedrock), in m^(1-2m) / year  (times dt).
     let k_fb = // erosion_base as f64 + 2.244 / mmaxh as f64 * max_uplift as f64;
         // 2.244*(5.010e-4)/512*5- (1.097e-5)
@@ -653,28 +663,40 @@ fn erode(
         // 2e-5 * dt;
         // 2.244/2048*5*32/(250000/4)*10^6
         // ln(tan(30/360*2*pi))-ln(tan(6/360*2*pi))*1500 = 3378
-        erosion_base as f64 + 2.244 / mmaxh as f64 * /*10.0*//*5.0*//*9.0*//*7.5*//*5.0*//*2.5*/1.5/*3.75*/ * max_uplift as f64;
-        // 1e-5 * dt;
+        //erosion_base as f64 + 2.244 / mmaxh as f64 * /*10.0*//*5.0*//*9.0*//*7.5*//*5.0*//*2.5*//*1.5*//*5.0*//*1.0*//*1.5*//*2.5*//*3.75*/ * max_uplift as f64;
+        // 2.5e-6 * dt;
+        2e-5 * dt;
+        // see http://geosci.uchicago.edu/~kite/doc/Whipple_and_Tucker_1999.pdf
+        //5e-6 * dt; // 2e-5 was designed for steady state uplift of 2mm / y whih would amount to 500 m / 250,000 y.
         // (2.244*(5.010e-4)/512)/(2.244*(5.010e-4)/2500) = 4.88...
         // 2.444 * 5
     // Stream power erosion constant (sediment), in m^(1-2m) / year (times dt).
-    let k_fs = k_fb * /*1.0*//*2.0*/2.0/*4.0*/;
+    let k_fs_mult = 2.0;/*1.5*/;
+    // let k_fs = k_fb * 1.0/*1.5*//*2.0*//*2.0*//*4.0*/;
     // u = k * h_max / 2.244
     // let uplift_scale = erosion_base as f64 + (k_fb * mmaxh / 2.244 / 5.010e-4 as f64 * mmaxh as f64) * dt;
-    let ((dh, indirection, newh, area), mut max_slopes) = rayon::join(
+    let ((dh, indirection, newh, maxh, area), (mut max_slopes, ht)) = rayon::join(
         || {
-            let mut dh = downhill(h, |posi| is_ocean(posi));
+            let mut dh = downhill(h, |posi| is_ocean(posi) && h[posi] <= 0.0);
             log::debug!("Computed downhill...");
-            let (boundary_len, indirection, newh) = get_lakes(&h, &mut dh);
+            let (boundary_len, indirection, newh, maxh) = get_lakes(&h, &mut dh);
             log::debug!("Got lakes...");
             let area = get_drainage(&newh, &dh, boundary_len);
             log::debug!("Got flux...");
-            (dh, indirection, newh, area)
+            (dh, indirection, newh, maxh, area)
         },
         || {
-            let max_slope = get_max_slope(h, rock_strength_nz);
-            log::debug!("Got max slopes...");
-            max_slope
+            rayon::join(
+                || {
+                    let max_slope = get_max_slope(h, rock_strength_nz);
+                    log::debug!("Got max slopes...");
+                    max_slope
+                },
+                || {
+                    // Store the elevation at t
+                    h.to_vec().into_boxed_slice()
+                },
+            )
         },
     );
 
@@ -683,140 +705,297 @@ fn erode(
     let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
     let chunk_area = neighbor_coef.x * neighbor_coef.y;
     // TODO: Make more principled.
-    let mid_slope = (30.0 / 360.0 * 2.0 * f64::consts::PI).tan();
-    // Iterate in ascending height order.
+    let mid_slope = (30.0 / 360.0 * 2.0 * f64::consts::PI).tan();//1.0;
+
+    let mut lake_water_volume = vec![/*-1i32*/0.0f64; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+    let mut elev = vec![/*-1i32*/0.0f64; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+    let mut hp = vec![/*-1i32*/0.0f64; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+    let mut deltah = vec![/*-1i32*/0.0f64; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+
+    // calculate the elevation / SPL, including sediment flux
+    let tol = 1.0e-4f64 * (maxh as f64 + 1.0);
+    let mut err = 2.0 * tol;
+
+    // Some variables for tracking statistics, currently only for debugging purposes.
     let mut maxh = 0.0;
     let mut nland = 0usize;
     let mut sums = 0.0;
     let mut sumh = 0.0;
-    for &posi in &*newh {
-        let posi = posi as usize;
-        let posj = dh[posi];
-        if posj < 0 {
-            if posj == -1 {
-                panic!("Disconnected lake!");
-            }
-            // wh for oceans is always 0.
-            wh[posi] = 0.0;
-            // max_slopes[posi] = kd(posi);
-        // Egress with no outgoing flows.
-        } else {
-            // *is_done.at(posi) = done_val;
-            let posj = posj as usize;
-            let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
+    let mut ntherm = 0usize;
 
-            // Has an outgoing flow edge (posi, posj).
-            // flux(i) = k * A[i]^m / ((p(i) - p(j)).magnitude()), and δt = 1
-            let neighbor_distance = (neighbor_coef * dxy).magnitude();
-            // Since the area is in meters^(2m) and neighbor_distance is in m, so long as m = 0.5,
-            // we have meters^(1) / meters^(1), so they should cancel out.  Otherwise, we would
-            // want to multiply neighbor_distance by height_scale and area[posi] by
-            // height_scale^2, to make sure we were working in the correct units for dz
-            // (which has height_scale height unit = 1.0 meters).
-            let old_h_i = h[posi] as f64;
-            let old_b_i = b[posi] as f64;
-            let uplift_i = uplift(posi) as f64;
-            let k = if (old_h_i - old_b_i as f64) > sediment_thickness {
-                // Sediment
-                k_fs
+    // Gauss-Seidel iteration
+
+    let mut lake_sediment = vec![/*-1i32*/0.0f64; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+    let mut lake_sill = vec![/*-1i32*/-1isize; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
+
+    let mut n_gs_stream_power_law = 0;
+    while err > tol && n_gs_stream_power_law < 99 {
+        log::info!("Stream power iteration #{:?}", n_gs_stream_power_law);
+
+        // Reset statistics in each loop.
+        maxh = 0.0;
+        nland = 0usize;
+        sums = 0.0;
+        sumh = 0.0;
+        ntherm = 0usize;
+
+        // Keep track of how many iterations we've gone to test for convergence.
+        n_gs_stream_power_law += 1;
+
+        rayon::join(
+            || {
+                // guess/update the elevation at t+Δt (k)
+                hp.par_iter_mut().zip(h.par_iter()).for_each(|(mut hp, h)| {
+                    *hp = *h as f64;
+                });
+            },
+            || {
+                // calculate erosion/deposition at each node
+                deltah.par_iter_mut().enumerate().for_each(|(posi, mut deltah)| {
+                    *deltah = (ht[posi] - h[posi]) as f64;
+                });
+            },
+        );
+
+        // sum the erosion in stack order
+        for &posi in newh.iter().rev() {
+            let posi = posi as usize;
+            let posj = dh[posi];
+            if posj < 0 {
+                lake_sediment[posi] = deltah[posi];
             } else {
-                // Bedrock
-                k_fb
-            };
-            // let k = k * uplift_i / max_uplift as f64;
-            let flux = k * (p * chunk_area * area[posi] as f64).sqrt() / neighbor_distance;
-            // h[i](t + dt) = (h[i](t) + δt * (uplift[i] + flux(i) * h[j](t + δt))) / (1 + flux(i) * δt).
-            // NOTE: posj has already been computed since it's downhill from us.
-            // Therefore, we can rely on wh being set to the water height for that node.
-            let h_j = h[posj] as f64;
-            let wh_j = wh[posj] as f64;
-            // hi(t + ∂t) = (hi(t) + ∂t(ui + kp^mAi^m(hj(t + ∂t)/||pi - pj||))) / (1 + ∂t * kp^mAi^m / ||pi - pj||)
-            let mut new_h_i = old_h_i + uplift_i;
-            // Only perform erosion if we are above the water level of the previous node.
-            if new_h_i > wh_j {
-                new_h_i = (new_h_i + flux * h_j) / (1.0 + flux);
-                // If we dipped below the receiver's water level, set our height to the receiver's
-                // water level.
-                if new_h_i <= wh_j {
-                    new_h_i = wh_j;
-                } else {
-                    let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
-                    let mag_slope = dz/*.abs()*/ / neighbor_distance;
-
-                    nland += 1;
-                    sumh += new_h_i;
-                    sums += mag_slope;
+                let posj = posj as usize;
+                deltah[posi] -= (ht[posi] as f64 - hp[posi]);
+                let lposi = lake_sill[posi];
+                if lposi == posi as isize {
+                    if deltah[posi] <= 0.0 {
+                        lake_sediment[posi] = 0.0;
+                    } else {
+                        lake_sediment[posi] = deltah[posi];
+                    }
                 }
-            }
-            // Set max_slope to this node's water height (max of receiver's water height and
-            // this node's height).
-            wh[posi] = wh_j.max(new_h_i) as f32;
-            // Prevent erosion from dropping us below our receiver, unless we were already below it.
-            // new_h_i = h_j.min(old_h_i + uplift_i).max(new_h_i);
-            // Find out if this is a lake bottom.
-            /* let indirection_idx = indirection[posi];
-            let is_lake_bottom = indirection_idx < 0;
-            let _fake_neighbor = is_lake_bottom && dxy.x.abs() > 1.0 && dxy.y.abs() > 1.0;
-            // Test the slope.
-            let max_slope = max_slopes[posi] as f64;
-            // Hacky version of thermal erosion: only consider lowest neighbor, don't redistribute
-            // uplift to other neighbors.
-            let (posk, h_k) = /* neighbors(posi)
-                .filter(|&posk| *is_done.at(posk) == done_val)
-                // .filter(|&posk| *is_done.at(posk) == done_val || is_ocean(posk))
-                .map(|posk| (posk, h[posk] as f64))
-                // .filter(|&(posk, h_k)| *is_done.at(posk) == done_val || h_k < 0.0)
-                .min_by(|&(_, a), &(_, b)| a.partial_cmp(&b).unwrap())
-                .unwrap_or((posj, h_j)); */
-                (posj, h_j);
-                // .max(h_j);
-            let (posk, h_k) = if h_k < h_j {
-                (posk, h_k)
-            } else {
-                (posj, h_j)
-            };
-            let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posk)).map(|e| e as f64);
-            let neighbor_distance = (neighbor_coef * dxy).magnitude();
-            let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
-            let mag_slope = dz/*.abs()*/ / neighbor_distance; */
-            // If you're on the lake bottom and not right next to your neighbor, don't compute a
-            // slope.
-            if
-            /* !is_lake_bottom */ /* !fake_neighbor */
-            true {
-                /* if
-                /* !is_lake_bottom && */
-                mag_slope > max_slope {
-                    // println!("old slope: {:?}, new slope: {:?}, dz: {:?}, h_j: {:?}, new_h_i: {:?}", mag_slope, max_slope, dz, h_j, new_h_i);
-                    // Thermal erosion says this can't happen, so we reduce dh_i to make the slope
-                    // exactly max_slope.
-                    // max_slope = (old_h_i + dh - h_j) / height_scale/* * CONFIG.mountain_scale */ / NEIGHBOR_DISTANCE
-                    // dh = max_slope * NEIGHBOR_DISTANCE * height_scale/* / CONFIG.mountain_scale */ + h_j - old_h_i.
-                    let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
-                    new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l * (mag_slope - max_slope));
-                    let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
-                    let slope = dz/*.abs()*/ / neighbor_distance;
-                    sums += slope;
-                    /* max_slopes[posi] = /*(mag_slope - max_slope) * */kd(posi);
-                    sums += mag_slope; */
-                // let slope = dz.signum() * max_slope;
-                // new_h_i = slope * neighbor_distance * height_scale /* / CONFIG.mountain_scale as f64 */ + h_j;
-                // sums += max_slope;
-                } else {
-                    // max_slopes[posi] = 0.0;
-                    sums += mag_slope;
-                    // Just use the computed rate.
-                } */
-                h[posi] = new_h_i as f32;
-                // Make sure to update the basement as well!
-                b[posi] = (old_b_i + uplift_i).min(new_h_i) as f32;
+                deltah[posi] += ht[posi] as f64 - hp[posi];
+                deltah[posj] += deltah[posi];
             }
         }
-        // *is_done.at(posi) = done_val;
-        maxh = h[posi].max(maxh);
+        // do ij=nn,1,-1
+        //   ijk=stack(ij)
+        //   ijr=rec(ijk)
+        //   if (ijr.ne.ijk) then
+        //     dh(ijk)=dh(ijk)-(ht(ijk)-hp(ijk))
+        //     if (lake_sill(ijk).eq.ijk) then
+        //       if (dh(ijk).le.0.d0) then
+        //         lake_sediment(ijk)=0.d0
+        //       else
+        //         lake_sediment(ijk)=dh(ijk)
+        //       endif
+        //     endif
+        //     dh(ijk)=dh(ijk)+(ht(ijk)-hp(ijk))
+        //     dh(ijr)=dh(ijr)+dh(ijk)
+        //   else
+        //     lake_sediment(ijk)=dh(ijk)
+        //   endif
+        // enddo
+
+        elev.par_iter_mut().enumerate().for_each(|(posi, mut elev)| {
+            if dh[posi] < 0 {
+                *elev = ht[posi] as f64;
+            } else {
+                *elev = ht[posi] as f64 + (deltah[posi] - (ht[posi] as f64 - hp[posi])) * g(posi) as f64 / area[posi] as f64;
+            }
+        });
+
+        // Iterate in ascending height order.
+        let mut sum_err = 0.0f64;
+        for &posi in &*newh {
+            let posi = posi as usize;
+            let posj = dh[posi];
+            if posj < 0 {
+                if posj == -1 {
+                    panic!("Disconnected lake!");
+                }
+                // wh for oceans is always 0.
+                wh[posi] = 0.0;
+                lake_sill[posi] = posi as isize;
+                lake_water_volume[posi] = 0.0;
+                // max_slopes[posi] = kd(posi);
+                // Egress with no outgoing flows.
+            } else {
+                // *is_done.at(posi) = done_val;
+                let posj = posj as usize;
+                let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
+
+                // Has an outgoing flow edge (posi, posj).
+                // flux(i) = k * A[i]^m / ((p(i) - p(j)).magnitude()), and δt = 1
+                let neighbor_distance = (neighbor_coef * dxy).magnitude();
+                // Since the area is in meters^(2m) and neighbor_distance is in m, so long as m = 0.5,
+                // we have meters^(1) / meters^(1), so they should cancel out.  Otherwise, we would
+                // want to multiply neighbor_distance by height_scale and area[posi] by
+                // height_scale^2, to make sure we were working in the correct units for dz
+                // (which has height_scale height unit = 1.0 meters).
+                let old_h_i = /*h*/elev[posi] as f64;
+                let old_b_i = b[posi] as f64;
+                let uplift_i = uplift(posi) as f64;
+                assert!(uplift_i.is_normal() && uplift_i == 0.0 || uplift_i.is_positive());
+                // h[i](t + dt) = (h[i](t) + δt * (uplift[i] + flux(i) * h[j](t + δt))) / (1 + flux(i) * δt).
+                // NOTE: posj has already been computed since it's downhill from us.
+                // Therefore, we can rely on wh being set to the water height for that node.
+                let h_j = h[posj] as f64;
+                let wh_j = wh[posj] as f64;
+                let mut new_h_i = old_h_i + uplift_i;
+                // Only perform erosion if we are above the water level of the previous node.
+                if new_h_i > wh_j {
+                    // hi(t + ∂t) = (hi(t) + ∂t(ui + kp^mAi^m(hj(t + ∂t)/||pi - pj||))) / (1 + ∂t * kp^mAi^m / ||pi - pj||)
+                    let k = if (old_h_i - old_b_i as f64) > sediment_thickness {
+                        // Sediment
+                        // k_fs
+                        k_fs_mult * kf(posi) as f64
+                    } else {
+                        // Bedrock
+                        // k_fb
+                        kf(posi) as f64
+                    } * dt;
+                    // let k = k * uplift_i / max_uplift as f64;
+                    let flux = k * (p * chunk_area * area[posi] as f64).powf(m) / neighbor_distance;
+                    assert!(flux.is_normal() && flux.is_positive());
+                    new_h_i = (new_h_i + flux * h_j) / (1.0 + flux);
+                    lake_sill[posi] = posi as isize;
+                    lake_water_volume[posi] = 0.0;
+                    // If we dipped below the receiver's water level, set our height to the receiver's
+                    // water level.
+                    if new_h_i <= wh_j {
+                        new_h_i = wh_j;
+                    } else {
+                        let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
+                        let mag_slope = dz/*.abs()*/ / neighbor_distance;
+
+                        nland += 1;
+                        sumh += new_h_i;
+                        sums += mag_slope;
+                    }
+                } else {
+                    let lposj = lake_sill[posj];
+                    lake_sill[posi] = lposj;
+                    if lposj >= 0 {
+                        let lposj = lposj as usize;
+                        lake_water_volume[lposj] += wh_j - new_h_i;
+                    }
+                }
+                // Set max_slope to this node's water height (max of receiver's water height and
+                // this node's height).
+                wh[posi] = wh_j.max(new_h_i) as f32;
+                // Prevent erosion from dropping us below our receiver, unless we were already below it.
+                // new_h_i = h_j.min(old_h_i + uplift_i).max(new_h_i);
+                // Find out if this is a lake bottom.
+                /* let indirection_idx = indirection[posi];
+                let is_lake_bottom = indirection_idx < 0;
+                let _fake_neighbor = is_lake_bottom && dxy.x.abs() > 1.0 && dxy.y.abs() > 1.0;
+                // Test the slope.
+                let max_slope = max_slopes[posi] as f64;
+                // Hacky version of thermal erosion: only consider lowest neighbor, don't redistribute
+                // uplift to other neighbors.
+                let (posk, h_k) = /* neighbors(posi)
+                    .filter(|&posk| *is_done.at(posk) == done_val)
+                    // .filter(|&posk| *is_done.at(posk) == done_val || is_ocean(posk))
+                    .map(|posk| (posk, h[posk] as f64))
+                    // .filter(|&(posk, h_k)| *is_done.at(posk) == done_val || h_k < 0.0)
+                    .min_by(|&(_, a), &(_, b)| a.partial_cmp(&b).unwrap())
+                    .unwrap_or((posj, h_j)); */
+                    (posj, h_j);
+                    // .max(h_j);
+                let (posk, h_k) = if h_k < h_j {
+                    (posk, h_k)
+                } else {
+                    (posj, h_j)
+                };
+                let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posk)).map(|e| e as f64);
+                let neighbor_distance = (neighbor_coef * dxy).magnitude();
+                let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
+                let mag_slope = dz/*.abs()*/ / neighbor_distance; */
+                // If you're on the lake bottom and not right next to your neighbor, don't compute a
+                // slope.
+                if
+                /* !is_lake_bottom */ /* !fake_neighbor */
+                true {
+                    /* if
+                    /* !is_lake_bottom && */
+                    mag_slope > max_slope {
+                        // println!("old slope: {:?}, new slope: {:?}, dz: {:?}, h_j: {:?}, new_h_i: {:?}", mag_slope, max_slope, dz, h_j, new_h_i);
+                        // Thermal erosion says this can't happen, so we reduce dh_i to make the slope
+                        // exactly max_slope.
+                        // max_slope = (old_h_i + dh - h_j) / height_scale/* * CONFIG.mountain_scale */ / NEIGHBOR_DISTANCE
+                        // dh = max_slope * NEIGHBOR_DISTANCE * height_scale/* / CONFIG.mountain_scale */ + h_j - old_h_i.
+                        let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
+                        new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l * (mag_slope - max_slope));
+                        let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
+                        let slope = dz/*.abs()*/ / neighbor_distance;
+                        sums += slope;
+                        /* max_slopes[posi] = /*(mag_slope - max_slope) * */kd(posi);
+                        sums += mag_slope; */
+                    // let slope = dz.signum() * max_slope;
+                    // new_h_i = slope * neighbor_distance * height_scale /* / CONFIG.mountain_scale as f64 */ + h_j;
+                    // sums += max_slope;
+                    } else {
+                        // max_slopes[posi] = 0.0;
+                        sums += mag_slope;
+                        // Just use the computed rate.
+                    } */
+                    h[posi] = new_h_i as f32;
+                    // Make sure to update the basement as well!
+                    // b[posi] = (old_b_i + uplift_i).min(new_h_i) as f32;
+                }
+            }
+            // *is_done.at(posi) = done_val;
+            maxh = h[posi].max(maxh);
+
+            // Add sum of squares of errors.
+            sum_err += (h[posi] as f64 - hp[posi]).powi(2);
+        }
+
+        err = (sum_err / newh.len() as f64).sqrt();
+        if max_g == 0.0 {
+            err = 0.0;
+        }
+        if n_gs_stream_power_law == 99 {
+            log::warn!("Beware: Gauss-Siedel scheme not convergent");
+        }
     }
+
+    //b=min(h,b)
+
+    // update the basement
+    b.par_iter_mut().zip(h.par_iter()).enumerate().for_each(|(posi, (mut b, h))| {
+        let old_b_i = *b;
+        let uplift_i = uplift(posi);
+
+        *b = (old_b_i + uplift_i).min(*h);
+    });
+
+    // update the height to reflect sediment flux.
+    h.par_iter_mut().enumerate().for_each(|(posi, mut h)| {
+        let lposi = lake_sill[posi];
+        if lposi >= 0 {
+            let lposi = lposi as usize;
+            if lake_water_volume[lposi] > 0.0 {
+                // +max(0.d0,min(lake_sediment(lake_sill(ij)),lake_water_volume(lake_sill(ij))))/
+                // lake_water_volume(lake_sill(ij))*(water(ij)-h(ij))
+                *h +=
+                    (0.0.max(lake_sediment[lposi].min(lake_water_volume[lposi])) /
+                    lake_water_volume[lposi] *
+                    (wh[posi] - *h) as f64) as f32;
+            }
+        }
+    });
+    // do ij=1,nn
+    //   if (lake_sill(ij).ne.0) then
+    //     if (lake_water_volume(lake_sill(ij)).gt.0.d0) h(ij)=h(ij) &
+    //     +max(0.d0,min(lake_sediment(lake_sill(ij)),lake_water_volume(lake_sill(ij))))/ &
+    //     lake_water_volume(lake_sill(ij))*(water(ij)-h(ij))
+    //   endif
+    // enddo
+
     log::info!(
-        "Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?})",
+        "Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?}) (num land: {:?}) (num thermal: {:?})",
         maxh,
         if nland == 0 {
             f64::INFINITY
@@ -828,30 +1007,38 @@ fn erode(
         } else {
             sums / nland as f64
         },
+        nland,
+        ntherm,
     );
 
     // Apply thermal erosion.
     maxh = 0.0;
     sumh = 0.0;
     sums = 0.0;
+    nland = 0usize;
+    ntherm = 0usize;
     for &posi in &*newh {
         let posi = posi as usize;
         let posj = dh[posi];
-        let old_h_i = h[posi] as f64;
+        let old_h_i = /*h*/b[posi] as f64;
         let old_b_i = b[posi] as f64;
         let max_slope = max_slopes[posi] as f64;
         // Remember k_d for this chunk in max_slopes.
-        max_slopes[posi] = if (old_h_i - old_b_i as f64) > sediment_thickness {
+        let kd_factor =
+            // 1.0;
+            (1.0 / (max_slope / mid_slope/*.sqrt()*//*.powf(0.03125)*/).powf(/*2.0*/1.0))/*.min(kdsed)*/;
+        max_slopes[posi] = if (old_h_i - old_b_i as f64) > sediment_thickness && kdsed > 0.0 {
             // Sediment
-            kdsed
+            kdsed/* * kd_factor*/
         } else {
             // Bedrock
-            kd(posi) / (max_slope / mid_slope/*.sqrt()*//*.powf(0.03125)*/)/*.sqrt()*/
+            kd(posi) as f64 / kd_factor
         };
         if posj < 0 {
             if posj == -1 {
                 panic!("Disconnected lake!");
             }
+            wh[posi] = 0.0;
             // Egress with no outgoing flows.
         } else {
             let posj = posj as usize;
@@ -860,13 +1047,15 @@ fn erode(
             let wh_j = wh[posj] as f64;
             // If you're on the lake bottom and not right next to your neighbor, don't compute a
             // slope.
+            let mut new_h_i = old_h_i;//old_b_i;
             if
             /* !is_lake_bottom */ /* !fake_neighbor */
             wh_j < old_h_i {
-                let mut new_h_i = old_h_i;//old_b_i;
-                // NOTE: Currently assuming that talus angle is the same once the substance is
-                // submerged in water, but actually that's probably not true.
-                let h_j = h[posj] as f64/*wh_j*/;
+                // NOTE: Currently assuming that talus angle is not eroded once the substance is
+                // totally submerged in water, and that talus angle if part of the substance is
+                // in water is 0 (or the same as the dry part, if this is set to wh_j), but
+                // actually that's probably not true.
+                let h_j = /*h[posj] as f64*/wh_j;
                 // let h_j = b[posj] as f64;
                 /* let indirection_idx = indirection[posi];
                 let is_lake_bottom = indirection_idx < 0;
@@ -902,30 +1091,56 @@ fn erode(
                     // dh = max_slope * NEIGHBOR_DISTANCE * height_scale/* / CONFIG.mountain_scale */ + h_j - old_h_i.
                     let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
                     new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l * (mag_slope - max_slope));
-                    let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
-                    let slope = dz/*.abs()*/ / neighbor_distance;
-                    sums += slope;
-                    /* max_slopes[posi] = /*(mag_slope - max_slope) * */kd(posi);
-                    sums += mag_slope; */
-                // let slope = dz.signum() * max_slope;
-                // new_h_i = slope * neighbor_distance * height_scale /* / CONFIG.mountain_scale as f64 */ + h_j;
-                // sums += max_slope;
+                    if new_h_i <= wh_j {
+                        new_h_i = wh_j;
+                    } else {
+                        let dz = (new_h_i - /*h_j*//*h_k*/wh_j).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
+                        let slope = dz/*.abs()*/ / neighbor_distance;
+                        sums += slope;
+                        // max_slopes[posi] = /*(mag_slope - max_slope) * */max_slopes[posi].max(kdsed);
+                        /* max_slopes[posi] = /*(mag_slope - max_slope) * */kd(posi);
+                        sums += mag_slope; */
+                        /* if kd_factor < 1.0 {
+                            max_slopes[posi] /= kd_factor;
+                        } else {
+                            max_slopes[posi] *= kd_factor;
+                        } */
+                        // max_slopes[posi] *= kd_factor;
+                        nland += 1;
+                        sumh += new_h_i;
+                        // let slope = dz.signum() * max_slope;
+                        // new_h_i = slope * neighbor_distance * height_scale /* / CONFIG.mountain_scale as f64 */ + h_j;
+                        // sums += max_slope;
+                    }
+                    ntherm += 1;
                 } else {
+                    /*if kd_factor < 1.0 {
+                        max_slopes[posi] *= kd_factor;
+                    }*/
+                    /* if (old_h_i - old_b_i as f64) <= sediment_thickness {
+                        max_slopes[posi] *= kd_factor;
+                    } */
+                    // max_slopes[posi] *= kd_factor;
                     sums += mag_slope;
                     // Just use the computed rate.
+                    nland += 1;
+                    sumh += new_h_i;
                 }
-                h/*b*/[posi] = new_h_i as f32;
+                // h/*b*/[posi] = old_h_i.min(new_h_i) as f32;
                 // Make sure to update the basement as well!
                 // b[posi] = old_b_i.min(new_h_i) as f32;
                 b[posi] = old_b_i.min(old_b_i + (new_h_i - old_h_i)) as f32;
-                sumh += new_h_i;
+                // sumh += new_h_i;
             }
+            // Set wh to this node's water height (max of receiver's water height and
+            // this node's height).
+            wh[posi] = wh_j.max(new_h_i) as f32;
         }
-        *is_done.at(posi) = done_val;
+        // *is_done.at(posi) = done_val;
         maxh = h[posi].max(maxh);
     }
     log::debug!(
-        "Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (avg slope: {:?})",
+        "Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (avg slope: {:?}) (num land: {:?}) (num thermal: {:?})",
         maxh,
         if nland == 0 {
             f64::INFINITY
@@ -937,6 +1152,8 @@ fn erode(
         } else {
             sums / nland as f64
         },
+        nland,
+        ntherm,
     );
 
     // Apply hillslope diffusion.
@@ -949,20 +1166,8 @@ fn erode(
               |posi| max_slopes[posi]/*kd*/,
               /* kdsed */-1.0,
     );
-    log::debug!(
-        "Done eroding (max height: {:?}) (avg height: {:?}) (avg slope: {:?})",
-        maxh,
-        if nland == 0 {
-            f64::INFINITY
-        } else {
-            sumh / nland as f64
-        },
-        if nland == 0 {
-            f64::INFINITY
-        } else {
-            sums / nland as f64
-        },
-    );
+    log::debug!("Done applying diffusion.");
+    log::debug!("Done eroding.");
 }
 
 /// The Planchon-Darboux algorithm for extracting drainage networks.
@@ -990,6 +1195,7 @@ pub fn fill_sinks(
         .map(|(posi, &h)| {
             let is_near_edge = is_ocean(posi);
             if is_near_edge {
+                debug_assert!(h <= 0.0);
                 h
             } else {
                 infinity
@@ -1051,7 +1257,7 @@ pub fn fill_sinks(
 /// - A list of chunks on the boundary (non-lake egress points).
 /// - The second indirection vector (associating chunk indices with their lake's adjacency list).
 /// - The adjacency list (stored in a single vector), indexed by the second indirection vector.
-pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[u32]>) {
+pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[u32]>, f32) {
     // Associates each lake index with its root node (the deepest one in the lake), and a list of
     // adjacent lakes.  The list of adjacent lakes includes the lake index of the adjacent lake,
     // and a node index in the adjacent lake which has a neighbor in this lake.  The particular
@@ -1124,6 +1330,7 @@ pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[
     log::debug!("Old lake roots: {:?}", lake_roots.len());
 
     let newh = newh.into_boxed_slice();
+    let mut maxh = -f32::INFINITY;
     // Now, we know that the sum of all the indirection nodes will be the same as the number of
     // nodes.  We can allocate a *single* vector with 8 * nodes entries, to be used as storage
     // space, and augment our indirection vector with the starting index, resulting in a vector of
@@ -1134,6 +1341,8 @@ pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[
         let lake_idx_ = indirection_[chunk_idx];
         let lake_idx = lake_idx_ as usize;
         let height = h[chunk_idx_ as usize];
+        // While we're here, compute the max elevation difference from zero among all nodes.
+        maxh = maxh.max(height.abs());
         // For every neighbor, check to see whether it is already set; if the neighbor is set,
         // its height is ≤ our height.  We should search through the edge list for the
         // neighbor's lake to see if there's an entry; if not, we insert, and otherwise we
@@ -1429,7 +1638,7 @@ pub fn get_lakes(h: &[f32], downhill: &mut [isize]) -> (usize, Box<[i32]>, Box<[
             }
         });
     assert_eq!(newh.len(), downhill.len());
-    (boundary.len(), indirection, newh.into_boxed_slice())
+    (boundary.len(), indirection, newh.into_boxed_slice(), maxh)
 }
 
 /// Perform erosion n times.
@@ -1443,6 +1652,9 @@ pub fn do_erosion(
     oldb: impl Fn(usize) -> f32 + Sync,
     is_ocean: impl Fn(usize) -> bool + Sync,
     uplift: impl Fn(usize) -> f32 + Sync,
+    kf: impl Fn(usize) -> f32 + Sync,
+    kd: impl Fn(usize) -> f32 + Sync,
+    g: impl Fn(usize) -> f32 + Sync,
 ) -> (Box<[f32]>, Box<[f32]>) {
     let oldh_ = (0..WORLD_SIZE.x * WORLD_SIZE.y)
         .into_par_iter()
@@ -1455,6 +1667,24 @@ pub fn do_erosion(
         .map(|posi| oldb(posi))
         .collect::<Vec<_>>()
         .into_boxed_slice();
+    // Stream power law erodability constant for fluvial erosion (bedrock)
+    let kf = (0..WORLD_SIZE.x * WORLD_SIZE.y)
+        .into_par_iter()
+        .map(|posi| kf(posi))
+        .collect::<Vec<_>>()
+        .into_boxed_slice();
+    // Stream power law erodability constant for hillslope diffusion (bedrock)
+    let kd = (0..WORLD_SIZE.x * WORLD_SIZE.y)
+        .into_par_iter()
+        .map(|posi| kd(posi))
+        .collect::<Vec<_>>()
+        .into_boxed_slice();
+    // Deposition coefficient
+    let g = (0..WORLD_SIZE.x * WORLD_SIZE.y)
+        .into_par_iter()
+        .map(|posi| g(posi))
+        .collect::<Vec<_>>()
+        .into_boxed_slice();
     let mut wh = vec![0.0; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
     // TODO: Don't do this, maybe?
     // (To elaborate, maybe we should have varying uplift or compute it some other way).
@@ -1475,6 +1705,11 @@ pub fn do_erosion(
         .cloned()
         .max_by(|a, b| a.partial_cmp(&b).unwrap())
         .unwrap();
+    let max_g = g
+        .into_par_iter()
+        .cloned()
+        .max_by(|a, b| a.partial_cmp(&b).unwrap())
+        .unwrap();
     log::debug!("Max uplift: {:?}", max_uplift);
     // Height of terrain, including sediment.
     let mut h = oldh_;
@@ -1488,12 +1723,15 @@ pub fn do_erosion(
     let k_fb = /*(erosion_base as f64 + 2.244 / mmaxh as f64 * /*10.0*//*5.0*//*9.0*//*7.5*//*5.0*//*2.5*//*1.5*/4.0/*1.0*//*3.75*/ * max_uplift as f64) / dt;*/
         2.0e-5 * dt;
     let kd_bedrock =
-        /*1e-2*//*0.25e-2*/1e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
+        /*1e-2*//*0.25e-2*/1e-2 / 1.0 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
         /* * k_fb / 2e-5 */;
     let kdsed =
-        /*1.5e-2*//*1e-4*//*1.25e-2*/1.5e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
+        /*1.5e-2*//*1e-4*//*1.25e-2*/1.5e-2 / 1.0 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
         /* * k_fb / 2e-5 */;
-    let kd = |posi: usize| kd_bedrock; // if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
+    // let kd = |posi: usize| kd_bedrock; // if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
+    let kd = |posi: usize| kd[posi]; // if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
+    let kf = |posi: usize| kf[posi];
+    let g = |posi: usize| g[posi];
     // Hillslope diffusion coefficient for sediment.
     let mut is_done = bitbox![0; WORLD_SIZE.x * WORLD_SIZE.y];
     for i in 0..n {
@@ -1509,11 +1747,15 @@ pub fn do_erosion(
             i & 1 == 0,
             erosion_base,
             max_uplift,
-            kdsed,
+            max_g,
+            -1.0,
+            // kdsed,
             seed,
             rock_strength_nz,
             |posi| uplift[posi],
+            kf,
             kd,
+            g,
             |posi| is_ocean(posi),
         );
     }
diff --git a/world/src/sim/mod.rs b/world/src/sim/mod.rs
index 95462eb56c..fe85b86312 100644
--- a/world/src/sim/mod.rs
+++ b/world/src/sim/mod.rs
@@ -51,8 +51,8 @@ use vek::*;
 // but I think that is probably less important since I don't think we actually cast a chunk id to
 // float, just coordinates... could be wrong though!
 pub const WORLD_SIZE: Vec2<usize> = Vec2 {
-    x: 1024 * 2,
-    y: 1024 * 2,
+    x: 1024,
+    y: 1024,
 };
 
 /// A structure that holds cached noise values and cumulative distribution functions for the input
@@ -107,6 +107,19 @@ pub(crate) struct GenCtx {
     pub uplift_nz: Worley,
 }
 
+pub struct WorldOpts {
+    /// Set to false to disable seeding elements during worldgen.
+    pub seed_elements: bool,
+}
+
+impl Default for WorldOpts {
+    fn default() -> Self {
+        Self {
+            seed_elements: true,
+        }
+    }
+}
+
 pub struct WorldSim {
     pub seed: u32,
     pub(crate) chunks: Vec<SimChunk>,
@@ -117,10 +130,10 @@ pub struct WorldSim {
 }
 
 impl WorldSim {
-    pub fn generate(seed: u32) -> Self {
+    pub fn generate(seed: u32, opts: WorldOpts) -> Self {
         let mut rng = ChaChaRng::from_seed(seed_expan::rng_state(seed));
         let continent_scale = 5_000.0f64/*32768.0*/.div(32.0).mul(TerrainChunkSize::RECT_SIZE.x as f64);
-        let rock_lacunarity = 0.5/*HybridMulti::DEFAULT_LACUNARITY*/;
+        let rock_lacunarity = 0.5/*2.0*//*HybridMulti::DEFAULT_LACUNARITY*/;
 
         let gen_ctx = GenCtx {
             turb_x_nz: SuperSimplex::new().set_seed(rng.gen()),
@@ -128,6 +141,7 @@ impl WorldSim {
             chaos_nz: RidgedMulti::new()
                 .set_octaves(/*7*//*3*/ /*7*//*3*/7)
                 .set_frequency(
+                    // RidgedMulti::DEFAULT_FREQUENCY * (5_000.0 / continent_scale)
                     /*RidgedMulti::DEFAULT_FREQUENCY **/ 3_000.0 * 8.0 / continent_scale,
                 )
                 // .set_persistence(RidgedMulti::DEFAULT_LACUNARITY.powf(-(1.0 - 0.5)))
@@ -144,10 +158,11 @@ impl WorldSim {
                 // .set_frequency(1.0 / ((1 << 0) as f64))
                 // .set_lacunarity(1.0)
                 // persistence = lacunarity^(-(1.0 - fractal increment))
-                .set_persistence(HybridMulti_::DEFAULT_LACUNARITY.powf(-(1.0 - /*0.75*/0.75)))
+                .set_lacunarity(HybridMulti_::DEFAULT_LACUNARITY)
+                .set_persistence(HybridMulti_::DEFAULT_LACUNARITY.powf(-(1.0 - /*0.75*/0.0)))
                 // .set_persistence(/*0.5*//*0.5*/0.5 + 1.0 / ((1 << 6) as f64))
                 // .set_offset(/*0.7*//*0.5*//*0.75*/0.7)
-                // .set_offset(/*0.7*//*0.5*//*0.75*/0.0)
+                .set_offset(/*0.7*//*0.5*//*0.75*/0.0)
                 .set_seed(rng.gen()),
             //temp_nz: SuperSimplex::new().set_seed(rng.gen()),
             temp_nz: Fbm::new()
@@ -192,9 +207,11 @@ impl WorldSim {
             town_gen: StructureGen2d::new(rng.gen(), 2048, 1024),
             river_seed: RandomField::new(rng.gen()),
             rock_strength_nz: /*Fbm*/HybridMulti_/*BasicMulti*//*Fbm*/::new()
-            .set_octaves(/*6*//*5*//*4*//*5*//*4*/6)
+            .set_octaves(/*6*//*5*//*4*//*5*//*4*//*6*/10)
+            .set_lacunarity(rock_lacunarity)
             // persistence = lacunarity^(-(1.0 - fractal increment))
             // NOTE: In paper, fractal increment is roughly 0.25.
+            .set_offset(0.0)
             .set_persistence(/*0.9*/ /*2.0*//*1.5*//*HybridMulti::DEFAULT_LACUNARITY*/rock_lacunarity.powf(-(1.0 - 0.25/*0.9*/)))
             // 256*32/2^4
             // (0.5^(-(1.0-0.9)))^4/256/32*2^4*16*32
@@ -202,7 +219,7 @@ impl WorldSim {
             // (0.5^(-(1.0-0.9)))^1/256/32*2^4*256*4
             // (2^(-(1.0-0.9)))^4
             // 16.0
-            .set_frequency(/*0.9*/ /*Fbm*/HybridMulti_::DEFAULT_FREQUENCY / (64.0/*256.0*//*1.0*//*16.0*/ * 32.0/* TerrainChunkSize::RECT_SIZE.x as f64 */))
+            .set_frequency(/*0.9*/ /*Fbm*//*HybridMulti_::DEFAULT_FREQUENCY*/1.0 / (8.0/*8.0*//*256.0*//*1.0*//*16.0*/ * 32.0/* TerrainChunkSize::RECT_SIZE.x as f64 */))
             // .set_persistence(/*0.9*/ /*2.0*/0.67)
             // .set_frequency(/*0.9*/ Fbm::DEFAULT_FREQUENCY / (2.0 * 32.0))
             // .set_lacunarity(0.5)
@@ -210,9 +227,10 @@ impl WorldSim {
             uplift_nz: Worley::new()
                 .set_seed(rng.gen())
                 .set_frequency(1.0 / (TerrainChunkSize::RECT_SIZE.x as f64 * 256.0))
-                .set_displacement(0.5)
-                .set_range_function(RangeFunction::Chebyshev)
-                .enable_range(true),
+                // .set_displacement(/*0.5*/0.0)
+                .set_displacement(/*0.5*/1.0)
+                .set_range_function(RangeFunction::Euclidean)
+                // .enable_range(true),
         };
 
         let river_seed = &gen_ctx.river_seed;
@@ -229,12 +247,30 @@ impl WorldSim {
             .set_scale(1.0 / rock_strength_scale); */
 
         let height_scale = 1.0f64; // 1.0 / CONFIG.mountain_scale as f64;
-        let max_erosion_per_delta_t = /*32.0*//*128.0*/128.0 * height_scale;
+        let max_erosion_per_delta_t = 128.0/*128.0*//*32.0*/ * height_scale;
         let erosion_pow_low = /*0.25*//*1.5*//*2.0*//*0.5*//*4.0*//*0.25*//*1.0*//*2.0*//*1.5*//*1.5*//*0.35*//*0.43*//*0.5*//*0.45*//*0.37*/1.002;
         let erosion_pow_high = /*1.5*//*1.0*//*0.55*//*0.51*//*2.0*/1.002;
         let erosion_center = /*0.45*//*0.75*//*0.75*//*0.5*//*0.75*/0.5;
-        let n_steps = /*100*//*50*//*100*/100/*100*//*37*/;//150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200;
-        let n_small_steps = 8;//8;//8; // 8
+        let n_steps = /*200*/50; // /*100*//*50*//*100*//*100*//*50*//*25*/25/*100*//*37*/;//150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200;
+        let n_small_steps = 8;//50;//8;//8;//8;//8;//8; // 8
+
+
+        // Logistic regression.  Make sure x ∈ (0, 1).
+        let logit = |x: f64| x.ln() - (-x).ln_1p();
+        // 0.5 + 0.5 * tanh(ln(1 / (1 - 0.1) - 1) / (2 * (sqrt(3)/pi)))
+        let logistic_2_base = 3.0f64.sqrt() * f64::consts::FRAC_2_PI;
+        let logistic_base = /*3.0f64.sqrt() * f64::consts::FRAC_1_PI*/1.0f64;
+        // Assumes μ = 0, σ = 1
+        let logistic_cdf = |x: f64| (x / logistic_2_base).tanh() * 0.5 + 0.5;
+
+        let exp_inverse_cdf = |x: f64/*, pow: f64*/| -(-x).ln_1p()/* / ln(pow)*/;
+        // 2 / pi * ln(tan(pi/2 * p))
+        let hypsec_inverse_cdf = |x: f64| f64::consts::FRAC_2_PI * ((x * f64::consts::FRAC_PI_2).tan().ln());
+
+        let min_epsilon =
+            1.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64).max(f64::EPSILON as f64 * 0.5);
+        let max_epsilon = (1.0 - 1.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64))
+            .min(1.0 - f64::EPSILON as f64 * 0.5);
 
         // fractal dimension should be between 0 and 0.9999...
         // (lacunarity^octaves)^(-H) = persistence^(octaves)
@@ -293,9 +329,9 @@ impl WorldSim {
                             /* .mul(0.25)
                             .add(0.125) */)
                             // .add(0.5)
-                        // .sub(0.05)
+                        .sub(0.05)
                         // .add(0.05)
-                        .add(0.075)
+                        // .add(0.075)
                         .mul(0.35), /*-0.0175*/
                     )
                 })
@@ -458,15 +494,70 @@ impl WorldSim {
         let is_ocean = get_oceans(old_height);
         let is_ocean_fn = |posi: usize| is_ocean[posi];
 
+        let uplift_nz_dist = gen_ctx.uplift_nz
+            .clone()
+            .enable_range(true);
         // Recalculate altitudes without oceans.
         // NaNs in these uniform vectors wherever pure_water() returns true.
-        let (alt_old_no_ocean, alt_old_inverse) = uniform_noise(|posi, _| {
-            if is_ocean_fn(posi) {
-                None
-            } else {
-                Some(old_height(posi) /*.abs()*/)
-            }
-        });
+        let ((alt_old_no_ocean, alt_old_inverse), (uplift_uniform, _)) = rayon::join(
+            || {
+                uniform_noise(|posi, _| {
+                    if is_ocean_fn(posi) {
+                        None
+                    } else {
+                        Some(old_height(posi) /*.abs()*/)
+                    }
+                })
+            },
+            || {
+                uniform_noise(|posi, wposf| {
+                    if is_ocean_fn(posi) {
+                        None
+                    } else {
+                        let turb_wposf =
+                            wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
+                        let turb = Vec2::new(
+                            gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
+                            gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
+                        ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+                        // let turb = Vec2::zero();
+                        let turb_wposf = wposf + turb;
+                        let turb_wposi = turb_wposf
+                            .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64)
+                            .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0));
+                        let turb_posi = vec2_as_uniform_idx(turb_wposi);
+                        let udist = uplift_nz_dist.get(turb_wposf.into_array())
+                            .min(1.0)
+                            .max(-1.0)
+                            .mul(0.5)
+                            .add(0.5);
+                        let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array())
+                            /* .min(0.5)
+                            .max(-0.5)*/
+                            .min(1.0)
+                            .max(-1.0)
+                            .mul(0.5)
+                            .add(0.5);
+                        let oheight = alt_old[/*(turb_posi / 64) * 64*/posi].0 as f64 - 0.5;
+                        assert!(udist >= 0.0);
+                        assert!(udist <= 1.0);
+                        let uheight_1 = uheight;//.powf(2.0);
+                        let udist_1 = (0.5 - udist).mul(2.0).max(0.0);
+                        let udist_2 = udist.mul(2.0).min(1.0);
+                        let height =
+                            // uheight_1;
+                            // uheight_1 * (/*udist_2*/udist.powf(2.0) * (f64::consts::PI * 2.0 * (1.0 / (1.0 - udist).max(f64::EPSILON)).min(2.5)/*udist * 5.0*/ * 2.0).cos().mul(0.5).add(0.5));
+                            // uheight * udist_ * (udist_ * 4.0 * 2 * f64::consts::PI).sin()
+                            uheight /* * 0.8*/ * udist_1.powf(2.0) +
+                            /*exp_inverse_cdf*/(oheight.max(0.0).min(max_epsilon).abs())/*.powf(2.0)*/ * 0.2 * udist_2.powf(2.0);
+                            // * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/
+                        // let height = uheight * (0.5 - udist) * 0.8 + (oheight.signum() * oheight.max(0.0).abs().powf(2.0)) * 0.2;// * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/
+
+                        Some(height)
+                    }
+                })
+            },
+        );
 
         let old_height_uniform = |posi: usize| alt_old_no_ocean[posi].0;
         let alt_old_min_uniform = 0.0;
@@ -494,17 +585,6 @@ impl WorldSim {
 
         // Perform some erosion.
 
-        // Logistic regression.  Make sure x ∈ (0, 1).
-        let logit = |x: f64| x.ln() - (-x).ln_1p();
-        // 0.5 + 0.5 * tanh(ln(1 / (1 - 0.1) - 1) / (2 * (sqrt(3)/pi)))
-        let logistic_2_base = 3.0f64.sqrt() * f64::consts::FRAC_2_PI;
-        let logistic_base = /*3.0f64.sqrt() * f64::consts::FRAC_1_PI*/1.0f64;
-        // Assumes μ = 0, σ = 1
-        let logistic_cdf = |x: f64| (x / logistic_2_base).tanh() * 0.5 + 0.5;
-
-        let exp_inverse_cdf = |x: f64/*, pow: f64*/| -(-x).ln_1p()/* / ln(pow)*/;
-        // 2 / pi * ln(tan(pi/2 * p))
-        let hypsec_inverse_cdf = |x: f64| f64::consts::FRAC_2_PI * ((x * f64::consts::FRAC_PI_2).tan().ln());
         // 2^((2^10-2)/256) = 15.91...
         // -ln(1-(1-(2^(-22)*0.5)))
         // -ln(1-(1-(2^(-53)*0.5)))
@@ -518,12 +598,8 @@ impl WorldSim {
         // ((-ln(1-((1-2^(-10*2)))))/ln(1.002))/((-ln(1-((1 - 2^(-10*2)))))/ln(1+2^(-9)))
         // ((-ln(1-((2^(-10*2)))))/ln(1.002))/((-ln(1-((1 - 2^(-10*2)))))/ln(1+2^(-9)))
         // ((-ln(1-((1-2^(-10*2)))))/ln(1.002))/((-ln(1-((1 - 2^(-10*2)))))/ln(1.002))
-        let min_epsilon =
-            1.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64).max(f64::EPSILON as f64 * 0.5);
-        let max_epsilon = (1.0 - 1.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64))
-            .min(1.0 - f64::EPSILON as f64 * 0.5);
         // ((ln(0.6)-ln(1-0.6)) - (ln(1/(2048*2048))-ln((1-1/(2048*2048)))))/((ln(1-1/(2048*2048))-ln(1-(1-1/(2048*2048)))) - (ln(1/(2048*2048))-ln((1-1/(2048*2048)))))
-        let inv_func = /*|x: f64| x*//*exp_inverse_cdf*/logit/*hypsec_inverse_cdf*/;
+        let inv_func = |x: f64| x/*exp_inverse_cdf*//*logit*//*hypsec_inverse_cdf*/;
         let alt_exp_min_uniform = /*exp_inverse_cdf*//*logit*/inv_func(min_epsilon);
         let alt_exp_max_uniform = /*exp_inverse_cdf*//*logit*/inv_func(max_epsilon);
 
@@ -540,6 +616,82 @@ impl WorldSim {
         ).powf(erosion_pow).mul(0.5))*/; */
         let erosion_factor = |x: f64| (/*if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * */(/*exp_inverse_cdf*//*logit*/inv_func(x) - alt_exp_min_uniform) / (alt_exp_max_uniform - alt_exp_min_uniform))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
 ).powf(erosion_pow).mul(0.5))*//*.powf(0.5)*//*.powf(1.5)*//*.powf(2.0)*/;
+        let kf_func = {
+            |posi| {
+                if is_ocean_fn(posi) {
+                    return 1.0e-4;
+                }
+                let wposf = (uniform_idx_as_vec2(posi)
+                    * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
+                .map(|e| e as f64);
+                let turb_wposf =
+                    wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
+                let turb = Vec2::new(
+                    gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
+                    gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
+                ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+                // let turb = Vec2::zero();
+                let turb_wposf = wposf + turb;
+                let turb_wposi = turb_wposf
+                    .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64)
+                    .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0));
+                let turb_posi = vec2_as_uniform_idx(turb_wposi);
+                let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array())
+                    /* .min(0.5)
+                    .max(-0.5)*/
+                    .min(1.0)
+                    .max(-1.0)
+                    .mul(0.5)
+                    .add(0.5);
+                // kf = 1.5e-4: high-high (plateau [fan sediment])
+                // kf = 1e-4: high (plateau)
+                // kf = 2e-5: normal (dyke [unexposed])
+                // kf = 1e-6: normal-low (dyke [exposed])
+                // kf = 2e-6: low (mountain)
+                ((1.0 - uheight) * (1.5e-4 - 2.0e-6) + 2.0e-6) as f32
+            }
+        };
+        let kd_func = {
+            |posi| {
+                if is_ocean_fn(posi) {
+                    return 1.0e-2;
+                }
+                let wposf = (uniform_idx_as_vec2(posi)
+                    * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
+                .map(|e| e as f64);
+                let turb_wposf =
+                    wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
+                let turb = Vec2::new(
+                    gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
+                    gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
+                ) * /*64.0*/32.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+                // let turb = Vec2::zero();
+                let turb_wposf = wposf + turb;
+                let turb_wposi = turb_wposf
+                    .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64)
+                    .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0));
+                let turb_posi = vec2_as_uniform_idx(turb_wposi);
+                let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array())
+                    /* .min(0.5)
+                    .max(-0.5)*/
+                    .min(1.0)
+                    .max(-1.0)
+                    .mul(0.5)
+                    .add(0.5);
+                // kd = 1e-1: high (mountain, dyke)
+                // kd = 1.5e-2: normal-high (plateau [fan sediment])
+                // kd = 1e-2: normal (plateau)
+                1.0e-2
+                // (uheight * (1.0e-1 - 1.0e-2) + 1.0e-2) as f32
+            }
+        };
+        let g_func =
+            |posi| {
+                if /*is_ocean_fn(posi)*/map_edge_factor(posi) == 0.0 {
+                    return 0.0;
+                }
+                1.0
+            };
         let uplift_fn =
             |posi| {
                 if is_ocean_fn(posi) {
@@ -577,7 +729,7 @@ impl WorldSim {
                             .mul(0.045)*/)
                 };
                 let height =
-                    ((old_height_uniform(posi) - alt_old_min_uniform) as f64
+                    ((/*old_height_uniform*/uplift_uniform[posi].0/*1*/ - alt_old_min_uniform) as f64
                     / (alt_old_max_uniform - alt_old_min_uniform) as f64)
                     /*((old_height(posi) - alt_old_min) as f64
                     / (alt_old_max - alt_old_min) as f64)*/
@@ -619,8 +771,8 @@ impl WorldSim {
                 assert!(height <= 1.0);
                 // assert!(alt_main >= 0.0);
                 let (bump_factor, bump_max) = if
-                /*height < f32::EPSILON as f64 * 0.5/*false*/*/
-                true {
+                /*height < f32::EPSILON as f64 * 0.5*//*false*/
+                /*true*/false {
                     (
                         /*(alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64*/
                         (alt_main / CONFIG.mountain_scale as f64 * 128.0).mul(0.1).powf(1.2) * /*(1.0 / CONFIG.mountain_scale as f64)*/(f32::EPSILON * 0.5) as f64,
@@ -635,20 +787,34 @@ impl WorldSim {
                 // tan(pi/6)*32 ~ 18
                 // tan(54/360*2*pi)*32
                 // let height = 1.0f64;
-                let height = gen_ctx.uplift_nz.get(wposf.into_array())
-                    .min(0.5)
-                    .max(-0.5)
-                    .add(0.5);
+                /* let turb_wposf =
+                    wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(64.0);
+                let turb = Vec2::new(
+                    gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
+                    gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
+                ) * 64.0 * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+                let turb_wposf = wposf + turb;
+                let height = gen_ctx.uplift_nz.get(turb_wposf.into_array())
+                    /* .min(0.5)
+                    .max(-0.5)*/
+                    .min(1.0)
+                    .max(-1.0)
+                    .mul(0.5)
+                    .add(0.5); */
+
                 // let height = 1.0 / 7.0f64;
+                // let height = 0.0 / 31.0f64;
                 let bfrac = /*erosion_factor(0.5);*/0.0;
                 let height = (height - bfrac).abs().div(1.0 - bfrac);
                 let height = height
+                    .mul(31.0 / 32.0)
+                    .add(1.0 / 32.0)
                     /* .mul(15.0 / 16.0)
                     .add(1.0 / 16.0) */
                     /* .mul(5.0 / 8.0)
                     .add(3.0 / 8.0) */
-                    .mul(7.0 / 8.0)
-                    .add(1.0 / 8.0)
+                    /* .mul(7.0 / 8.0)
+                    .add(1.0 / 8.0) */
                     .mul(max_erosion_per_delta_t)
                     .sub(/*1.0 / CONFIG.mountain_scale as f64*/ bump_max)
                     .add(bump_factor);
@@ -698,7 +864,10 @@ impl WorldSim {
                         /* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
                             .mul(0.045)*/)
                 };
-                old_height_uniform(posi)/*.powf(2.0)*/
+
+                /* 0.0 */
+                old_height_uniform(posi)/*.powf(2.0)*/ - 0.5
+                // uplift_fn(posi) * (CONFIG.mountain_scale / max_erosion_per_delta_t as f32)
                 // 0.0
                 /* // 0.0
                 // -/*CONFIG.sea_level / CONFIG.mountain_scale*//* 0.75 */1.0
@@ -722,6 +891,9 @@ impl WorldSim {
             |posi| alt_func(posi),// if is_ocean_fn(posi) { old_height(posi) } else { 0.0 },
             is_ocean_fn,
             uplift_fn,
+            |posi| kf_func(posi),
+            |posi| kd_func(posi),
+            |posi| g_func(posi),
         );
         // Quick "small scale" erosion cycle in order to lower extreme angles.
         let (alt, basement) = do_erosion(
@@ -734,12 +906,15 @@ impl WorldSim {
             |posi| basement[posi],
             is_ocean_fn,
             |posi| uplift_fn(posi) * (1.0 * height_scale / max_erosion_per_delta_t) as f32,
+            |posi| kf_func(posi),
+            |posi| kd_func(posi),
+            |posi| g_func(posi),
         );
 
         let is_ocean = get_oceans(|posi| alt[posi]);
         let is_ocean_fn = |posi: usize| is_ocean[posi];
         let mut dh = downhill(&alt, /*old_height*/ is_ocean_fn);
-        let (boundary_len, indirection, water_alt_pos) = get_lakes(&/*water_alt*/alt, &mut dh);
+        let (boundary_len, indirection, water_alt_pos, _) = get_lakes(&/*water_alt*/alt, &mut dh);
         let flux_old = get_drainage(&water_alt_pos, &dh, boundary_len);
 
         let water_height_initial = |chunk_idx| {
@@ -787,6 +962,11 @@ impl WorldSim {
         };
 
         let water_alt = fill_sinks(water_height_initial, is_ocean_fn);
+        /* let water_alt = (0..WORLD_SIZE.x * WORLD_SIZE.y)
+            .into_par_iter()
+            .map(|posi| water_height_initial(posi))
+            .collect::<Vec<_>>(); */
+
         let rivers = get_rivers(&water_alt_pos, &water_alt, &dh, &indirection, &flux_old);
 
         let water_alt = indirection
@@ -939,7 +1119,9 @@ impl WorldSim {
             rng,
         };
 
-        this.seed_elements();
+        if opts.seed_elements {
+            this.seed_elements();
+        }
 
         this
     }
@@ -1565,7 +1747,7 @@ impl SimChunk {
             is_cliffs: cliff > 0.5 && !is_underwater,
             near_cliffs: cliff > 0.2,
             tree_density,
-            forest_kind: if temp > 0.0 {
+            forest_kind: if temp > CONFIG.temperate_temp {
                 if temp > CONFIG.desert_temp {
                     if humidity > CONFIG.jungle_hum {
                         // Forests in desert temperatures with extremely high humidity
diff --git a/world/src/sim/util.rs b/world/src/sim/util.rs
index b96ca14e99..046cb05581 100644
--- a/world/src/sim/util.rs
+++ b/world/src/sim/util.rs
@@ -281,13 +281,19 @@ pub fn get_oceans(oldh: impl Fn(usize) -> f32 + Sync) -> BitBox {
     // the sides are connected to it, and any subsequent ocean tiles must be connected to it.
     let mut is_ocean = bitbox![0; WORLD_SIZE.x * WORLD_SIZE.y];
     let mut stack = Vec::new();
+    let mut do_push = |pos| {
+        let posi = vec2_as_uniform_idx(pos);
+        if oldh(posi) <= 0.0 {
+            stack.push(posi);
+        }
+    };
     for x in 0..WORLD_SIZE.x as i32 {
-        stack.push(vec2_as_uniform_idx(Vec2::new(x, 0)));
-        stack.push(vec2_as_uniform_idx(Vec2::new(x, WORLD_SIZE.y as i32 - 1)));
+        do_push(Vec2::new(x, 0));
+        do_push(Vec2::new(x, WORLD_SIZE.y as i32 - 1));
     }
     for y in 1..WORLD_SIZE.y as i32 - 1 {
-        stack.push(vec2_as_uniform_idx(Vec2::new(0, y)));
-        stack.push(vec2_as_uniform_idx(Vec2::new(WORLD_SIZE.x as i32 - 1, y)));
+        do_push(Vec2::new(0, y));
+        do_push(Vec2::new(WORLD_SIZE.x as i32 - 1, y));
     }
     while let Some(chunk_idx) = stack.pop() {
         // println!("Ocean chunk {:?}: {:?}", uniform_idx_as_vec2(chunk_idx), oldh(chunk_idx));
@@ -532,6 +538,7 @@ impl NoiseFn<Point2<f64>> for HybridMulti {
         // Offset and bias to scale into [offset - 1.0, 1.0 + offset] range.
         let bias = 1.0;
         let mut result = (self.sources[0].get(point) + self.offset) * bias * self.persistence;
+        let mut exp_scale = 1.0;
         let mut scale = self.persistence;
         let mut weight = result;
 
@@ -547,7 +554,7 @@ impl NoiseFn<Point2<f64>> for HybridMulti {
             let mut signal = (self.sources[x].get(point) + self.offset) * bias;
 
             // Scale the amplitude appropriately for this frequency.
-            let exp_scale = self.persistence.powi(x as i32);
+            exp_scale *= self.persistence;
             signal *= exp_scale;
 
             // Add it in, weighted by previous octave's noise value.
@@ -571,6 +578,7 @@ impl NoiseFn<Point3<f64>> for HybridMulti {
         // Offset and bias to scale into [offset - 1.0, 1.0 + offset] range.
         let bias = 1.0;
         let mut result = (self.sources[0].get(point) + self.offset) * bias * self.persistence;
+        let mut exp_scale = 1.0;
         let mut scale = self.persistence;
         let mut weight = result;
 
@@ -586,7 +594,7 @@ impl NoiseFn<Point3<f64>> for HybridMulti {
             let mut signal = (self.sources[x].get(point) + self.offset) * bias;
 
             // Scale the amplitude appropriately for this frequency.
-            let exp_scale = self.persistence.powi(x as i32);
+            exp_scale *= self.persistence;
             signal *= exp_scale;
 
             // Add it in, weighted by previous octave's noise value.
@@ -610,7 +618,7 @@ impl NoiseFn<Point4<f64>> for HybridMulti {
         // Offset and bias to scale into [offset - 1.0, 1.0 + offset] range.
         let bias = 1.0;
         let mut result = (self.sources[0].get(point) + self.offset) * bias * self.persistence;
-        let mut exp_scale = self.persistence;
+        let mut exp_scale = 1.0;
         let mut scale = self.persistence;
         let mut weight = result;