0-wave-render.lua

require('protectenv')()
local ffi = require('ffi')
local complex = require('complex')
local png = require('png')

local wlen = 4
local wnum = 2*math.pi / wlen

local map = ffi.new('complex_t[512][512]')

local function clamp(x)
    return 0.5 * (math.abs(x) - math.abs(x-1)) + 0.5
end

local function hankel01(z)
    local a = math.pi / 4
    local b = math.pi / 4
    local r = b / math.sqrt(z)
    return complex.frompolar(r, z - a)
end

local function hankel11(z)
    local a = 3 * math.pi / 4
    local b = math.pi / 4
    local r = b / math.sqrt(z)
    return complex.frompolar(r, z - a)
end

local function planesource(l, cx, cy, nx, ny)
    l = complex(l)
    return
        function(x, y)
            local dx = x - cx
            local dy = y - cy
            local r = dx * nx + dy * ny
            return l * complex.frompolar(1, wnum * r)
        end
end

local function pointsource(l, cx, cy)
    l = complex(l)
    return
        function(x, y)
            local dx = x - cx
            local dy = y - cy
            local r = math.sqrt(dx^2 + dy^2)
            if r < 1e-20 then
                r = 1e-20
            end
            return l * hankel01(wnum*r)
        end
end

local function dirsource(l, cx, cy, nx, ny)
    l = complex(l)
    return
        function(x, y)
            local dx = x - cx
            local dy = y - cy
            local r = math.sqrt(dx^2 + dy^2)
            if r < 1e-20 then
                r = 1e-20
            end
            local cosa = (dx*nx + dy*ny) / r
            -- if cosa < 0 then
                -- cosa = 0
            -- end
            return l * cosa * hankel11(wnum*r)
        end
end

local function sourcepair(sa, sb)
    return
        function(x, y)
            return sa(x, y) + sb(x, y)
        end
end

local sources = {}

local mainsource =
    sourcepair(
        pointsource(complex(1), 120.5, 240.5),
        pointsource(complex(1), 160.5, 305.5))
-- local angle = - 0.4 * math.pi
-- local mainsource = planesource(complex(0.07), pcx, pcy, math.cos(angle), math.sin(angle))

local function ampat(x, y)
    local z = complex()
    for i, source in ipairs(sources) do
        z = z + source(x, y)
    end
    return z
end

for i = -150, 0 do
    local cx = 200.5
    local cy = 255.5 + i
    local zm = mainsource(cx, cy)
    local zd = zm * complex(0, 0.5 * wnum)
    sources[#sources + 1] = dirsource(- 0.5 * zd, cx, cy, 1, 0)
    local dzmdx = 1000*(mainsource(cx+0.001, cy) - zm)
    local zp = dzmdx * complex(0, -0.5)
    sources[#sources + 1] = pointsource(- 0.5 * zp, cx, cy)
end

sources[#sources + 1] = mainsource

for y = 0, 511 do
    for x = 0, 511 do
        map[y][x] = ampat(x, y)
    end
end

local aimg = ffi.new('uint8_t[512][512][4]')
local bimg = ffi.new('uint8_t[512][512][4]')

local scale = 5

for y = 0, 511 do
    for x = 0, 511 do
        local z = map[y][x]
        local r = clamp(scale * (z.re))
        local g = clamp(scale * (-0.5 * z.re + math.sqrt(3/4) * z.im))
        local b = clamp(scale * (-0.5 * z.re - math.sqrt(3/4) * z.im))
        aimg[y][x][0] = math.sqrt(r) * 255
        aimg[y][x][1] = math.sqrt(g) * 255
        aimg[y][x][2] = math.sqrt(b) * 255
        aimg[y][x][3] = 255
        local lum = clamp((scale * z):magsqr())
        bimg[y][x][0] = math.sqrt(lum) * 255
        bimg[y][x][1] = math.sqrt(lum) * 255
        bimg[y][x][2] = math.sqrt(lum) * 255
        bimg[y][x][3] = 255
    end
end

png.writefile('a.png', 512, 512, aimg)
png.writefile('b.png', 512, 512, bimg)

1-png.lua

local png = {}

local ffi = require('ffi')

local libpng = ffi.load('libpng16.dll')

ffi.cdef[[
typedef struct png_control *png_controlp;
typedef struct
{
    png_controlp opaque;    /* Initialize to NULL, free with png_image_free */
    uint32_t  version;   /* Set to PNG_IMAGE_VERSION */
    uint32_t  width;     /* Image width in pixels (columns) */
    uint32_t  height;    /* Image height in pixels (rows) */
    uint32_t  format;    /* Image format as defined below */
    uint32_t  flags;     /* A bit mask containing informational flags */
    uint32_t  colormap_entries;
                            /* Number of entries in the color-map */
    uint32_t  warning_or_error;

    char       message[64];
} png_image, *png_imagep;

typedef struct png_color_struct
{
    uint8_t* red;
    uint8_t* green;
    uint8_t* blue;
} png_color;
typedef png_color * png_colorp;
typedef const png_color * png_const_colorp;

int png_image_begin_read_from_file(
    png_imagep image, const char *file_name);

int png_image_finish_read(
    png_imagep image,
    png_const_colorp background,
    void *buffer, int32_t row_stride,
    void *colormap);

void png_image_free(png_imagep image);

int png_image_write_to_file(
    png_imagep image,
    const char *file, int convert_to_8bit, const void *buffer,
    int32_t row_stride, const void *colormap);
]]

local PNG_IMAGE_VERSION = 1

local PNG_IMAGE_WARNING = 1
local PNG_IMAGE_ERROR = 2

local PNG_FORMAT_RGBA = 3

local function newimage()
    local i = ffi.gc(ffi.new('png_image'), libpng.png_image_free)
    i.version = PNG_IMAGE_VERSION
    return i
end

local function checkerror(i)
    local woe = bit.band(i.warning_or_error, 3) ~= 0
    if woe == 2 or woe == 3 then
        error(ffi.string(i.message))
    elseif woe == 1 then
        printf('png warning: ' .. ffi.string(i.message))
    end
end

function png.readfile(path)
    local i = newimage()
    libpng.png_image_begin_read_from_file(i, path)
    checkerror(i)
    local bm = ffi.new(ffi.typeof('uint8_t[$][$][4]', i.height, i.width))
    i.format = PNG_FORMAT_RGBA
    i.flags = 0
    libpng.png_image_finish_read(i, nil, bm, i.width * 4, nil)
    checkerror(i)
    libpng.png_image_free(i)
    return i.width, i.height, bm
end

function png.writefile(path, width, height, bm)
    local i = newimage()
    i.width = width
    i.height = height
    i.format = PNG_FORMAT_RGBA
    i.flags = 0
    libpng.png_image_write_to_file(i, path, 0, bm, width * 4, nil)
    checkerror(i)
end

return png

2-complex.lua

local ffi = require('ffi')

ffi.cdef[[
typedef struct complex_t {
    float re, im;
} complex_t;
]]

local new

local mtable = {}

mtable.__index = {}

function mtable.__index:magsqr()
    return self.re^2 + self.im^2
end

function mtable.__index:mag()
    return math.sqrt(self:magsqr())
end

function mtable.__index:arg()
    return math.atan2(self.im, self.re)
end

function mtable.__index:topolar()
    return self:mag(), self:arg()
end

function mtable.__index:conj()
    return new(self.re, -self.im)
end

function mtable:__tostring()
    return 'complex(' .. self.re .. ',' .. self.im .. ')'
end

function mtable.__unm(z1)
    return new(-z1.re, -z1.im)
end

function mtable.__add(z1, z2)
    if type(z1) == 'number' then
        return new(z1 + z2.re, z2.im)
    elseif type(z2) == 'number' then
        return new(z1.re + z2, z1.im)
    else
        return new(z1.re + z2.re, z1.im + z2.im)
    end
end

function mtable.__sub(z1, z2)
    if type(z1) == 'number' then
        return new(z1 - z2.re, -z2.im)
    elseif type(z2) == 'number' then
        return new(z1.re - z2, z1.im)
    else
        return new(z1.re - z2.re, z1.im - z2.im)
    end
end

function mtable.__mul(z1, z2)
    if type(z1) == 'number' then
        return new(z1*z2.re, z1*z2.im)
    elseif type(z2) == 'number' then
        return new(z1.re*z2, z1.im*z2)
    else
        return new(z1.re*z2.re - z1.im*z2.im, z1.im*z2.re + z1.re*z2.im)
    end
end

function mtable.__div(z1, z2)
    --[[
        (z1.re + z1.im I) / (z2.re + z2.im I)
        ((z1.re*z2.re + z1.im*z2.im) + (z1.im*z2.re - z1.re*z2.im) I) / (z2.re^2 + z2.im^2)

        z1.im == 0
        (z1.re*z2.re - z1.re*z2.im I) / (z2.re^2 + z2.im^2)

        z2.im == 0
        (z1.re + z1.im I) / z2.re
    ]]
    if type(z1) == 'number' then
        local d = z2:magsqr()
        local nre = z1*z2.re
        local nim = - z1*z2.im
        return new(nre/d, nim/d)
    elseif type(z2) == 'number' then
        return new(z1.re/z2, z1.im/z2)
    else
        local d = z2:magsqr()
        local nre = z1.re*z2.re + z1.im*z2.im
        local nim = z1.im*z2.re - z1.re*z2.im
        return new(nre/d, nim/d)
    end
end

new = ffi.metatype('complex_t', mtable)

local complex = setmetatable({},
    {__call = function(_, ...) return new(...) end}
)

function complex.frompolar(m, p)
    return new(m * math.cos(p), m * math.sin(p))
end

return complex

3-protectenv.lua

return function()
  local oldenv = getfenv(2)
  local newenv = setmetatable({}, {
    __index = oldenv,
    __newindex = function(t, k, v) error('attempt to set global: ' .. tostring(k)) end})
  setfenv(2, newenv)
end