BeginPackage["CSM`Aggregate`"]

initialGrid::usage = "initialGrid[n] gives the initial configuration
	with a radius n."
occupy::usage = "occupy[grid, {x, y}] occupies the grid position {x, y}."
walk::usage = "walk[grid] aggregates one particle."
gridGraphics::usage = "gridGraphics[g] generates a graphics of grid g."
aggregate::usage = "aggregate[grid, n] performs n experiments (with walk[])."
orthogonal::usage = "orthogonal is the orthogonal neighborhood."
diagonal::usage = "diagonal is the diagonal neighborhood."
neighborhood::usage = "neighborhood is the current neighborhood."

Begin["`Private`"]

orthogonal = {{1,0}, {0,1}, {-1,0}, {0,-1}}
diagonal = Complement[Flatten[Outer[List, {-1,0,1}, {-1,0,1}], 1], {{0, 0}}]
neighborhood = orthogonal (* default *)
randomDirection :=
    neighborhood[[ Random[Integer, {1, Length[neighborhood]}] ]]

occupied = 2; neighboring = 1; empty = 0;

occupy[grid_, {x_, y_}] :=
    Module[{f = grid},
      Apply[ Function[{dx,dy}, f[[x+dx, y+dy]] = Max[f[[x+dx, y+dy]], neighboring]],
             neighborhood, {1} ];
      f[[x, y]] = occupied;
      f ]

walk[grid_] :=
    Module[{x, y, n = Length[grid], c = (Length[grid]-1)/2},
    While[ True,
        {x, y} = Table[Random[Integer, {2, n-1}], {2}];
        If[ grid[[x, y]] == occupied, Continue[] ]; (* already occupied *)
        While[1 < x < n && 1 < y < n && grid[[x, y]] == empty,
              {x, y} = {x, y} + randomDirection ];
        If[ 1 < x < n && 1 < y < n,
            Return[occupy[grid, {x, y}]] ]        (* ok *)
    ]]

initialGrid[n_] := occupy[ Table[0, {2n+1}, {2n+1}], {n+1, n+1} ]
aggregate[grid_, m_] := Nest[ walk, grid, m ]

coloring = {occupied -> 0, neighboring -> 1, empty -> 1};
gridGraphics[ grid_, opts___ ] := Graphics[Raster[grid /. coloring],
    opts, AspectRatio->1, Frame->True, FrameTicks->None]

End[]

Protect[ initialGrid, occupy, walk, gridGraphics, aggregate ]
Protect[ orthogonal, diagonal ]

EndPackage[]
