% Simulation of ant foraging model (Chapter 3, Box A)
% Written by David Sumpter

%Parameters, alpha and k.
%To simulate the bees change alpha=1.
alpha=2;
k=2;


%Make an example 
figure(1)
phi=6.3;
phi2=6.3*0.95;
alpha=2
k=2;
initpot1=0.3;
initpot2=0.8;
time=[0:0.01:40];
[T Y] =  ode45('FN',time,[initpot1,initpot2],[],phi,phi2,alpha,k);
plot(T,Y(:,1),'k')
hold on
plot(T,Y(:,2),'k:')
hlx=xlabel('Time (minutes)')
hly=ylabel('Number of ants crossing bridge per minute')
set(hlx,'FontSize',14);
set(hly,'FontSize',14);


%Range of flow (phi) values
vals=[0.01:0.02:10];

ci=0;
for phi=vals
    
    %Difference between phi1 and phi2. This is equivalent to qx=1 and qy=0.95
    phi2=phi*0.95;
    
    %Range of initial values
    invals=[0.1:0.2:phi-0.01];

    ci=ci+1;
    phi
    cj=0;
    for initpot1=invals;
        initpot2=phi-initpot1;
        cj=cj+1;
   
        time=[0:0.5:1000];

        %Calls the differential equation simulator. You need the file FN.
        [T Y] =  ode45('FN',time,[initpot1,initpot2],[],phi,phi2,alpha,k);

        finalY1(ci,cj)=Y(length(Y),1);
        finalY2(ci,cj)=Y(length(Y),2);
        
        startY1(ci,cj)=initpot1;
        param(ci,cj)=phi;
        directionY1(ci,cj)=finalY1(ci,cj)-startY1(ci,cj);

    end
end

%Make the bifurcation diagram
figure(2)
hold off
finalY1r=round(finalY1*10)/10;
      
for ci=1:1:size(param,1)
    shift=1;
    
    uf = unique(finalY1r(ci,:))
    for u=uf
        if ((u>0.001) & (sum(finalY1r(ci,:)==u)>1))
            for cj=1:size(param,2)          
                if ((finalY1r(ci,cj)==u) & (finalY1(ci,cj)~=0))
                    plot(vals(ci),finalY1(ci,cj),'k.')
                    hold on
                    if (mod(ci,40)==0)
                        
                        if ((finalY1r(ci,cj)==0) & (finalY2(ci,cj)>0.01))
                            %plot(param(ci,cj)+shift,startY1(ci,cj),'b.')
                            plot([param(ci,cj)+shift param(ci,cj)+shift],[startY1(ci,cj) finalY1(ci,cj)],'b')
                            plot(param(ci,cj)+shift,finalY1(ci,cj)+0.18,'bv')
                        else 
                            if startY1(ci,cj)>finalY1r(ci,cj)
                                %plot(param(ci,cj),startY1(ci,cj),'r.')
                                plot([param(ci,cj) param(ci,cj)],[startY1(ci,cj) finalY1(ci,cj)],'k')
                                plot(param(ci,cj),finalY1(ci,cj)+0.18,'kv')
                            else
                                %plot(param(ci,cj),startY1(ci,cj),'r.')
                                plot([param(ci,cj) param(ci,cj)],[startY1(ci,cj) finalY1(ci,cj)],'k')
                                plot(param(ci,cj),finalY1(ci,cj)-0.18,'k^')
                            end
                        end
                    end                    
                end
             end
            
        end
    end
end
hlx=xlabel('Flow rate out of nest (\phi): ants per minute')
hly=ylabel('Equilibrium number of ants crossing bridge X per minute')
set(hlx,'FontSize',14);
set(hly,'FontSize',14);
axis([0 phi 0 phi])