% Simulation of Niwa's model (Chapter 5, Box A)
% Written by David Sumpter

% Number of sites 
s=200;
k=1;

%Initialise sites, one individual per site.
x=zeros(1,s);
x(1,1:s)=1;

%Split probability.
p=0.05;

%Number of time steps
maxt=1000;
%Keeps track of how many individuals are at random site.
track=zeros(maxt,1);
%Keeps track of all sites.
store=zeros(maxt,s);


for t=1:maxt
    
    store(t,:)=x;
    y=zeros(1,s);
    for site=1:s

        r=rand;

    
        if ((x(site)>0) & (r<p))
            %Split has occured
            n1=floor(rand*x(site))+1;
            %Some remain at the site
            y(site)=x(site)+y(site)-n1;  
            site1=site;
            %The others move to a random neighbouring site.
            while (site1==site)
                site1=mod(site-1+floor(rand*3),s);
                if (site1==0)
                    site1=s;
                end
            end
            y(site1)=y(site1)+n1;
        else
            %No split just local movement.
            site1=site;
            while (site1==site)
                %Move to a random neighbouring site.
                site1=mod(site-1+floor(rand*3),s);
                if (site1==0)
                    site1=s;
                end
            end
            y(site1)=y(site1)+x(site);
        end
    end 
    
    x=y;
    track(t)=x(floor(rand*s)+1);
    
end

rangeN=[0:0.5:4.5]
range=[1:1:4000]
rangeN=range;

h=histc(track(1:maxt),range)
NPe=sum(h.*range')/sum(h)

NP=NPe*3.08;
fitw=exp(-rangeN/NP.*(1-exp(-rangeN/NP)/2))./rangeN
fitw=fitw/sum(fitw);

%Plot frequency distribution
figure(2)
subplot(1,1,1)
hpl=plot(log(range),log(h/sum(h)),'kx')
set(hpl,'MarkerSize',10);
hold on
hpl=plot(log(range),log(fitw),'k')
set(hpl,'MarkerSize',10);
hold off
axis([0 4.5 -7 0])
hlx=xlabel('log(Group size)');
hly=ylabel('log(Frequency)');
set(hlx,'FontSize',12);
set(hly,'FontSize',12);

%Plot time series through space
figure(1)
subplot(1,1,1)
imagesc(store,[0 10])
axis([0 200 0 1000]) 
colormap([1 [0.9:-1/63:0] ; 1 [0.9:-1/63:0] ; 1 [0.9:-1/63:0]]')
hlx=ylabel('Time')
hly=xlabel('Position')
set(hlx,'FontSize',12);
set(hly,'FontSize',12);