% Simulation of Rands et al. state-based synchronisation model (Chapter 6, Box C)
% Written by David Sumpter

%Set up parameter values
b=10;
c=5;
e=4;
r=0.5;
inc=1;
inc2=0;

endtime=1000;

%Set up initial states
s=zeros(endtime,2)
forage=zeros(endtime,2)
s(1,1)=1;
s(1,2)=5;
forage(1,1)=0;
forage(1,2)=0;


for t=1:endtime
    
    %Decide new state for individual 1. 
    if s(t,1)<=b/c
        forage(t+1,1)=1;
    else
        if s(t,1)<=b/(c-e)
            %In this case wait to see what the other does before deciding
            forage(t+1,1)=2;
        else
            forage(t+1,1)=0;
        end
    end
    if s(t,2)<=b/c
        forage(t+1,2)=1;
        if forage(t+1,1)==2
            forage(t+1,1)=1;
        end
    else
        if s(t,2)<=b/(c-e)
            if forage(t+1,1)<2
                %Do what the other individual does.
                forage(t+1,2)=forage(t+1,1);
            else
                %Continue to do what you did before.
                forage(t+1,1)=forage(t,1);
                forage(t+1,2)=forage(t,2);
            end
        else
            forage(t+1,2)=0;
            if forage(t+1,1)==2
                forage(t+1,1)=0;
            end
        end
    end
    %Update state including random component.
    s(t+1,:)=s(t,:)+forage(t+1,:).*(inc+inc2*(s(t,:)>mean(s(t,:))))-(1-forage(t+1,:))*r+0.5*(rand(1,2)-2)/2;
end

%Plot simulation outcome
figure(1)
plot(s(:,1),'k:')
hold on
plot(s(:,2),'k')
axis([1 1000 0 13])
xlabel('Day (t)')
ylabel('Nutritional State (s)')
hold off

%Make phase diagram
figure(2)
plot([b/(c-e) b/(c-e)],[0 b/(c-e)],'k')
hold on
plot([b/c b/c],[b/c 13],'k')
plot([b/c b/(c-e)],[b/c b/c],'k')
plot([b/c b/(c-e)],[b/(c-e) b/(c-e)],'k')
axis([0 13 0 13])
hold off
text(1,1,'Forage')
text(11,11,'Rest')
text(4,6,'Forage if partner foragers and rest if partner rests')
xlabel('Focal individualīs nutritional state (s_1)')
ylabel('Partnerīs nutritional state (s_2)')

%gtext('b/c')
%gtext('b/c')
%gtext('b/(c-e)')
%gtext('b/(c-e)')