%
%  math495/stat490 -- 03 nov 03 -- djm
%
%  w10stones.m:  

%  clear workspace & plots
clear;  %close all

%  experimental parameters
Nsize = 2*16;  Nstates=2;  Nsteps = 50;  

%  output parameters
Ncycle = 1;  Nfigs = 0;

%  initialize random number generator
phone = 2914814;
rand('state',phone);  %  initialize with my office phone number

%  initial data
ichoice = 2;
switch ichoice
    case{1}
        half = floor(Nsize/2);
        stones = zeros(Nsize,Nsize);
        stones(1:half,1:half) = ones(half,half);
    case{2}
        stones = floor(Nstates*rand(Nsize,Nsize)) + 1;
end

figure(1);
j = 1:Nsize;
image(j,j,stones,'CDataMapping','scaled');  hold on
title(['\bf stepping stones']);
xlabel(['\bf step # 0'])
axis equal;  axis image
pause(0)
print -djpeg90 -r60 file0.jpg

%  Nsteps generations
for kk = 1:Nsteps
    %  choose a random element (indexed on 0:Nsize-1 for mod to work)
    %  note transpose indexing for "image" command!!!
    ran_st = floor(Nsize*rand(2,1));
    
    %  neighbour coordinates - with periodicity!
    neigh = [mod(ran_st(1)+1,Nsize), mod(ran_st(2)+1,Nsize) ;
             mod(ran_st(1)+1,Nsize), mod(ran_st(2)+0,Nsize) ;
             mod(ran_st(1)+1,Nsize), mod(ran_st(2)-1,Nsize) ;
             mod(ran_st(1)+0,Nsize), mod(ran_st(2)-1,Nsize) ;
             mod(ran_st(1)-1,Nsize), mod(ran_st(2)-1,Nsize) ;
             mod(ran_st(1)-1,Nsize), mod(ran_st(2)+0,Nsize) ;
             mod(ran_st(1)-1,Nsize), mod(ran_st(2)+1,Nsize) ;
             mod(ran_st(1)+0,Nsize), mod(ran_st(2)+1,Nsize) ];
    centre = stones(ran_st(2)+1,ran_st(1)+1);

    if(mod(kk,Ncycle)==0)
        image(j,j,stones,'CDataMapping','scaled');  caxis([1 Nstates]);  hold on
        plot(ran_st(1)+1,ran_st(2)+1,'w*')
        plot(neigh(:,1)+1,neigh(:,2)+1,'wx')
        title(['\bf stepping stones']);
        xlabel(['\bf step # ',num2str(kk)])
        axis equal;  axis image
        hold off
        pause(1)
    end
    
    %  set the 8 neighbours
    stones(neigh(:,2)+1,neigh(:,1)+1) = centre;
    
    if(mod(kk,Ncycle)==0)
        caxis([1 Nstates])
        image(j,j,stones,'CDataMapping','scaled');  caxis([1 Nstates]);  hold on
        title(['\bf stepping stones']);
        xlabel(['\bf step # ',num2str(kk)])
        axis equal;  axis image
        
        if(Nfigs==1)
            filename = ['file' num2str(kk) '.jpg'];
            print(gcf,'-djpeg90','-r60',filename)
        end
        hold off
        pause(0)
    end
end
xlabel(['\bf step # ',num2str(kk)])