

 This page uses content from Wikipedia. The original article was at Abelian sandpile model. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)



Implement the Abelian sandpile model also known as Bak–Tang–Wiesenfeld model. Its history, mathematical definition and properties can be found under its wikipedia article.

The task requires the creation of a 2D grid of arbitrary size on which "piles of sand" can be placed. Any "pile" that has 4 or more sand particles on it collapses, resulting in four particles being subtracted from the pile and distributed among its neighbors.

It is recommended to display the output in some kind of image format, as terminal emulators are usually too small to display images larger than a few dozen characters tall. As an example of how to accomplish this, see the Bitmap/Write a PPM file task.

Examples up to 2^30, wow!
javascript running on web
Examples:

0 0 0 0 0    0 0 0 0 0
0 0 0 0 0    0 0 1 0 0
0 0 4 0 0 -> 0 1 0 1 0
0 0 0 0 0    0 0 1 0 0
0 0 0 0 0    0 0 0 0 0

0 0 0 0 0    0 0 0 0 0
0 0 0 0 0    0 0 1 0 0
0 0 6 0 0 -> 0 1 2 1 0
0 0 0 0 0    0 0 1 0 0
0 0 0 0 0    0 0 0 0 0

0  0 0  0  0    0 0 1 0 0
0  0 0  0  0    0 2 1 2 0
0  0 16 0  0 -> 1 1 0 1 1
0  0 0  0  0    0 2 1 2 0
0  0 0  0  0    0 0 1 0 0

-- demo\rosetta\Abelian_sandpile_model.exw
include pGUI.e
 
Ihandle dlg, canvas
cdCanvas cddbuffer
 
sequence board = {{0,0,0},
                  {0,0,0},
                  {0,0,0}}
 
procedure drop(integer y, x)
    sequence moves = {}
    while true do
        board[y,x] += 1
        if board[y,x]>=4 then
            board[y,x] -= 4
            moves &= {{y,x-1},{y,x+1},{y-1,x},{y+1,x}}
        end if
        -- extend board if rqd (maintain a border of zeroes)
        if x=1 then                             -- extend left
            for i=1 to length(board) do
                board[i] = prepend(board[i],0)
            end for
            for i=1 to length(moves) do
                moves[i][2] += 1
            end for
        elsif x=length(board[1]) then           -- extend right
            for i=1 to length(board) do
                board[i] = append(board[i],0)
            end for
        end if
        -- (copy the all-0 lines from the other end...)
        if y=1 then                             -- extend up
            board = prepend(board,board[$])
            for i=1 to length(moves) do
                moves[i][1] += 1
            end for
        elsif y=length(board) then              -- extend down
            board = append(board,board[1])
        end if
        if length(moves)=0 then exit end if
        {y,x} = moves[$]
        moves = moves[1..$-1]
    end while   
    IupUpdate(canvas)
end procedure
 
function timer_cb(Ihandle /*ih*/)
    integer y = floor(length(board)/2)+1,
            x = floor(length(board[1])/2)+1
    drop(y,x)
    return IUP_DEFAULT
end function
 
function redraw_cb(Ihandle ih, integer /*posx*/, integer /*posy*/)
    IupGLMakeCurrent(ih)
    cdCanvasActivate(cddbuffer)
    cdCanvasClear(cddbuffer)
    for y=1 to length(board) do
        for x=1 to length(board[1]) do 
            integer c = board[y][x]
            if c!=0 then
                integer colour = {CD_VIOLET,CD_RED,CD_BLUE}[c]
                cdCanvasPixel(cddbuffer, x, y, colour)
            end if
        end for
    end for
    cdCanvasFlush(cddbuffer)
    return IUP_DEFAULT
end function
 
function map_cb(Ihandle ih)
    IupGLMakeCurrent(ih)
    atom res = IupGetDouble(NULL, "SCREENDPI")/25.4
    cddbuffer = cdCreateCanvas(CD_GL, "300x100 %g", {res})
    cdCanvasSetBackground(cddbuffer, CD_PARCHMENT)
    return IUP_DEFAULT
end function
 
procedure main()
    IupOpen()
    canvas = IupGLCanvas("RASTERSIZE=300x100")
    IupSetCallbacks({canvas}, {"ACTION", Icallback("redraw_cb"),
                               "MAP_CB", Icallback("map_cb")})
    dlg = IupDialog(canvas,"TITLE=\"Abelian sandpile model\"")
    IupShow(dlg)
    Ihandle timer = IupTimer(Icallback("timer_cb"), 10)
    IupMainLoop()
    IupClose()
end procedure
 
main()
