diff --git a/Cargo.lock b/Cargo.lock
index 0b3576a..4ac4907 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -1,6 +1,6 @@
 # This file is automatically @generated by Cargo.
 # It is not intended for manual editing.
-version = 3
+version = 4
 
 [[package]]
 name = "addr2line"
@@ -35,6 +35,37 @@ version = "1.0.86"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "b3d1d046238990b9cf5bcde22a3fb3584ee5cf65fb2765f454ed428c7a0063da"
 
+[[package]]
+name = "argh"
+version = "0.1.14"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7f384d96bfd3c0b3c41f24dae69ee9602c091d64fc432225cf5295b5abbe0036"
+dependencies = [
+ "argh_derive",
+ "argh_shared",
+]
+
+[[package]]
+name = "argh_derive"
+version = "0.1.14"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "938e5f66269c1f168035e29ed3fb437b084e476465e9314a0328f4005d7be599"
+dependencies = [
+ "argh_shared",
+ "proc-macro2",
+ "quote",
+ "syn 2.0.87",
+]
+
+[[package]]
+name = "argh_shared"
+version = "0.1.14"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5127f8a5bc1cfb0faf1f6248491452b8a5b6901068d8da2d47cbb285986ae683"
+dependencies = [
+ "serde",
+]
+
 [[package]]
 name = "autocfg"
 version = "1.3.0"
@@ -114,6 +145,15 @@ version = "1.6.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "514de17de45fdb8dc022b1a7975556c53c86f9f0aa5f534b98977b171857c2c9"
 
+[[package]]
+name = "castaway"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dec551ab6e7578819132c713a93c022a05d60159dc86e7a7050223577484c55a"
+dependencies = [
+ "rustversion",
+]
+
 [[package]]
 name = "cc"
 version = "1.0.98"
@@ -222,6 +262,7 @@ name = "ewo3"
 version = "0.1.0"
 dependencies = [
  "anyhow",
+ "argh",
  "bincode",
  "euclid",
  "fastrand",
@@ -230,6 +271,8 @@ dependencies = [
  "image",
  "indexmap",
  "lazy_static",
+ "ndarray",
+ "ndarray-conv",
  "noise-functions",
  "rayon",
  "seahash",
@@ -287,7 +330,7 @@ checksum = "87750cf4b7a4c0625b1529e4c543c2182106e4dedc60a2a6455e00d212c489ac"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.66",
+ "syn 2.0.87",
 ]
 
 [[package]]
@@ -443,6 +486,16 @@ version = "0.4.21"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "90ed8c1e510134f979dbc4f070f87d4313098b704861a105fe34231c70a3901c"
 
+[[package]]
+name = "matrixmultiply"
+version = "0.3.10"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a06de3016e9fae57a36fd14dba131fccf49f74b40b7fbdb472f96e361ec71a08"
+dependencies = [
+ "autocfg",
+ "rawpointer",
+]
+
 [[package]]
 name = "memchr"
 version = "2.7.2"
@@ -470,12 +523,106 @@ dependencies = [
  "windows-sys 0.48.0",
 ]
 
+[[package]]
+name = "ndarray"
+version = "0.17.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "520080814a7a6b4a6e9070823bb24b4531daac8c4627e08ba5de8c5ef2f2752d"
+dependencies = [
+ "matrixmultiply",
+ "num-complex",
+ "num-integer",
+ "num-traits",
+ "portable-atomic",
+ "portable-atomic-util",
+ "rawpointer",
+ "rayon",
+]
+
+[[package]]
+name = "ndarray-conv"
+version = "0.6.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a502892118d1e27f9c514f61c44dd6fe46821ef53cb12b5af9043e8ca97a786b"
+dependencies = [
+ "castaway",
+ "ndarray",
+ "num",
+ "realfft",
+ "rustfft",
+ "thiserror 2.0.18",
+]
+
 [[package]]
 name = "noise-functions"
 version = "0.2.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "822a69eedf004ac2f492119af7a8203790b1c9115b9a9ef6bcd0cde5d6783565"
 
+[[package]]
+name = "num"
+version = "0.4.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "35bd024e8b2ff75562e5f34e7f4905839deb4b22955ef5e73d2fea1b9813cb23"
+dependencies = [
+ "num-bigint",
+ "num-complex",
+ "num-integer",
+ "num-iter",
+ "num-rational",
+ "num-traits",
+]
+
+[[package]]
+name = "num-bigint"
+version = "0.4.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a5e44f723f1133c9deac646763579fdb3ac745e418f2a7af9cd0c431da1f20b9"
+dependencies = [
+ "num-integer",
+ "num-traits",
+]
+
+[[package]]
+name = "num-complex"
+version = "0.4.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-integer"
+version = "0.1.46"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-iter"
+version = "0.1.45"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1429034a0490724d0075ebb2bc9e875d6503c3cf69e235a8941aa757d83ef5bf"
+dependencies = [
+ "autocfg",
+ "num-integer",
+ "num-traits",
+]
+
+[[package]]
+name = "num-rational"
+version = "0.4.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824"
+dependencies = [
+ "num-bigint",
+ "num-integer",
+ "num-traits",
+]
+
 [[package]]
 name = "num-traits"
 version = "0.2.19"
@@ -558,12 +705,36 @@ dependencies = [
  "miniz_oxide",
 ]
 
+[[package]]
+name = "portable-atomic"
+version = "1.13.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c33a9471896f1c69cecef8d20cbe2f7accd12527ce60845ff44c153bb2a21b49"
+
+[[package]]
+name = "portable-atomic-util"
+version = "0.2.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7a9db96d7fa8782dd8c15ce32ffe8680bbd1e978a43bf51a34d39483540495f5"
+dependencies = [
+ "portable-atomic",
+]
+
 [[package]]
 name = "ppv-lite86"
 version = "0.2.17"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "5b40af805b3121feab8a3c29f04d8ad262fa8e0561883e7653e024ae4479e6de"
 
+[[package]]
+name = "primal-check"
+version = "0.3.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dc0d895b311e3af9902528fbb8f928688abbd95872819320517cc24ca6b2bd08"
+dependencies = [
+ "num-integer",
+]
+
 [[package]]
 name = "proc-macro2"
 version = "1.0.85"
@@ -575,9 +746,9 @@ dependencies = [
 
 [[package]]
 name = "quote"
-version = "1.0.36"
+version = "1.0.44"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "0fa76aaf39101c457836aec0ce2316dbdc3ab723cdda1c6bd4e6ad4208acaca7"
+checksum = "21b2ebcf727b7760c461f091f9f0f539b77b8e87f2fd88131e7f1b433b3cece4"
 dependencies = [
  "proc-macro2",
 ]
@@ -612,6 +783,12 @@ dependencies = [
  "getrandom",
 ]
 
+[[package]]
+name = "rawpointer"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
+
 [[package]]
 name = "rayon"
 version = "1.10.0"
@@ -632,6 +809,15 @@ dependencies = [
  "crossbeam-utils",
 ]
 
+[[package]]
+name = "realfft"
+version = "3.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f821338fddb99d089116342c46e9f1fbf3828dba077674613e734e01d6ea8677"
+dependencies = [
+ "rustfft",
+]
+
 [[package]]
 name = "redox_syscall"
 version = "0.5.1"
@@ -647,6 +833,26 @@ version = "0.1.24"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "719b953e2095829ee67db738b3bfa9fa368c94900df327b3f07fe6e794d2fe1f"
 
+[[package]]
+name = "rustfft"
+version = "6.4.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "21db5f9893e91f41798c88680037dba611ca6674703c1a18601b01a72c8adb89"
+dependencies = [
+ "num-complex",
+ "num-integer",
+ "num-traits",
+ "primal-check",
+ "strength_reduce",
+ "transpose",
+]
+
+[[package]]
+name = "rustversion"
+version = "1.0.22"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b39cdef0fa800fc44525c84ccb54a029961a8215f9619753635a9c0d2538d46d"
+
 [[package]]
 name = "ryu"
 version = "1.0.18"
@@ -682,7 +888,7 @@ checksum = "500cbc0ebeb6f46627f50f3f5811ccf6bf00643be300b4c3eabc0ef55dc5b5ba"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.66",
+ "syn 2.0.87",
 ]
 
 [[package]]
@@ -753,6 +959,12 @@ version = "0.9.8"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "6980e8d7511241f8acf4aebddbb1ff938df5eebe98691418c4468d0b72a96a67"
 
+[[package]]
+name = "strength_reduce"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "fe895eb47f22e2ddd4dabc02bce419d2e643c8e3b585c78158b349195bc24d82"
+
 [[package]]
 name = "syn"
 version = "1.0.109"
@@ -766,9 +978,9 @@ dependencies = [
 
 [[package]]
 name = "syn"
-version = "2.0.66"
+version = "2.0.87"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c42f3f41a2de00b01c0aaad383c5a45241efc8b2d1eda5661812fda5f3cdcff5"
+checksum = "25aa4ce346d03a6dcd68dd8b4010bcb74e54e62c90c573f394c46eae99aba32d"
 dependencies = [
  "proc-macro2",
  "quote",
@@ -781,7 +993,16 @@ version = "1.0.61"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "c546c80d6be4bc6a00c0f01730c08df82eaa7a7a61f11d656526506112cc1709"
 dependencies = [
- "thiserror-impl",
+ "thiserror-impl 1.0.61",
+]
+
+[[package]]
+name = "thiserror"
+version = "2.0.18"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4288b5bcbc7920c07a1149a35cf9590a2aa808e0bc1eafaade0b80947865fbc4"
+dependencies = [
+ "thiserror-impl 2.0.18",
 ]
 
 [[package]]
@@ -792,7 +1013,18 @@ checksum = "46c3384250002a6d5af4d114f2845d37b57521033f30d5c3f46c4d70e1197533"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.66",
+ "syn 2.0.87",
+]
+
+[[package]]
+name = "thiserror-impl"
+version = "2.0.18"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ebc4ee7f67670e9b64d05fa4253e753e016c6c95ff35b89b7941d6b856dec1d5"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.87",
 ]
 
 [[package]]
@@ -833,7 +1065,7 @@ checksum = "5f5ae998a069d4b5aba8ee9dad856af7d520c3699e6159b185c2acd48155d39a"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.66",
+ "syn 2.0.87",
 ]
 
 [[package]]
@@ -848,6 +1080,16 @@ dependencies = [
  "tungstenite",
 ]
 
+[[package]]
+name = "transpose"
+version = "0.2.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1ad61aed86bc3faea4300c7aee358b4c6d0c8d6ccc36524c96e4c92ccf26e77e"
+dependencies = [
+ "num-integer",
+ "strength_reduce",
+]
+
 [[package]]
 name = "tungstenite"
 version = "0.23.0"
@@ -862,7 +1104,7 @@ dependencies = [
  "log",
  "rand",
  "sha1",
- "thiserror",
+ "thiserror 1.0.61",
  "utf-8",
 ]
 
@@ -1058,5 +1300,5 @@ checksum = "15e934569e47891f7d9411f1a451d947a60e000ab3bd24fbb970f000387d1b3b"
 dependencies = [
  "proc-macro2",
  "quote",
- "syn 2.0.66",
+ "syn 2.0.87",
 ]
diff --git a/Cargo.toml b/Cargo.toml
index c30c20b..45f6b74 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -23,10 +23,13 @@ image = { version = "0.25", default-features = false, features = ["png"] }
 rayon = "1"
 bincode = { version = "2.0.0-rc.3", features = ["serde"] }
 smallvec = "1"
+argh = "0.1"
+ndarray-conv = "0.6.0"
+ndarray = "0.17.2"
 
 [[bin]]
 name = "worldgen"
 path = "src/worldgen_test.rs"
 
 [profile.release]
-debug = true
\ No newline at end of file
+debug = true
diff --git a/src/main.rs b/src/main.rs
index 6913fe2..b51d4ce 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -1,3 +1,6 @@
+#![feature(test)]
+extern crate test;
+
 use hecs::{CommandBuffer, Entity, With, World};
 use futures_util::{stream::TryStreamExt, SinkExt, StreamExt};
 use indexmap::IndexMap;
@@ -88,9 +91,25 @@ struct GameState {
     clients: Slab<Client>,
     ticks: u64,
     map: worldgen::GeneratedWorld,
+    baseline_soil_nutrients: Map<f32>,
+    baseline_water: Map<f32>,
+    baseline_salt: Map<f32>,
+    baseline_temperature: Map<f32>,
+    dynamic_soil_nutrients: Map<f32>,
+    dynamic_water: Map<f32>,
     positions: PositionIndex
 }
 
+impl GameState {
+    fn actual_water(&self, pos: Coord) -> f32 {
+        self.baseline_water[pos] + self.dynamic_water[pos]
+    }
+
+    fn actual_soil_nutrients(&self, pos: Coord) -> f32 {
+        self.baseline_soil_nutrients[pos] + self.dynamic_soil_nutrients[pos]
+    }
+}
+
 #[derive(Debug, Clone, PartialEq, Eq, Hash)]
 enum Item {
     Dirt,
@@ -264,7 +283,10 @@ struct DespawnOnImpact;
 struct Inventory(indexmap::IndexMap<Item, u64>);
 
 #[derive(Debug, Clone)]
-struct Plant(plant::Genome);
+struct Plant {
+    genome: plant::Genome,
+    current_size: f32
+}
 
 #[derive(Debug, Clone)]
 struct NewlyAdded; // ugly hack to work around ECS deficiencies
@@ -381,6 +403,10 @@ impl StochasticNumber {
     }
 }
 
+const PLANT_TICK_DELAY: u64 = 128;
+const FIELD_DECAY_DELAY: u64 = 100;
+const PLANT_GROWTH_SCALE: f32 = 0.01;
+
 async fn game_tick(state: &mut GameState) -> Result<()> {
     let mut buffer = hecs::CommandBuffer::new();
 
@@ -391,6 +417,16 @@ async fn game_tick(state: &mut GameState) -> Result<()> {
 
     buffer.run_on(&mut state.world);
 
+    if state.ticks % FIELD_DECAY_DELAY == 0 {
+        state.baseline_soil_nutrients.for_each_mut(|nutrients| *nutrients *= 0.999);
+    } else if state.ticks % FIELD_DECAY_DELAY == 1 {
+        state.baseline_water.for_each_mut(|water| *water *= 0.999);
+    } else if state.ticks % FIELD_DECAY_DELAY == 2 {
+        state.dynamic_soil_nutrients = smooth(&state.dynamic_soil_nutrients, 3);
+    } else if state.ticks % FIELD_DECAY_DELAY == 3 {
+        state.dynamic_water = smooth(&state.dynamic_water, 3);
+    }
+
     // Spawn enemies
     for (_entity, (pos, EnemyTarget { spawn_range, spawn_density, spawn_rate_inv, .. })) in state.world.query::<(&Position, &EnemyTarget)>().iter() {
         let pos = pos.head();
@@ -445,6 +481,22 @@ async fn game_tick(state: &mut GameState) -> Result<()> {
         }
     }
 
+    // Run plant simulations.
+    for (entity, (pos, plant)) in state.world.query::<(&Position, &mut Plant)>().iter() {
+        if (entity.id() as u64) % PLANT_TICK_DELAY == state.ticks % PLANT_TICK_DELAY {
+            let pos = pos.head();
+            let water = state.actual_soil_nutrients(pos);
+            let soil_nutrients = state.actual_soil_nutrients(pos);
+            let salt = state.baseline_salt[pos];
+            let temperature = state.baseline_temperature[pos];
+            let base_growth_rate = plant.genome.effective_growth_rate(soil_nutrients, water, temperature, salt);
+            plant.current_size += base_growth_rate * PLANT_GROWTH_SCALE * plant.current_size.powf(-0.25); // allometric scaling law
+            plant.current_size = plant.current_size.min(plant.genome.max_size);
+            let can_reproduce = plant.current_size >= plant.genome.max_size * plant.genome.reproductive_size_fraction();
+            //state.dynamic_soil_nutrients;
+        }
+    }
+
     // Process enemy motion and ranged attacks
     for (entity, (pos, ranged, energy, jump)) in state.world.query::<hecs::With<(&Position, Option<&mut RangedAttack>, Option<&mut Energy>, Option<&Jump>), &Enemy>>().iter() {
         let pos = pos.head();
@@ -826,12 +878,25 @@ async fn main() -> Result<()> {
         }
     };
 
+    let baseline_soil_nutrients = world.soil_nutrients.clone();
+    let baseline_water = world.groundwater.clone();
+    let baseline_salt = world.salt.clone();
+    let baseline_temperature = world.temperature.clone();
+    let dynamic_soil_nutrients = baseline_soil_nutrients.clone();
+    let dynamic_water = baseline_water.clone();
+
     let state = Arc::new(Mutex::new(GameState {
         world: World::new(),
         clients: Slab::new(),
         ticks: 0,
         positions: PositionIndex::new(world.radius),
-        map: world
+        map: world,
+        baseline_soil_nutrients,
+        baseline_water,
+        baseline_salt,
+        baseline_temperature,
+        dynamic_soil_nutrients,
+        dynamic_water
     }));
 
     {
@@ -844,7 +909,7 @@ async fn main() -> Result<()> {
                 Render('+'),
                 Health(100.0, 100.0),
                 //ShrinkOnDeath,
-                Plant(plant::Genome::random()),
+                Plant { genome: plant::Genome::random(), current_size: 0.0 },
                 NewlyAdded
             ));
         }
@@ -862,9 +927,11 @@ async fn main() -> Result<()> {
         interval.set_missed_tick_behavior(tokio::time::MissedTickBehavior::Skip);
         loop {
             let mut state = state.lock().await;
+            let time = std::time::Instant::now();
             if let Err(e) = game_tick(&mut state).await {
                 println!("{:?}", e);
             }
+            println!("Tick time: {:?}", time.elapsed());
             interval.tick().await;
         }
     });
diff --git a/src/map.rs b/src/map.rs
index 1ee5441..0c6f1e3 100644
--- a/src/map.rs
+++ b/src/map.rs
@@ -1,6 +1,9 @@
 use euclid::{Point3D, Point2D, Vector2D};
 use serde::{Deserialize, Serialize};
 use std::ops::{Index, IndexMut};
+use std::marker::PhantomData;
+use ndarray::prelude::*;
+use ndarray_conv::{ConvExt, ConvMode, PaddingMode};
 
 pub struct AxialWorldSpace;
 pub struct CubicWorldSpace;
@@ -60,22 +63,30 @@ pub fn count_hexes(x: i32) -> i32 {
     x*(x+1)*3+1
 }
 
-struct CoordsIndexIterator {
+struct CoordIterator {
     radius: i32,
-    index: usize,
+    count: usize,
     r: i32,
     q: i32,
     max: usize
 }
 
-impl Iterator for CoordsIndexIterator {
-    type Item = (Coord, usize);
+struct CoordIterMut<'a, T> {
+    coords: CoordIterator,
+    data: *mut T,
+    radius: i32,
+    side_length: usize,
+    _marker: PhantomData<&'a mut T>
+}
+
+impl Iterator for CoordIterator {
+    type Item = Coord;
     fn next(&mut self) -> Option<Self::Item> {
-        if self.index == self.max {
+        if self.count == self.max {
             return None;
         }
-        let result = (Coord::new(self.q, self.r), self.index);
-        self.index += 1;
+        let result = Coord::new(self.q, self.r);
+        self.count += 1;
         self.q += 1;
         if self.r < 0 && self.q == self.radius + 1 {
             self.r += 1;
@@ -89,81 +100,143 @@ impl Iterator for CoordsIndexIterator {
     }
 }
 
+// blame OpenAI for this, and also Mozilla
+impl<'a, T> Iterator for CoordIterMut<'a, T> {
+    type Item = (Coord, &'a mut T);
+
+    fn next(&mut self) -> Option<Self::Item> {
+        let coord = self.coords.next()?;
+        let q = (coord.x + self.radius) as usize;
+        let r = (coord.y + self.radius) as usize;
+        let i = q + r * self.side_length;
+        unsafe {
+            // CoordIterator yields each valid map coordinate once, so each backing index is yielded once.
+            Some((coord, &mut *self.data.add(i)))
+        }
+    }
+}
+
 #[derive(Debug, Clone, Serialize, Deserialize)]
 pub struct Map<T> {
-    pub data: Vec<T>,
-    pub radius: i32
+    data: Vec<T>,
+    pub radius: i32,
+    side_length: usize
 }
 
 impl<T> Map<T> {
     pub fn new(radius: i32, fill: T) -> Map<T> where T: Clone {
-        let size = count_hexes(radius) as usize;
+        // We represent worlds using axial coordinates (+q is right, +r is down-right).
+        // This results in a parallelogram shape.
+        // However, the map is a hexagon for the purposes of gameplay, so the top-left and bottom-right corners are invalid.
+        let side_length = (radius * 2 + 1) as usize;
+        let size = side_length.pow(2);
         Map {
             data: vec![fill; size],
-            radius
-        }
-    }
-
-    pub fn from_fn<S, F: FnMut(Coord) -> S>(mut f: F, radius: i32) -> Map<S> {
-        let size = count_hexes(radius) as usize;
-        Map {
             radius,
-            data: Vec::from_iter(CoordsIndexIterator {
-                radius,
-                index: 0,
-                max: size,
-                r: -radius,
-                q: 0
-            }.map(|(c, _i)| f(c)))
+            side_length
         }
     }
 
-    pub fn map<S, F: FnMut(&T) -> S>(mut f: F, other: &Self) -> Map<S> {
+    pub fn size(&self) -> usize {
+        count_hexes(self.radius) as usize
+    }
+
+    pub fn from_fn<S: Default + Clone, F: FnMut(Coord) -> S>(mut f: F, radius: i32) -> Map<S> {
+        let mut map = Map::new(radius, S::default());
+        for coord in map.iter_coords() {
+            map[coord] = f(coord);
+        }
+        map
+    }
+
+    pub fn map<S: Default + Clone, F: FnMut(&T) -> S>(mut f: F, other: &Self) -> Map<S> {
         Map::<S>::from_fn(|c| f(&other[c]), other.radius)
     }
 
     pub fn coord_to_index(&self, c: Coord) -> usize {
-        let r = c.y + self.radius;
-        let fh = r.min(self.radius);
-        let mut coords_above = fh*(self.radius+1) + fh*(fh-1)/2;
-        if fh < r {
-            let d = r - fh;
-            coords_above += d*(2*self.radius+1) - d*(d-1)/2;
-        }
-        let q_start = if r < self.radius { -r } else { -self.radius };
-        (coords_above + (c.x - q_start)) as usize
+        let q = (c.x + self.radius) as usize;
+        let r = (c.y + self.radius) as usize;
+        q * self.side_length + r
     }
 
     pub fn in_range(&self, coord: Coord) -> bool {
         hex_distance(coord, Coord::origin()) <= self.radius
     }
 
-    pub fn iter_coords(&self) -> impl Iterator<Item=(Coord, usize)> {
-        CoordsIndexIterator {
+    pub fn iter_coords(&self) -> impl Iterator<Item=Coord> {
+        CoordIterator {
             radius: self.radius,
-            index: 0,
-            max: self.data.len(),
+            count: 0,
+            max: self.size(),
             r: -self.radius,
             q: 0
         }
     }
 
     pub fn iter(&self) -> impl Iterator<Item=(Coord, &T)> {
-        self.iter_coords().map(|(c, i)| (c, &self.data[i]))
+        self.iter_coords().map(|c| (c, &self[c]))
+    }
+
+    pub fn iter_data(&self) -> impl Iterator<Item=&T> {
+        self.iter_coords().map(|c| &self[c])
+    }
+
+    pub fn iter_mut(&mut self) -> impl Iterator<Item=(Coord, &mut T)> {
+        CoordIterMut {
+            coords: CoordIterator {
+                radius: self.radius,
+                count: 0,
+                max: self.size(),
+                r: -self.radius,
+                q: 0
+            },
+            data: self.data.as_mut_ptr(),
+            radius: self.radius,
+            side_length: self.side_length,
+            _marker: PhantomData
+        }
+    }
+
+    pub fn for_each_mut(&mut self, mut f: impl FnMut(&mut T)) {
+        for (_, value) in self.iter_mut() {
+            f(value);
+        }
+    }
+}
+
+// 2D hex convolution
+pub fn smooth(map: &Map<f32>, radius: i32) -> Map<f32> {
+    //let mut map = map.clone();
+    //map[Coord::new(1, 0)] = 3.0;
+    let data: ArrayBase<ndarray::ViewRepr<&f32>, Dim<[usize; 2]>, f32> = ArrayView2::from_shape((map.side_length, map.side_length), map.data.as_slice()).unwrap();
+    //println!("{:?}", data[(map.radius as usize + 1, map.radius as usize)]);
+    let mut kernel = Array2::from_elem((radius as usize * 2 + 1, radius as usize * 2 + 1), 0.0);
+
+    for (_, offset) in hex_range(radius) {
+        kernel[(offset.x as usize + radius as usize, offset.y as usize + radius as usize)] = 1.0 / count_hexes(radius) as f32;
+    }
+
+    let result = ConvExt::conv(&data, &kernel, ConvMode::Same, PaddingMode::Zeros).unwrap();
+
+    Map {
+        radius: map.radius,
+        side_length: map.side_length,
+        data: result.into_raw_vec(), // TODO fix
     }
 }
 
 impl<T> Index<Coord> for Map<T> {
     type Output = T;
     fn index(&self, index: Coord) -> &Self::Output {
-        //println!("{:?}", index);
+        debug_assert!(self.in_range(index), "invalid coord: {:?}", index);
         &self.data[self.coord_to_index(index)]
     }
 }
 
 impl<T> IndexMut<Coord> for Map<T> {
     fn index_mut(&mut self, index: Coord) -> &mut Self::Output {
+        debug_assert!(self.in_range(index), "invalid coord: {:?}", index);
         let i = self.coord_to_index(index);
         &mut self.data[i]
     }
-}
\ No newline at end of file
+}
diff --git a/src/plant.rs b/src/plant.rs
index b902d4c..98aac81 100644
--- a/src/plant.rs
+++ b/src/plant.rs
@@ -9,15 +9,17 @@ pub enum CropType {
 #[derive(Debug, Clone, Copy)]
 pub struct Genome {
     crop_type: CropType,
-    // polygenic traits; parameterized as N(0,1)
+    // polygenic traits; parameterized as N(0,1) (allegedly)
     growth_rate: f32,
-    nitrogen_fixation_rate: f32,
+    nutrient_addition_rate: f32,
     optimal_water_level: f32,
     optimal_temperature: f32,
     reproduction_rate: f32,
+    reproductive_size_fraction_param: f32,
     temperature_tolerance: f32,
     water_tolerance: f32,
-    max_size: f32
+    salt_tolerance: f32,
+    pub max_size: f32
     // TODO color trait
 }
 
@@ -37,16 +39,24 @@ fn normal_scaled(mu: f32, sigma: f32) -> f32 {
 }
 
 impl Genome {
-    pub fn effective_growth_rate(&self, water: f32, temperature: f32) -> f32 {
+    pub fn reproductive_size_fraction(&self) -> f32 {
+        sigmoid(self.reproductive_size_fraction_param)
+    }
+
+    pub fn effective_growth_rate(&self, nutrients: f32, water: f32, temperature: f32, salt: f32) -> f32 {
         let water_diff = (water - self.optimal_water_level).abs();
         let temperature_diff = (temperature - self.optimal_temperature).abs();
-        1.5f32.powf(self.growth_rate)
+        let salt_excess = (salt - sigmoid(self.salt_tolerance)).max(0.0);
+        let base = 1.5f32.powf(self.growth_rate)
             - self.reproduction_rate * 0.1 // faster reproduction trades off slightly against growth
-            - self.nitrogen_fixation_rate.max(0.0) * 0.16 // same for nitrogen fixation
+            - self.nutrient_addition_rate.max(0.0) * 0.16 // nutrient enrichment has a growth tradeoff
             - (water_diff - sigmoid(self.water_tolerance)).max(0.0) // penalize plants when far from optimal environmental range
             - (temperature_diff - sigmoid(self.temperature_tolerance)).max(0.0) // same for temperature
+            - salt_excess
             - self.water_tolerance * 0.2
             - self.temperature_tolerance * 0.2
+            - self.salt_tolerance * 0.2;
+        base * (-nutrients.min(0.0)).exp()
     }
 
     pub fn random() -> Genome {
@@ -58,23 +68,25 @@ impl Genome {
             _ => unreachable!()
         };
 
-        let (nitrogen_fixation_rate, optimal_water_level, optimal_temperature, max_size) = match crop_type {
-            CropType::Grass => (-10.0, 0.0, 0.0, 0.0),
-            CropType::EucalyptusTree => (-10.0, 2.0, 1.0, 5.0),
-            CropType::BushTomato => (-10.0, -1.0, 1.5, 1.0),
-            CropType::GoldenWattleTree => (2.0, 1.5, 1.0, 3.0),
+        let (nutrient_addition_rate, optimal_water_level, optimal_temperature, reproductive_size_fraction_param, salt_tolerance, max_size) = match crop_type {
+            CropType::Grass => (-10.0, 0.0, 0.0, -1.0, 0.2, 0.0),
+            CropType::EucalyptusTree => (-10.0, 2.0, 1.0, 1.0, 1.3, 5.0),
+            CropType::BushTomato => (-10.0, -1.0, 1.5, -0.3, 1.6, 1.0),
+            CropType::GoldenWattleTree => (2.0, 1.5, 1.0, 0.5, 0.9, 3.0),
 
         };
 
         Genome {
             crop_type: crop_type,
             growth_rate: normal(),
-            nitrogen_fixation_rate,
+            nutrient_addition_rate,
             optimal_water_level,
             optimal_temperature,
             reproduction_rate: normal(),
+            reproductive_size_fraction_param: normal() + reproductive_size_fraction_param,
             temperature_tolerance: normal(),
             water_tolerance: normal(),
+            salt_tolerance: normal() + salt_tolerance,
             max_size
         }
     }
@@ -84,12 +96,14 @@ impl Genome {
         Some(Genome {
             crop_type: self.crop_type,
             growth_rate: (self.growth_rate + other.growth_rate) / 2.0 + normal_scaled(0.0, 0.1),
-            nitrogen_fixation_rate: (self.nitrogen_fixation_rate + other.nitrogen_fixation_rate) / 2.0 + normal_scaled(0.0, 0.03),
+            nutrient_addition_rate: (self.nutrient_addition_rate + other.nutrient_addition_rate) / 2.0 + normal_scaled(0.0, 0.03),
             optimal_water_level: (self.optimal_water_level + other.optimal_water_level) / 2.0 + normal_scaled(0.0, 0.03),
             optimal_temperature: (self.optimal_temperature + other.optimal_temperature) / 2.0 + normal_scaled(0.0, 0.03),
             reproduction_rate: (self.reproduction_rate + other.reproduction_rate) / 2.0 + normal_scaled(0.0, 0.5),
+            reproductive_size_fraction_param: (self.reproductive_size_fraction_param + other.reproductive_size_fraction_param) / 2.0 + normal_scaled(0.0, 0.2),
             temperature_tolerance: (self.temperature_tolerance + other.temperature_tolerance) / 2.0 + normal_scaled(0.0, 0.2),
             water_tolerance: (self.water_tolerance + other.water_tolerance) / 2.0 + normal_scaled(0.0, 0.2),
+            salt_tolerance: (self.salt_tolerance + other.salt_tolerance) / 2.0 + normal_scaled(0.0, 0.2),
             max_size: (self.max_size + other.max_size) / 2.0 + normal_scaled(0.0, 0.3)
         })
     }
diff --git a/src/worldgen.rs b/src/worldgen.rs
index 18ce10d..1c2fb92 100644
--- a/src/worldgen.rs
+++ b/src/worldgen.rs
@@ -1,6 +1,7 @@
 use std::{cmp::Ordering, collections::{hash_map::Entry, BinaryHeap, HashMap, HashSet, VecDeque}, hash::{Hash, Hasher}};
 
 use noise_functions::Sample3;
+use rayon::prelude::*;
 use serde::{Deserialize, Serialize};
 use crate::map::*;
 
@@ -9,38 +10,46 @@ const NOISE_SCALE: f32 = 0.0005;
 const HEIGHT_EXPONENT: f32 = 0.3;
 const WATER_SOURCES: usize = 40;
 
-pub fn height_baseline(pos: Coord) -> f32 {
-    let w_frac = (hex_distance(pos, Coord::origin()) as f32 / WORLD_RADIUS as f32).powf(3.0);
+pub fn height_baseline_with_radius(pos: Coord, radius: i32) -> f32 {
+    let w_frac = (hex_distance(pos, Coord::origin()) as f32 / radius as f32).powf(3.0);
     let pos = to_cubic(pos);
-    let noise = 
+    let noise =
         noise_functions::OpenSimplex2s.ridged(6, 0.7, 1.55).seed(406).sample3([10.0 + pos.x as f32 * NOISE_SCALE, pos.y as f32 * NOISE_SCALE, pos.z as f32 * NOISE_SCALE]);
     let range = 1.0 - 2.0 * (w_frac - 0.5).powf(2.0);
     noise * range - w_frac
 }
 
+pub fn height_baseline(pos: Coord) -> f32 {
+    height_baseline_with_radius(pos, WORLD_RADIUS)
+}
+
 fn percentilize<F: Fn(f32) -> f32>(raw: &mut Map<f32>, postprocess: F) {
-    let mut xs: Vec<(usize, f32)> = raw.data.iter().copied().enumerate().collect();
+    let mut xs: Vec<(Coord, f32)> = raw.iter().map(|(coord, value)| (coord, *value)).collect();
     xs.sort_unstable_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
     for (j, (i, _x)) in xs.into_iter().enumerate() {
-        let percentile = j as f32 / raw.data.len() as f32;
-        raw.data[i] = postprocess(percentile);
+        let percentile = j as f32 / raw.size() as f32;
+        raw[i] = postprocess(percentile);
     }
 }
 
-fn normalize<F: Fn(f32) -> f32>(raw: &mut Map<f32>, postprocess: F) {
-    let mut min = raw.data.iter().copied().min_by(|a, b| a.partial_cmp(b).unwrap()).unwrap();
-    let mut max = raw.data.iter().copied().max_by(|a, b| a.partial_cmp(b).unwrap()).unwrap();
+fn normalize<F: Fn(f32) -> f32 + Sync>(raw: &mut Map<f32>, postprocess: F) {
+    let mut min = raw.iter_data().copied().min_by(|a, b| a.partial_cmp(b).unwrap()).unwrap();
+    let mut max = raw.iter_data().copied().max_by(|a, b| a.partial_cmp(b).unwrap()).unwrap();
     if min == max {
         min = 0.0;
         max = 1.0;
     }
-    for x in raw.data.iter_mut() {
+    for (_c, x) in raw.iter_mut() {
         *x = postprocess((*x - min) / (max - min));
     }
 }
 
 pub fn generate_heights() -> Map<f32> {
-    let mut raw = Map::<f32>::from_fn(height_baseline, WORLD_RADIUS);
+    generate_heights_with_radius(WORLD_RADIUS)
+}
+
+pub fn generate_heights_with_radius(radius: i32) -> Map<f32> {
+    let mut raw = Map::<f32>::from_fn(|pos| height_baseline_with_radius(pos, radius), radius);
     percentilize(&mut raw, |x| {
         let s = 1.0 - (1.0 - x).powf(HEIGHT_EXPONENT);
         s * 2.0 - 1.0
@@ -49,30 +58,35 @@ pub fn generate_heights() -> Map<f32> {
 }
 
 pub fn generate_contours(field: &Map<f32>, interval: f32) -> Vec<(Coord, f32, f32, CoordVec)> {
-    let mut v = vec![];
     // Starting at the origin, we want to detect contour lines in any of the six directions.
     // Go in one of the perpendicular directions to generate base directions then scan until the edge is reached starting from any of those points.
-    for scan_direction in DIRECTIONS {
-        // "perpendicular" doesn't really work in axial coordinates, but this apparently does so whatever
-        let slit_direction = rotate_60(*scan_direction);
-        for x in 0..field.radius {
-            let base = Coord::zero() + slit_direction * x;
-            let mut last: Option<f32> = None;
-            for y in 0..field.radius {
-                let point = base + *scan_direction * y;
-                if hex_distance(point, Coord::zero()) <= field.radius {
-                    let sample = field[point] / interval;
-                    if let Some(last) = last {
-                        if last.trunc() != sample.trunc() {
-                            v.push((point, last, sample, *scan_direction));
+    DIRECTIONS.par_iter()
+        .map(|scan_direction| {
+            let mut out = vec![];
+            // "perpendicular" doesn't really work in axial coordinates, but this apparently does so whatever
+            let slit_direction = rotate_60(*scan_direction);
+            for x in 0..field.radius {
+                let base = Coord::zero() + slit_direction * x;
+                let mut last: Option<f32> = None;
+                for y in 0..field.radius {
+                    let point = base + *scan_direction * y;
+                    if hex_distance(point, Coord::zero()) <= field.radius {
+                        let sample = field[point] / interval;
+                        if let Some(last) = last {
+                            if last.trunc() != sample.trunc() {
+                                out.push((point, last, sample, *scan_direction));
+                            }
                         }
+                        last = Some(sample);
                     }
-                    last = Some(sample);
                 }
             }
-        }
-    }
-    v
+            out
+        })
+        .reduce(Vec::new, |mut acc, mut chunk| {
+            acc.append(&mut chunk);
+            acc
+        })
 }
 
 #[derive(Clone, Copy, Debug)]
@@ -109,7 +123,7 @@ impl<T: Hash + Eq + PartialEq> PartialOrd for PointWrapper<T> {
 
 pub fn generate_separated_high_points(n: usize, sep: i32, map: &Map<f32>) -> Vec<Coord> {
     let mut points = vec![];
-    let mut priority = BinaryHeap::with_capacity(map.data.len());
+    let mut priority = BinaryHeap::with_capacity(map.size());
     for (coord, height) in map.iter() {
         priority.push(PointWrapper(-*height, coord));
     }
@@ -192,12 +206,13 @@ pub fn distance_map<I: Iterator<Item=Coord>>(radius: i32, sources: I) -> Map<f32
 
 pub fn compute_groundwater(water: &Map<f32>, rain: &Map<f32>, heightmap: &Map<f32>) -> Map<f32> {
     let mut groundwater = distance_map(
-        water.radius, 
+        water.radius,
         water.iter().filter_map(|(c, i)| if *i > 0.0 { Some(c) } else { None }));
     percentilize(&mut groundwater, |x| (1.0 - x).powf(0.3));
-    for (coord, h) in heightmap.iter() {
-        groundwater[coord] -= *h * 0.05;
-        groundwater[coord] += rain[coord] * 0.15;
+
+    for (coord, gw) in groundwater.iter_mut() {
+        *gw -= heightmap[coord] * 0.05;
+        *gw += rain[coord] * 0.15;
     }
     percentilize(&mut groundwater, |x| x.powf(0.7));
     groundwater
@@ -225,18 +240,18 @@ fn floodfill(src: Coord, all: &HashSet<Coord>) -> HashSet<Coord> {
 }
 
 pub fn get_sea(heightmap: &Map<f32>) -> (HashSet<Coord>, HashSet<Coord>) {
-    let sinks = heightmap.iter_coords().filter(|(c, _)| heightmap[*c] <= SEA_LEVEL).map(|(c, _)| c).collect::<HashSet<_>>();
-    let sea = floodfill(Coord::new(0, WORLD_RADIUS), &sinks);
+    let sinks = heightmap.iter_coords().filter(|c| heightmap[*c] <= SEA_LEVEL).collect::<HashSet<_>>();
+    let sea = floodfill(Coord::new(0, heightmap.radius), &sinks);
     (sinks, sea)
 }
 
 const SALT_REMOVAL: f32 = 0.13;
 const SALT_RANGE: f32 = 0.33;
 
-pub fn simulate_water(heightmap: &mut Map<f32>, rain_map: &Map<f32>, sea: &HashSet<Coord>, sinks: &HashSet<Coord>) -> (Map<f32>, Map<f32>) {  
+pub fn simulate_water(heightmap: &mut Map<f32>, rain_map: &Map<f32>, sea: &HashSet<Coord>, sinks: &HashSet<Coord>) -> (Map<f32>, Map<f32>) {
     let mut watermap = Map::<f32>::new(heightmap.radius, 0.0);
 
-    let sources = generate_separated_high_points(WATER_SOURCES, WORLD_RADIUS / 10, &rain_map);
+    let sources = generate_separated_high_points(WATER_SOURCES, (heightmap.radius / 10).max(1), rain_map);
     let mut remainder = sinks.clone();
 
     for s in sea.iter() {
@@ -258,13 +273,12 @@ pub fn simulate_water(heightmap: &mut Map<f32>, rain_map: &Map<f32>, sea: &HashS
         }
     }
 
-    let mut salt = distance_map(watermap.radius, sea.iter().copied());
+    let mut salt = distance_map(heightmap.radius, sea.iter().copied());
     normalize(&mut salt, |x| (SALT_RANGE - x).max(0.0) / SALT_RANGE);
-
-    for (coord, rain) in rain_map.iter() {
-        if *rain > 0.0 {
-            salt[coord] -= *rain * 0.3;
-            salt[coord] = salt[coord].max(0.0);
+    for (coord, salt) in salt.iter_mut() {
+        let rain = rain_map[coord];
+        if rain > 0.0 {
+            *salt = (*salt - rain * 0.3).max(0.0);
         }
     }
 
@@ -334,9 +348,13 @@ const NUTRIENT_NOISE_SCALE: f32 = 0.0015;
 // As a handwave, define soil nutrients to be partly randomized and partly based on water.
 // This kind of sort of makes sense because nitrogen is partly fixed by plants, which would have grown in water-having areas.
 pub fn soil_nutrients(groundwater: &Map<f32>) -> Map<f32> {
-    let mut soil_nutrients = Map::<f32>::from_fn(|cr| {
-        let c = to_cubic(cr);
-        noise_functions::OpenSimplex2s.seed(406).sample3([10.0 + c.x as f32 * NUTRIENT_NOISE_SCALE, c.y as f32 * NUTRIENT_NOISE_SCALE, c.z as f32 * NUTRIENT_NOISE_SCALE]) + groundwater[cr]
+    let mut soil_nutrients = Map::<f32>::from_fn(|d| {
+        let c = to_cubic(d);
+        noise_functions::OpenSimplex2s.seed(406).sample3([
+            10.0 + c.x as f32 * NUTRIENT_NOISE_SCALE,
+            c.y as f32 * NUTRIENT_NOISE_SCALE,
+            c.z as f32 * NUTRIENT_NOISE_SCALE
+        ]) + groundwater[d]
     }, groundwater.radius);
     percentilize(&mut soil_nutrients, |x| x.powf(0.4));
     soil_nutrients
@@ -356,22 +374,6 @@ fn august_roche_magnus(temperature: f32) -> f32 {
     6.1094 * f32::exp((17.625 * temperature) / (243.04 + temperature))
 }
 
-// 2D hex convolution
-fn smooth(map: &Map<f32>, radius: i32) -> Map<f32> {
-    let mut out = Map::<f32>::new(map.radius, 0.0);
-    for (coord, index) in map.iter_coords() {
-        let mut sum = map.data[index];
-        for (_, offset) in hex_range(radius) {
-            let neighbor = coord + offset;
-            if map.in_range(neighbor) {
-                sum += map[neighbor];
-            }
-        }
-        out.data[index] = sum / (count_hexes(radius) as f32);
-    }
-    out
-}
-
 const BASE_TEMPERATURE: f32 = 30.0; // degrees
 const HEIGHT_SCALE: f32 = 1e3; // unrealistic but makes world more interesting; m
 //const SEA_LEVEL_AIR_PRESSURE: f32 = 1013.0; // hPa
@@ -379,12 +381,13 @@ const HEIGHT_SCALE: f32 = 1e3; // unrealistic but makes world more interesting;
 const AIR_SPECIFIC_HEAT_CAPACITY: f32 = 1012.0; // J kg^-1 K^-1
 const EARTH_GRAVITY: f32 = 9.81; // m s^-2
 
-pub fn simulate_air(heightmap: &Map<f32>, sea: &HashSet<Coord>, scan_dir: CoordVec, perpendicular_dir: CoordVec) -> (Map<f32>, Map<f32>, Map<f32>) {   
-    let start_pos = Coord::origin() + -scan_dir * WORLD_RADIUS;
+pub fn simulate_air(heightmap: &Map<f32>, sea: &HashSet<Coord>, scan_dir: CoordVec, perpendicular_dir: CoordVec) -> (Map<f32>, Map<f32>, Map<f32>) {
+    let radius = heightmap.radius;
+    let start_pos = Coord::origin() + -scan_dir * radius;
     let mut rain_map = Map::<f32>::new(heightmap.radius, 0.0);
     let mut temperature_map = Map::<f32>::new(heightmap.radius, 0.0);
     let mut atmo_humidity = Map::<f32>::new(heightmap.radius, 0.0); // relative humidity
-    let mut frontier = (-WORLD_RADIUS..=WORLD_RADIUS).map(|x| WindSlice {
+    let mut frontier = (-radius..=radius).map(|x| WindSlice {
         coord: start_pos + perpendicular_dir * x,
         humidity: 0.0,
         temperature: 0.0,
@@ -396,6 +399,8 @@ pub fn simulate_air(heightmap: &Map<f32>, sea: &HashSet<Coord>, scan_dir: CoordV
         // We treat it as a line advancing and gaining/losing water and temperature.
         // Water is lost when the partial pressure of water in the air is greater than the saturation vapour pressure.
         // Temperature changes with height based on a slightly dubious equation I derived.
+        // TODO: recalculate this; Wikipedia says temperature change is actually
+        // because of air pressure changes and not directly derived from GPE.
         for slice in frontier.iter_mut() {
             if heightmap.in_range(slice.coord) {
                 any_in_range = true;
@@ -408,7 +413,7 @@ pub fn simulate_air(heightmap: &Map<f32>, sea: &HashSet<Coord>, scan_dir: CoordV
                     slice.temperature = BASE_TEMPERATURE;
                     slice.humidity = max_water * 0.9;
                     slice.last_height = SEA_LEVEL;
-                    
+
                 } else {
                     let excess = (slice.humidity - max_water).max(0.0);
                     slice.humidity -= excess;
@@ -431,16 +436,26 @@ pub fn simulate_air(heightmap: &Map<f32>, sea: &HashSet<Coord>, scan_dir: CoordV
         let mut next_temperature = vec![0.0; frontier.len()];
         let mut next_humidity = vec![0.0; frontier.len()];
         // Smooth out temperature and humidity.
-        for i in 1..(frontier.len()-1) { 
-            next_temperature[i] = 1.0/3.0 * (frontier[i-1].temperature + frontier[i+1].temperature + frontier[i].temperature);
-            next_humidity[i] = 1.0/3.0 * (frontier[i-1].humidity + frontier[i+1].humidity + frontier[i].humidity);
-        }
+        // TODO at some point: replace with generalized convolution
+        let frontier_len = frontier.len();
+        next_temperature.par_iter_mut().enumerate().for_each(|(i, out)| {
+            if i == 0 || i + 1 >= frontier_len {
+                return;
+            }
+            *out = 1.0/3.0 * (frontier[i-1].temperature + frontier[i+1].temperature + frontier[i].temperature);
+        });
+        next_humidity.par_iter_mut().enumerate().for_each(|(i, out)| {
+            if i == 0 || i + 1 >= frontier_len {
+                return;
+            }
+            *out = 1.0/3.0 * (frontier[i-1].humidity + frontier[i+1].humidity + frontier[i].humidity);
+        });
 
-        for (i, slice) in frontier.iter_mut().enumerate() {
+        frontier.par_iter_mut().enumerate().for_each(|(i, slice)| {
             slice.temperature = next_temperature[i];
             slice.humidity = next_humidity[i];
-        }
-        
+        });
+
         if !any_in_range { break; }
     }
 
@@ -463,13 +478,13 @@ pub enum TerrainType {
 
 #[derive(Debug, Clone, Serialize, Deserialize)]
 pub struct GeneratedWorld {
-    heightmap: Map<f32>,
-    terrain: Map<TerrainType>,
-    groundwater: Map<f32>,
-    salt: Map<f32>,
-    atmo_humidity: Map<f32>,
-    temperature: Map<f32>,
-    soil_nutrients: Map<f32>,
+    pub heightmap: Map<f32>,
+    pub terrain: Map<TerrainType>,
+    pub groundwater: Map<f32>,
+    pub salt: Map<f32>,
+    pub atmo_humidity: Map<f32>,
+    pub temperature: Map<f32>,
+    pub soil_nutrients: Map<f32>,
     pub radius: i32
 }
 
@@ -514,11 +529,11 @@ pub fn generate_world() -> GeneratedWorld {
 }
 
 impl TerrainType {
-    pub fn entry_cost(&self) -> Option<i64> { 
+    pub fn entry_cost(&self) -> Option<i64> {
         match *self {
             Self::Empty => Some(0),
             Self::Wall => None,
-            Self::ShallowWater => Some(10), 
+            Self::ShallowWater => Some(10),
             Self::DeepWater => None,
             Self::Contour => Some(1)
         }
@@ -538,15 +553,30 @@ impl TerrainType {
 impl GeneratedWorld {
     pub fn get_terrain(&self, pos: Coord) -> TerrainType {
         let distance = hex_distance(pos, Coord::origin());
-    
+
         if distance >= self.heightmap.radius {
             return TerrainType::Wall
         }
-        
+
         self.terrain[pos].clone()
     }
 
     pub fn radius(&self) -> i32 {
         self.heightmap.radius
     }
-}
\ No newline at end of file
+}
+
+#[cfg(test)]
+mod test {
+    use super::{Map, smooth, WORLD_RADIUS};
+    use test::bench::Bencher;
+
+    #[bench]
+    fn bench_smooth_map(b: &mut Bencher) {
+        use std::hint::black_box;
+        b.iter(|| {
+            let map = black_box(Map::new(WORLD_RADIUS, 0.0f32));
+            smooth(&map, 3);
+        });
+    }
+}
diff --git a/src/worldgen_test.rs b/src/worldgen_test.rs
index c09f29d..36995c0 100644
--- a/src/worldgen_test.rs
+++ b/src/worldgen_test.rs
@@ -1,25 +1,204 @@
+#![feature(test)]
+extern crate test;
+
 use anyhow::Result;
+use argh::FromArgs;
 use image::{ImageBuffer, Rgb};
-use std::collections::{HashMap, HashSet};
+use std::collections::HashMap;
+use std::time::Instant;
 
 mod worldgen;
 mod map;
 
-use worldgen::*;
 use map::*;
+use worldgen::*;
+
+#[derive(FromArgs)]
+/// Render worldgen debug fields to a PNG.
+struct Args {
+    /// world radius
+    #[argh(option, default = "WORLD_RADIUS")]
+    radius: i32,
+    /// output path
+    #[argh(option, default = "String::from(\"./out.png\")")]
+    output: String,
+    /// first channel field (defaults to 0 if omitted)
+    #[argh(option)]
+    c1: Option<String>,
+    /// second channel field (defaults to 0 if omitted)
+    #[argh(option)]
+    c2: Option<String>,
+    /// third channel field (defaults to 0 if omitted)
+    #[argh(option)]
+    c3: Option<String>,
+    /// color space: rgb or oklab
+    #[argh(option, default = "String::from(\"oklab\")")]
+    color_space: String,
+    /// percentile-normalize each selected channel
+    #[argh(switch)]
+    normalize: bool,
+}
+
+#[derive(Clone, Copy)]
+enum Field {
+    Height,
+    Rain,
+    Water,
+    Groundwater,
+    Salt,
+    Temperature,
+    Humidity,
+    Soil,
+    Contour,
+    SeaDistance,
+}
+
+impl Field {
+    fn parse(s: &str) -> Result<Self> {
+        match s {
+            "height" => Ok(Self::Height),
+            "rain" => Ok(Self::Rain),
+            "water" => Ok(Self::Water),
+            "groundwater" => Ok(Self::Groundwater),
+            "salt" => Ok(Self::Salt),
+            "temperature" => Ok(Self::Temperature),
+            "humidity" => Ok(Self::Humidity),
+            "soil" => Ok(Self::Soil),
+            "contour" => Ok(Self::Contour),
+            "sea_distance" => Ok(Self::SeaDistance),
+            _ => anyhow::bail!("unknown field: {s}"),
+        }
+    }
+}
+
+struct RenderData<'a> {
+    heightmap: &'a Map<f32>,
+    rain: &'a Map<f32>,
+    water: &'a Map<f32>,
+    groundwater: &'a Map<f32>,
+    salt: &'a Map<f32>,
+    temperature: &'a Map<f32>,
+    humidity: &'a Map<f32>,
+    soil: &'a Map<f32>,
+    contour_points: &'a HashMap<Coord, u8>,
+    sea_distance: &'a Map<f32>,
+}
+
+fn sample_field(field: Option<Field>, position: Coord, data: &RenderData) -> f32 {
+    match field {
+        None => 0.0,
+        Some(field) => match field {
+        Field::Height => ((data.heightmap[position] + 1.0) * 0.5).clamp(0.0, 1.0),
+        Field::Rain => data.rain[position].clamp(0.0, 1.0),
+        Field::Water => data.water[position].min(1.0),
+        Field::Groundwater => data.groundwater[position].clamp(0.0, 1.0),
+        Field::Salt => data.salt[position].clamp(0.0, 1.0),
+        Field::Temperature => data.temperature[position].clamp(0.0, 1.0),
+        Field::Humidity => data.humidity[position].clamp(0.0, 1.0),
+        Field::Soil => data.soil[position].clamp(0.0, 1.0),
+        Field::Contour => data.contour_points.get(&position).copied().unwrap_or_default() as f32 / 255.0,
+        Field::SeaDistance => data.sea_distance[position].clamp(0.0, 1.0),
+        },
+    }
+}
+
+fn field_range(field: Option<Field>, data: &RenderData) -> (f32, f32) {
+    if field.is_none() {
+        return (0.0, 1.0);
+    }
+    let mut min = f32::INFINITY;
+    let mut max = f32::NEG_INFINITY;
+    for (position, _) in data.heightmap.iter() {
+        let v = sample_field(field, position, data);
+        min = min.min(v);
+        max = max.max(v);
+    }
+    if (max - min).abs() < f32::EPSILON {
+        (0.0, 1.0)
+    } else {
+        (min, max)
+    }
+}
+
+fn to_rgb(c1: f32, c2: f32, c3: f32, color_space: &str) -> [u8; 3] {
+    fn linear_to_srgb(x: f32) -> f32 {
+        if x <= 0.0031308 {
+            12.92 * x
+        } else {
+            1.055 * x.powf(1.0 / 2.4) - 0.055
+        }
+    }
+
+    let (r, g, b) = match color_space {
+        "rgb" => (c1, c2, c3),
+        "oklab" => {
+            let l = c1.clamp(0.0, 1.0);
+            let a = c2 * 2.0 - 1.0;
+            let b = c3 * 2.0 - 1.0;
+
+            let l_ = l + 0.3963377774 * a + 0.2158037573 * b;
+            let m_ = l - 0.1055613458 * a - 0.0638541728 * b;
+            let s_ = l - 0.0894841775 * a - 1.2914855480 * b;
+
+            let l3 = l_ * l_ * l_;
+            let m3 = m_ * m_ * m_;
+            let s3 = s_ * s_ * s_;
+
+            let r_lin = 4.0767416621 * l3 - 3.3077115913 * m3 + 0.2309699292 * s3;
+            let g_lin = -1.2684380046 * l3 + 2.6097574011 * m3 - 0.3413193965 * s3;
+            let b_lin = -0.0041960863 * l3 - 0.7034186147 * m3 + 1.7076147010 * s3;
+
+            (
+                linear_to_srgb(r_lin).clamp(0.0, 1.0),
+                linear_to_srgb(g_lin).clamp(0.0, 1.0),
+                linear_to_srgb(b_lin).clamp(0.0, 1.0),
+            )
+        }
+        _ => (c1, c2, c3),
+    };
+    [(r * 255.0) as u8, (g * 255.0) as u8, (b * 255.0) as u8]
+}
 
 fn main() -> Result<()> {
-    let mut heightmap = generate_heights();
-    let (sinks, sea) = get_sea(&heightmap);
+    let total_start = Instant::now();
+    let args: Args = argh::from_env();
+    let f1 = args.c1.as_deref().map(Field::parse).transpose()?;
+    let f2 = args.c2.as_deref().map(Field::parse).transpose()?;
+    let f3 = args.c3.as_deref().map(Field::parse).transpose()?;
+    let radius = args.radius.max(1);
 
-    println!("wind...");
-    let (rain, temperature, atmo_humidity) = simulate_air(&heightmap, &sea, CoordVec::new(0, -1), CoordVec::new(1, 0));
+    let t = Instant::now();
+    let mut heightmap = generate_heights_with_radius(radius);
+    println!("heights: {:.3}s", t.elapsed().as_secs_f32());
 
     println!("hydro...");
+    let hydro_start = Instant::now();
+    let t = Instant::now();
+    let (sinks, sea) = get_sea(&heightmap);
+    println!("  sea/sinks: {:.3}s", t.elapsed().as_secs_f32());
+    let t = Instant::now();
+    let (rain, temperature, atmo_humidity) =
+        simulate_air(&heightmap, &sea, CoordVec::new(0, -1), CoordVec::new(1, 0));
+    println!("  air: {:.3}s", t.elapsed().as_secs_f32());
+    let t = Instant::now();
     let (water, salt) = simulate_water(&mut heightmap, &rain, &sea, &sinks);
+    println!("  water: {:.3}s", t.elapsed().as_secs_f32());
+    let t = Instant::now();
+    let groundwater = compute_groundwater(&water, &rain, &heightmap);
+    println!("  groundwater: {:.3}s", t.elapsed().as_secs_f32());
+    let t = Instant::now();
+    let mut sea_distance = distance_map(heightmap.radius, sea.iter().copied());
+    let radius = heightmap.radius as f32;
+    for (_, value) in sea_distance.iter_mut() {
+        *value = (*value / radius).clamp(0.0, 1.0);
+    }
+    println!("  sea distance: {:.3}s", t.elapsed().as_secs_f32());
+    println!("hydro total: {:.3}s", hydro_start.elapsed().as_secs_f32());
 
     println!("contours...");
+    let t = Instant::now();
     let contours = generate_contours(&heightmap, 0.15);
+    println!("contours: {:.3}s", t.elapsed().as_secs_f32());
     let mut contour_points = HashMap::new();
 
     for (point, x1, x2, _) in contours {
@@ -28,26 +207,51 @@ fn main() -> Result<()> {
         *entry = std::cmp::max(*entry, (steepness * 4000.0).abs() as u8);
     }
 
-    println!("groundwater...");
-    let groundwater = compute_groundwater(&water, &rain, &heightmap);
-
     println!("soil...");
+    let t = Instant::now();
     let soil_nutrients = soil_nutrients(&groundwater);
+    println!("soil: {:.3}s", t.elapsed().as_secs_f32());
 
     println!("rendering...");
-    let mut image = ImageBuffer::from_pixel((WORLD_RADIUS * 2 + 1) as u32, (WORLD_RADIUS * 2 + 1) as u32, Rgb::from([0u8, 0, 0]));
+    let t = Instant::now();
+    let image_radius = heightmap.radius;
+    let mut image = ImageBuffer::from_pixel((image_radius * 2 + 1) as u32, (image_radius * 2 + 1) as u32, Rgb::from([0u8, 0, 0]));
+    let render_data = RenderData {
+        heightmap: &heightmap,
+        rain: &rain,
+        water: &water,
+        groundwater: &groundwater,
+        salt: &salt,
+        temperature: &temperature,
+        humidity: &atmo_humidity,
+        soil: &soil_nutrients,
+        contour_points: &contour_points,
+        sea_distance: &sea_distance,
+    };
+    let r1 = field_range(f1, &render_data);
+    let r2 = field_range(f2, &render_data);
+    let r3 = field_range(f3, &render_data);
 
-    for (position, height) in heightmap.iter() {
-        let col = position.x + (position.y - (position.y & 1)) / 2 + WORLD_RADIUS;
-        let row = position.y + WORLD_RADIUS;
-        let height = *height * 0.5 + 1.0;
-        let green_channel = height;
-        let red_channel = soil_nutrients[position];
-        let blue_channel = water[position].min(1.0);
-        image.put_pixel(col as u32, row as u32, Rgb::from([(red_channel * 255.0) as u8, (green_channel * 255.0) as u8, (blue_channel * 255.0) as u8]));
+    for (position, _) in heightmap.iter() {
+        let col = position.x + (position.y - (position.y & 1)) / 2 + image_radius;
+        let row = position.y + image_radius;
+        let mut c1 = sample_field(f1, position, &render_data);
+        let mut c2 = sample_field(f2, position, &render_data);
+        let mut c3 = sample_field(f3, position, &render_data);
+        if args.normalize {
+            c1 = ((c1 - r1.0) / (r1.1 - r1.0)).clamp(0.0, 1.0);
+            c2 = ((c2 - r2.0) / (r2.1 - r2.0)).clamp(0.0, 1.0);
+            c3 = ((c3 - r3.0) / (r3.1 - r3.0)).clamp(0.0, 1.0);
+        }
+        let rgb = to_rgb(c1, c2, c3, &args.color_space);
+        image.put_pixel(col as u32, row as u32, Rgb::from(rgb));
     }
+    println!("render: {:.3}s", t.elapsed().as_secs_f32());
 
-    image.save("./out.png")?;
+    let t = Instant::now();
+    image.save(args.output)?;
+    println!("save: {:.3}s", t.elapsed().as_secs_f32());
+    println!("total: {:.3}s", total_start.elapsed().as_secs_f32());
 
     Ok(())
-}
\ No newline at end of file
+}