% Simulation of bird feeding (Chapter 3, Box C)
% Written by David Sumpter

clear
close all

% Spontaneous arrival in absence of other chickens.
Aspon=0.01;  
% Maximum probability of arrival.
Amax=0.4;     
% Threshold response level.
AT=2.3;
% Steepness of response, must be greater than 1.
Ak=4;       

%Probability per individual of leaving
Pleave = 0.2;
        
%Plot the response function
figure(10)
a=[0:0.05:5];
f=(5-a).*(Aspon+(Amax-Aspon)*(a.^Ak./(AT^Ak +a.^Ak)))-Pleave*a;
plot(a,f)
ylabel('Probability of arrival in next 5 seconds');
xlabel('Average number at feeder')


%This is the simple simulation where there is only one site with unlimited capacity.
spaces=1;
simtime=10000;

%Number of chickens feeding.
c=zeros(simtime,1);

%Number of chickens not feeding.
out=zeros(simtime,1);
out(1)=7;

for t=2:simtime
    c(t)=c(t-1);
    out(t)=out(t-1);
    onnow=c(t);
        
    Parrive=(Aspon+(Amax-Aspon)*(onnow.^Ak./(AT^Ak +onnow.^Ak)));
    
    %Arrival?
    if out(t-1)>0
        r=sum(rand(out(t-1),1)<Parrive);
        c(t)=c(t)+r;
        out(t)=out(t)-r;                
    end
    
    %Departure?
    if c(t-1)>0
        r=sum(rand(c(t-1),1)<Pleave);
        out(t)=out(t)+r;
        c(t)=c(t)-r;
    end
 
end

mf=mean(c(1200:simtime));
allfeed=c(1200:simtime);
    
%Make and plot distribution
infd=[0:1:7];
h=hist(allfeed,infd);
o=length(allfeed);

figure(1)
bar(infd,h/o,'w')
hold off
xlabel('Number at feeder')
ylabel('Frequency')
axis([-0.5 7+0.5 0 0.8])
hlx=xlabel('Number of birds at the feeder')
hly=ylabel('Frequency')
set(hlx,'FontSize',14);
set(hly,'FontSize',14);

figure(2)
a=[-0.5:0.05:7.5];
f=(7-a).*(Aspon+(Amax-Aspon)*(a.^Ak./(AT^Ak +a.^Ak)))
g=Pleave*a;
plot(a,f,'k')
hold on
plot(a,g,'k')
axis([0 7 0 1.3])
hlx=xlabel('Number of birds at the feeder')
hly=ylabel('Joining and leaving rates')
set(hlx,'FontSize',14);
set(hly,'FontSize',14);

%Plot time series
figure(3)
subplot(1,1,1)
plot([0:1:length(c(:,1))-1],c(:,1),'k')
hlx=xlabel('Time (minutes)')
hly=ylabel('Number of birds at the feeder')
axis([1000 2000 0 7])
set(hlx,'FontSize',14);
set(hly,'FontSize',14);

save birdsimout