diff --git a/init.lua b/init.lua
index 30cc6eb..3472e66 100644
--- a/init.lua
+++ b/init.lua
@@ -10,3 +10,9 @@ noordstar_caves = {}
 
 -- Load features to influence cave shapes
 load("shape")
+
+-- Start engine to generate caves
+load("engine")
+
+
+-- For show, the following code COULD in theory be run elsewhere
diff --git a/lua/engine.lua b/lua/engine.lua
new file mode 100644
index 0000000..59cea27
--- /dev/null
+++ b/lua/engine.lua
@@ -0,0 +1,267 @@
+
+-- Convert 3d relative coordinates to an index on a flat array
+local function from_3d_to_flat(dx, dy, dz, nx, ny)
+    return (nx * ny * dz) + (nx * dy) + dx + 1
+end
+
+-- Convert an index on a flat array to 3d relative coordinates
+local function from_flat_to_3d(i, nx, ny)
+    return {
+        dx = (i - 1) % nx,
+        dy = (i - 1) // nx % ny,
+        dz = (i - 1) // nx // ny
+    }
+end
+
+-- Iterate over a 3d area, both indexing by the index and the ansolute position
+local function iter_3d_area(minp, maxp, callback)
+    local nx = maxp.x - minp.x + 1
+    local ny = maxp.y - minp.y + 1
+    local nz = maxp.z - minp.z + 1
+
+    for i = 1, nx * ny * nz do
+        local dpos = from_flat_to_3d(i, nx, ny)
+        local pos = {
+            x = minp.x + dpos.x,
+            y = minp.y + dpos.y,
+            z = minp.z + dpos.z,
+        }
+
+        callback(i, pos)
+    end
+end
+
+-- Get an enhanced function from the def function that warns us when the
+-- function does not behave properly
+local function enhanced_func(def)
+    if type(def.name) ~= "string" then
+        error("Invalid nameless shape definition")
+    elseif type(def.func) ~= "function" then
+        error(
+            "Invalid shape definition misses an adjustment function"
+        )
+    else
+        return function(pos, n)
+            local out = def.func(pos, n)
+
+            if type(out) == "number" then
+                return out
+            elseif n == nil then
+                error(
+                    "Shape " .. def.name .. " function must return a number. The input `n` was nil. Perhaps your `noise_params` field is invalid?"
+                )
+            else
+                error("Shape " .. def.name .. " function must return a number.")
+            end
+        end
+    end
+end
+
+-- Get a flat array of cave shape noise values from a given cave shape def
+local function get_flat_from_shape_def(def, minp, maxp)
+    local f = enhanced_func(def)
+
+    local nx = maxp.x - minp.x + 1
+    local ny = maxp.y - minp.y + 1
+    local nz = maxp.z - minp.z + 1
+
+    local noise_flat_map = {}
+
+    -- If noise parameters have been defined, fill the table with noise
+    -- If not, all values remain nil
+    if def.noise_params ~= nil then
+        local p = PerlinNoiseMap(def.noise_params, { x = nx, y = ny, z = nz })
+        
+        if nz == 1 then
+            p:get_2d_map_flat(minp, noise_flat_map)
+        else
+            p:get_3d_map_flat(minp, noise_flat_map)
+        end
+    end
+    
+    iter_3d_area(minp, maxp, function(i, pos)
+        noise_flat_map[i] = f(pos, noise_flat_map[i])
+    end)
+
+    return noise_flat_map
+end
+
+local function get_flat_from_noise_params(minp, maxp, noise_params)
+    local nx = maxp.x - minp.x + 1
+    local ny = maxp.y - minp.y + 1
+    local nz = maxp.z - minp.z + 1
+
+    local buffer = {}
+
+    local p = PerlinNoiseMap(noise_params, { x = nx, y = ny, z = nz })
+
+    if nz == 1 then
+        p:get_2d_map_flat(minp, buffer)
+    else
+        p:get_3d_map_flat(minp, buffer)
+    end
+
+    return buffer
+end
+
+-- Based on the number of cave shapes, calculate how quickly connectivity is
+-- meant to change
+local function get_connectivity_noise_params(shape_size)
+    local factor = math.max(math.abs(shape_size) ^ 0.5, 1)
+
+    return {
+        offset = 50,
+        scale = 50,
+        spread = { x = factor * 250, y = factor * 100, z = factor * 250 },
+        seed = 297948,
+        octaves = 2,
+        persistence = 0.2,
+        lacunarity = 2.0,
+        flags = "eased"
+    }
+end
+
+-- Based on the number of cave shapes, calculate how quickly verticality is
+-- meant to change
+local function get_verticality_noise_params(shape_size)
+    local factor = math.max(math.abs(shape_size) ^ 0.5, 1)
+
+    return {
+        offset = 50,
+        scale = 50,
+        spread = { x = factor * 100, y = factor * 250, z = factor * 100 },
+        seed = 35644,
+        octaves = 2,
+        persistence = 0.2,
+        lacunarity = 2.0,
+        flags = "eased"
+    }
+end
+
+-- Get the distance of each cave shape to calculate their weight.
+local function shape_distance_weight(shape_size, dx, dy)
+    local factor = math.max(math.abs(shape_size) ^ 0.5, 1)
+    local max_distance = 100 / factor
+
+    if dx^2 + dy^2 > max_distance^2 then
+        return 0
+    else
+        return 1 / (dx^5 + dy^5)
+    end
+end
+
+-- Get a flat map containing all flat threshold values
+local function get_threshold_flat(minp, maxp)
+    local connectivity = get_flat_from_noise_params(
+        minp, maxp,
+        get_connectivity_noise_params(#noordstar_caves.registered_shapes)
+    )
+    local verticality = get_flat_from_noise_params(
+        minp, maxp,
+        get_verticality_noise_params(#noordstar_caves.registered_shapes)
+    )
+
+    local noise = {}
+
+    -- Get noise for all cave shapes
+    for key, def in pairs(noordstar_caves.registered_shapes) do
+        noise[key] = {
+            key = key,
+            noise = get_flat_from_shape_def(def, minp, maxp),
+            def = def
+        }
+    end
+
+    local thresholds = {}
+
+    -- Fill the table
+    iter_3d_area(minp, maxp, function(i, pos)
+        local total = 0
+        local count = 0
+
+        local x = connectivity[i]
+        local y = verticality[i]
+
+        for _, n in pairs(noise) do
+            local v = n.noise[i]
+
+            local dx = math.abs(x - n.def.connectivity)
+            local dy = math.abs(y - n.def.verticality)
+
+            local w = math.abs(shape_distance_weight(
+                #noordstar_caves.registered_shapes, dx, dy
+            ))
+
+            total = total + v * w
+            count = count + w
+        end
+
+        if count <= 0 then
+            thresholds[i] = -1000
+        else
+            thresholds[i] = total / count
+        end
+    end)
+
+    return thresholds
+end
+
+
+-- Old minimum height for the overworld
+local old_overworld_min = mcl_vars.mg_overworld_min
+
+-- If another mod doesn't override the maximum world depth, we will assume that
+-- the world depth is the following value.
+local world_depth = -800
+
+-- Otherwise, this variable can be changed using the following function,
+-- which will also update all the other necessary variables
+function noordstar_caves.set_world_depth(h)
+    -- Set world depth variable
+    world_depth = h
+end
+noordstar_caves.set_world_depth(world_depth)
+
+-- The `cave_vastness` is a variable that determines the size of caves.
+-- The function takes in an x, y and z variable, and returns a number between
+-- 0 and 1.
+-- Low numbers mean very rare and tiny caves, while high numbers mean massive
+-- caves, with massive open spaces.
+-- If you wish to overwrite this function, it is good to keep in mind to:
+--  - Make sure that the output changes VERY SLOWLY over time
+--  - This function will be run a LOT so it is very performance sensitive
+noordstar_caves.cave_vastness = function(pos)
+    return math.abs(pos.y) / world_depth
+end
+
+-- Secretly, we're using an internal function that also adds a safe layer for
+-- when we're approaching bedrock levels
+local function cave_vastness(pos)
+    if world_depth + 20 < pos.y then
+        return noordstar_caves.cave_vastness(pos)
+    elseif world_depth + 5 < pos.y then
+        return noordstar_caves.cave_vastness(pos) * math.abs(y - world_depth - 5) / 15
+    else
+        return -1000
+    end
+end
+
+minetest.register_on_generated(function(minp, maxp, blockseed)
+    -- Get voxelmanip
+    local vm = minetest.get_mapgen_object("voxelmanip")
+    local data = vm:get_data()
+    
+    -- Get threshold values
+    local thresholds = get_threshold_flat(minp, maxp)
+    local air = minetest.get_content_id("air")
+
+    iter_3d_area(minp, maxp, function(i, pos)
+        if thresholds[i] >= cave_vastness(pos) then
+            data[i] = air
+        end
+    end)
+
+    -- Write all changes to the Minetest world
+    vm:set_data(data)
+    vm:write_to_map()
+end)