Updated Ising model for critical temperature

-- Gotta modify the mod function because taking the mod of a negative number just returns that same number.
local function mod(x, y)
    return (x>=0) and math.fmod(x,y) or math.fmod(x,y)+y
end

local function initialize(width, height)
    local grid={}
    for i=0,width-1 do
        grid[i]={}
        for j=0,height-1 do
            --grid[i][j]=math.floor(2*math.random()-1)      -- Selecting random numbers in this way doesn't work and I have no idea why.
            grid[i][j]=2*math.random(2)-3
            --grid[i][j]=-1
        end
    end
    return grid
end

-- Copy a table.
local function tablecopy(t)
    local tcopy={}
    for k, v in pairs(t) do
        if type(v)=="table" then
            tcopy[k]=tablecopy(v)
        else
            tcopy[k]=v
        end
    end
    return tcopy
end

local function pickrandom(dimension)        -- dimension can be either width or height
    return math.random(dimension)-1         -- Subtract 1 because indexing starts at 0.
end

local function magnetization(grid)
    local mag=0
    local width, height = #grid+1, #grid[1]+1
    for i=0,#grid do
        for j=0,#grid[i] do
            mag=mag+grid[i][j]
        end
    end

    --mag=mag/(width*height)
    return mag
end

local function displaygrid(grid)
    local colorkey={"white", [-1]="black"}
    for i=0,#grid do
        for j=0,#grid[i] do
            gui.drawPixel(i, j, colorkey[grid[i][j]])
        end
    end
end

local function neighbors(grid, i, j, width, height)
    local l=mod(i-1, width)
    local r=mod(i+1, width)
    local u=mod(j-1, height)
    local d=mod(j+1, height)

    local nearest = {grid[l][j], grid[r][j], grid[i][u], grid[i][d]}

    return nearest
end

local function deltaU(grid, i, j, width, height)
    nearest=neighbors(grid, i, j, width, height)
    local sum=0
    for ii=1,#nearest do        -- Already using i as an argument, although Lua's scoping means we could use i if we wanted.
        sum=sum+nearest[ii]
    end

    return 2*grid[i][j]*sum
end

local function atrandom(grid, width, height, temp, iterations)
    for k=1,iterations do
        i, j = pickrandom(width), pickrandom(height)
        local Ediff = deltaU(grid, i, j, width, height)
        if Ediff<=0 or math.random()<math.exp(-Ediff/temp) then
            grid[i][j] = -grid[i][j]
        end
    end
    return grid
end

local width=240     -- Width and height are 240 and 160 (for GBA), but indexing starts at 0 so we subtract 1 when initializing the grid.  We still use these values for modular arithmetic.
local height=160
-- Critical temperature is ~2.2692
local kT=3
local rate=0.2*width*height     -- Number of iterations to perform before drawing all the pixels, which is computationally expensive.
local startframe=emu.framecount()
local grid = initialize(width, height)

while true do
    local mag = magnetization(grid)

    local t = emu.framecount()-startframe
    local A = 0.8
    local omega = math.pi/500
    local lambda = math.log(1/0.8)/500      -- omega*t=pi ==> exp(-lambda*t) = 0.8 ==> -lambda*pi/omega = ln(0.8) ==> lambda = omega/pi * ln(1/0.8)
    local temp = 2.2692 - A * math.exp(-lambda*t) * math.cos(omega*t)

    displaygrid(grid)
    local grid = atrandom(grid, width, height, temp, rate)
    --gui.pixelText(80, 20, kT,nil,"black")
    --gui.pixelText(80, 30, emu.framecount()-startframe,nil,"black")

    emu.frameadvance()
end