diff --git a/generate.lua b/generate.lua
new file mode 100644
index 0000000..ffbaf66
--- /dev/null
+++ b/generate.lua
@@ -0,0 +1,63 @@
+local EvolutionModel = dofile(mapgen_rivers.modpath .. '/erosion.lua')
+local twist = dofile(mapgen_rivers.modpath .. '/twist.lua')
+
+local size = {x=1000, y=1000}
+local blocksize = 12
+
+local np_base = {
+	offset = 0,
+	scale = 200,
+	seed = 2469,
+	octaves = 8,
+	spread = {x=4000/blocksize, y=4000/blocksize, z=5},
+	persist = 0.6,
+	lacunarity = 2,
+}
+
+local time = 10
+local time_step = 1
+local niter = math.ceil(time/time_step)
+time_step = time / niter
+
+local function generate()
+	local nobj_base = minetest.get_perlin_map(np_base, size)
+	nobj_base:calc_3d_map({x=0, y=0, z=0})
+
+	local dem = nobj_base:get_map_slice({z=1}, {z=1})
+	dem.X = size.x
+	dem.Y = size.y
+
+	local model = EvolutionModel()
+	model.dem = dem
+	local ref_dem = model:define_isostasy(dem)
+
+	for i=1, niter do
+		nobj_base:get_map_slice({z=i+1}, {z=1}, ref_dem)
+
+		model:diffuse(time_step)
+		model:flow()
+		model:erode(time_step)
+		model:isostasy()
+	end
+	model:flow()
+
+	local mfloor = math.floor
+	local mmin, mmax = math.min, math.max
+	local offset_x, offset_y = twist(model.dirs, model.rivers, 5)
+	for i=1, X*Y do
+		offset_x[i] = mmin(mmax(offset_x[i]*256, -128), 127)
+		offset_y[i] = mmin(mmax(offset_y[i]*256, -128), 127)
+	end
+
+	mapgen_rivers.write_map('dem', model.dem, 2)
+	mapgen_rivers.write_map('lakes', model.lakes, 2)
+	mapgen_rivers.write_map('dirs', model.dirs, 1)
+	mapgen_rivers.write_map('rivers', model.rivers, 4)
+	mapgen_rivers.write_map('offset_x', offset_x, 1)
+	mapgen_rivers.write_map('offset_y', offset_y, 1)
+	local sfile = io.open(mapgen_rivers.world_data_path .. 'size', "w")
+	sfile:write(X..'\n'..Y)
+	sfile:close()
+end
+
+return generate
diff --git a/load.lua b/load.lua
index 60b56d6..859c45f 100644
--- a/load.lua
+++ b/load.lua
@@ -1,6 +1,6 @@
 local worldpath = minetest.world_data_path
 
-local function load_map(filename, bytes, signed, size)
+function mapgen_rivers.load_map(filename, bytes, signed, size)
 	local file = io.open(worldpath .. filename, 'rb')
 	local data = file:read('*all')
 	if #data < bytes*size then
@@ -28,4 +28,28 @@ local function load_map(filename, bytes, signed, size)
 	return map
 end
 
-return load_map
+function mapgen_rivers.write_map(filename, data, bytes)
+    local size = #data
+    local file = io.open(worldpath .. filename, 'wb')
+    local mfloor = math.floor
+    local schar = string.char
+    local upack = unpack
+
+    local bytelist = {}
+    for j=1, bytes do
+        bytelist[j] = 0
+    end
+
+    for i=1, size do
+        local n = mfloor(data[i])
+        for j=bytes, 2, -1 do
+            bytelist[j] = n % 256
+            n = mfloor(n / 256)
+        end
+        bytelist[1] = n % 256
+
+        file:write(schar(upack(bytelist)))
+    end
+
+    file:close()
+end
diff --git a/polygons.lua b/polygons.lua
index 66b0b12..7ef9538 100644
--- a/polygons.lua
+++ b/polygons.lua
@@ -7,51 +7,39 @@ end
 local world_data_path = minetest.world_data_path
 minetest.mkdir(world_data_path)
 
-local load_map = dofile(modpath .. 'load.lua')
+dofile(modpath .. 'load.lua')
 
-local function copy_if_needed(filename)
-	local wfilename = world_data_path..filename
-	local wfile = io.open(wfilename, 'rb')
-	if wfile then
-		wfile:close()
-		return
-	end
-	local mfilename = mod_data_path..filename
-	local mfile = io.open(mfilename, 'rb')
-	local wfile = io.open(wfilename, 'wb')
-	wfile:write(mfile:read("*all"))
-	mfile:close()
-	wfile:close()
+-- Try to read file 'size'
+local sfile = io.open(world_data_path..'size', 'r')
+if not sfile then
+	-- Generate a map!!
+	local generate = dofile(mapgen_rivers.modpath .. '/generate.lua')
+	generate()
+	sfile = io.open(world_data_path..'size', 'r')
 end
 
-copy_if_needed('size')
-local sfile = io.open(world_data_path..'size', 'r')
+-- Read the map
+-- TODO: if data has just been generated, transfer it without reloading everything
 local X = tonumber(sfile:read('*l'))
 local Z = tonumber(sfile:read('*l'))
 sfile:close()
 
-copy_if_needed('dem')
-local dem = load_map('dem', 2, true, X*Z)
-copy_if_needed('lakes')
-local lakes = load_map('lakes', 2, true, X*Z)
-copy_if_needed('dirs')
-local dirs = load_map('dirs', 1, false, X*Z)
-copy_if_needed('rivers')
-local rivers = load_map('rivers', 4, false, X*Z)
+local dem = mapgen_rivers.load_map('dem', 2, true, X*Z)
+local lakes = mapgen_rivers.load_map('lakes', 2, true, X*Z)
+local dirs = mapen_rivers.load_map('dirs', 1, false, X*Z)
+local rivers = mapgen_rivers.load_map('rivers', 4, false, X*Z)
 
-copy_if_needed('offset_x')
 local offset_x = load_map('offset_x', 1, true, X*Z)
 for k, v in ipairs(offset_x) do
 	offset_x[k] = (v+0.5)/256
 end
 
-copy_if_needed('offset_y')
 local offset_z = load_map('offset_y', 1, true, X*Z)
 for k, v in ipairs(offset_z) do
 	offset_z[k] = (v+0.5)/256
 end
+-- Should have finished loading
 
--- To index a flat array representing a 2D map
 local function index(x, z)
 	return z*X+x+1
 end
diff --git a/terrainlib_lua/twist.lua b/terrainlib_lua/twist.lua
index e322c68..4674ed7 100644
--- a/terrainlib_lua/twist.lua
+++ b/terrainlib_lua/twist.lua
@@ -27,6 +27,7 @@ local function get_bounds(dirs, rivers)
 end
 
 local function twist(dirs, rivers, n)
+	n = n or 5
 	local X, Y = dirs.X, dirs.Y
 	local bounds_x, bounds_y = get_bounds(dirs, rivers)
 	local dn = 0.5 / n