duration = 10;
length = 200;

northdensity = 0;
southdensity = 0;
eastdensity = 0;
westdensity = 0;

east = 1;
north = 2;
west = 4;
south = 8;

northfield = rand(length, length)< northdensity;
southfield = rand(length, length)< southdensity;
eastfield = rand(length, length)< eastdensity;
westfield = rand(length, length)< westdensity;

field = northfield * north + southfield * south + westfield * west + eastfield * east;


lookuptable = [ 0, 1,2,3,4,10,6,7,8,9,5,11,12,13,14,15 ];

colormap(gray(5));

smallsquare = length / 10
field(length/2 - smallsquare:length/2 + smallsquare,length/2 - smallsquare:length/2 + smallsquare)=15;


for step=1:duration
    # Move all of the particles
    eastmovement = bitand(shift(field ,1, 2), east);
    northmovement = bitand(shift(field, -1), north);
    westmovement = bitand(shift(field, -1, 2), west)';
    southmovement = bitand(shift(field, 1), south);

    density = eastmovement + southmovement + westmovement + northmovement;
    field  = lookuptable(density+1);

    image((1+bitand(density, north)+bitand(density, west)+bitand(density, east)+bitand(density, south))>2*5);title(num2str(step));drawnow;
#    pause;
end