function init(w, h)
   width = w
   height = h

   R = width / 3
   R2 = R^2

   meshsize = 15 -- 2^n - 1 for nice look under antialiasing
   meshstepx = w / (meshsize + 1)
   meshstepy = h / (meshsize + 1)

   centerx = w / 2
   centery = h / 2

   range = math.min(width, height) / 2 * 0.7
end

-- C function reference:
-- setpixel(x, y, i) sets a grayscale pixel with intensity i
-- setpixel(x, y, i, r, g, b) sets an RGB pixel
-- drawline(x0, y0, x1, y1, ...) draws a line
-- cleartexture() clears the entire texture

function meshloc(x, y)
   return y * width + x
end

function makemesh()
   mesh = {}

   for y = 1, meshsize, 1 do
      for x = 1, meshsize, 1 do
         mesh[meshloc(x, y)] = {x * meshstepx, y * meshstepy, 100}
      end
   end

   return mesh
end

function displace(x, y, mesh)
   for i = 1, meshsize, 1 do
      for j = 1, meshsize, 1 do
         local dx = j * meshstepx - x
         local dy = i * meshstepy - y
         local r2 = dx^2 + dy^2
         if r2 <= R2 then
            local node = mesh[meshloc(j, i)]
            local r = math.sqrt(r2)
            local dpl = (R * r - r2) / 330

            node[1] = node[1] + (dx / r) * dpl
            node[2] = node[2] + (dy / r) * dpl
            node[3] = node[3] + 150 * (R2 - r2) / R2
         end
      end
   end
end

function plotmesh(mesh)
   -- plot grid
   for y = 1, meshsize - 1, 1 do
      for x = 1, meshsize - 1, 1 do
         local v1 = mesh[meshloc(x, y)]
         local v2 = mesh[meshloc(x + 1, y)]
         drawline(v1[1], v1[2], v2[1], v2[2], (v1[3] + v2[3]) / 2)
         
         local v2 = mesh[meshloc(x, y + 1)]
         drawline(v1[1], v1[2], v2[1], v2[2], (v1[3] + v2[3]) / 2)
      end
   end

   -- plot last edges
   for i = 1, meshsize - 1, 1 do
      local v1 = mesh[meshloc(meshsize, i)]
      local v2 = mesh[meshloc(meshsize, i + 1)]
      drawline(v1[1], v1[2], v2[1], v2[2], (v1[3] + v2[3]) / 2)

      local v1 = mesh[meshloc(i, meshsize)]
      local v2 = mesh[meshloc(i + 1, meshsize)]
      drawline(v1[1], v1[2], v2[1], v2[2], (v1[3] + v2[3]) / 2)
   end
end

function render(time_in_millis)
   local phase = 2 * math.pi * time_in_millis / 4000
   cleartexture()

   local mesh = makemesh()
   
   local ball1x = centerx + math.sin(phase * 0.5) * range
   local ball1y = centery - math.cos(phase * 1) * range
   displace(ball1x, ball1y, mesh)

   local ball2x = centerx + math.sin(phase * 1.2) * range
   local ball2y = centery - math.cos(phase * 0.6) * range
   displace(ball2x, ball2y, mesh)

   plotmesh(mesh)
end

