function PredatorPreySimulation

%The following model uses six matrices, N1OLD, N2OLD, N3OLD, N1NEW, N2NEW
%and N3NEW.  The OLD matrices hold the population densities.
%The NEW matrices hold the dN values each iteration and are distributed
%over the OLD matrices using a 5x5 Gaussian distribution.
%The resources growth rate(G) and the 5x5 Gaussian distribution are
%precalculated and loaded into matrices G and rd
%RK4 calculations are made for all populations and populations are required to
%be >0

%Flow Chart
%  Step 1   Define variables
%  Step 2   Zero matrices
%  Step 3   Precalulate G matrix
%  Step 4   Precalulate Gaussian 5x5 matrix
%  Step 5   Initial conditions:  Populate N1, N2 and N3 matrices
%  Step 6   Update time
%  Step 7   Calculate number of new individuals, place them into NEW arrays
%  Step 8   Use Gaussian model to distribute new individuals in OLD matrices
%  Step 9   Capture population information for (10,10).
%  Step 10  (Return to Step3)
%  Step 11  Set plotting variables and plot output

%Define Variables  Step 1  DT = 18/75 is stable

TMAX=20;           %TMAX is the maximum time (sec)
M=30;             %M is the number of time increments
DT=TMAX/M;         %DT is the time increment (sec)

N=19;              %NxN is the cell matrix
DX=10/N;           %DX is the length of the side of a cell

N10=1;             %Initial N1 population for populated cells
N20=.25;           %Initial N2 population for populated cells
N30=1;             %Initial N3 population for populated cells
g1=2;
g2=-0.5;
r1=2;
r2=1.2;
p12=0.3;
p21=-1.5;
p31=0.7;
p32=0.5;

%Zero matrices Step2

N1OLD=zeros(N,N);    %Step2.  Initialize the P matrices.  These will contain
N1NEW=zeros(N,N);    %the populations of each cell N1(x,y,t), N2(x,y,t)and
N2OLD=zeros(N,N);    %N3(x,y,t). Generally, the OLD matrices will hold
N2NEW=zeros(N,N);    %information for t,
N3OLD=zeros(N,N);    %while the NEW matrices will hold information for t+1
N3NEW=zeros(N,N);
N1PLOT=zeros(M);
N2PLOT=zeros(M);
N3PLOT=zeros(M);
T=zeros(M);

%Precalulate G Step 3
G=zeros(N,N);
for i=1:N
    for j=1:N
%         G(i,j)=0.5;                                     %for problem 2
        G(i,j)=cos((i-1)*pi/(2*N-2))*sin((j-1)*pi/(N-1)) ; %for Problem 3 etc
    end
end

%Precalculate Gaussian distribution for a 5x5 matrix constants Step 4
c=erf(1/sqrt(2));
b=0.5*(erf(3/sqrt(2))-erf(1/sqrt(2)));
a=0.5-(b+0.5*c);

vecta=[a b c b a]';
vectb=vecta';

rd=vecta*vectb;
%  end of Gaussian precalculation

%Step 5:  Initial Condition: Populate N1, N2 and N3 into OLD matrices

for i=1:N;          %Initial Population for Resource N3
    for j=1:N;
        N3OLD(i,j)=N30;
    end
end

for i=1:N;          %Initial Population for Prey N1
    for j=1:N;
        N1OLD(i,j)=N10;
    end
end

for i=1:3;          %Initial Population for Predator N2
    for j=N/2-.5:N/2+.5;
        N2OLD(i,j)=N20;
    end
end

%set plotting variables
xg=DX/2:10/N:10;
yg=DX/2:10/N:10;

for k=1:M;         %time loop  Step 6

    %Step 7.  Calculate new individuals

    for i=1:N
        for j=1:N
            N1=N1OLD(i,j);
            N2=N2OLD(i,j);
            N3=N3OLD(i,j);

            %Runge-kutta 4  for N1  (prey)
            k1=(g1*N1-N1*p12*N2-r1*N1*N1/N3);
            y=N1+k1*DT/2;
            k2=(g1*y-y*p12*N2-r1*y*y/N3);
            y=N1+k2*DT/2;
            k3=(g1*y-y*p12*N2-r1*y*y/N3);
            y=N1+k3*DT;
            k4=(g1*y-y*p12*N2-r1*y*y/N3);
            N1NEW(i,j)=DT/6*(k1+2*k2+2*k3+k4);

            %Runge-kutta 4  for N2  (predators)
            k1=(g2*N2-N2*p21*N1-r2*N2*N2/N3);
            y=N2+k1*DT/2;
            k2=(g2*y-y*p21*N1-r2*y*y/N3);
            y=N2+k2*DT/2;
            k3=(g2*y-y*p21*N1-r2*y*y/N3);
            y=N2+k3*DT;
            k4=(g2*y-y*p21*N1-r2*y*y/N3);
            N2NEW(i,j)=DT/6*(k1+2*k2+2*k3+k4);

            %Runge-kutta 4  for N3  (resources)
            k1=(G(i,j)-N3*(p31*N1+p32*N2));
            y=N3+k1*DT/2;
            k2=(G(i,j)-y*(p31*N1+p32*N2));
            y=N3+k2*DT/2;
            k3=(G(i,j)-y*(p31*N1+p32*N2));
            y=N3+k3*DT;
            k4=(G(i,j)-y*(p31*N1+p32*N2));
            N3NEW(i,j)=DT/6*(k1+2*k2+2*k3+k4);
        end
    end

    %Step 8
    %Use Gaussian model to redistribute new individuals into the
    %NEW matrices

    for l=1:N              %(l,m) is the index of the source cell being distributed
        for m=1:N          % (i,j) is the index of the destination cell
            for i=l-2:l+2;              %  (xi,yi) is the 5x5 Gaussian index
                for j=m-2:m+2       %we only need to consider cells +-2
                    xi=3-l+i;     %Calculate Gaussian indexes (for rd matrix)
                    yi=3-m+j;
                    if ((i<1)||(j<1)||(i>N)||(j>N))     %edge handling
                        N1OLD(l,m)=N1OLD(l,m)+ N1NEW(l,m)*rd(xi,yi);
                        if N1OLD(l,m)<0, N1OLD(l,m)=0; end
                        N2OLD(l,m)=N2OLD(l,m)+ N2NEW(l,m)*rd(xi,yi);       %Step 5
                        if N2OLD(l,m)<0, N2OLD(l,m)=0; end
                        N3OLD(l,m)=N3OLD(l,m)+ N3NEW(l,m)*rd(xi,yi);
                        if N3OLD(l,m)<0, N3OLD(l,m)=.001; end

                    else                                 %normal Gaussian distribution
                        N1OLD(i,j)=N1OLD(i,j)+ N1NEW(l,m)*rd(xi,yi);
                        if N1OLD(i,j)<0, N1OLD(i,j)=0; end
                        N2OLD(i,j)=N2OLD(i,j)+ N2NEW(l,m)*rd(xi,yi);       %Step 5
                        if N2OLD(i,j)<0, N2OLD(i,j)=0; end
                        N3OLD(i,j)=N3OLD(i,j)+ N3NEW(l,m)*rd(xi,yi);
                        if N3OLD(i,j)<0, N3OLD(i,j)=.001; end              %Prevent /0
                    end                                                     %if/else/end
                end
            end
        end
    end

    %  Step 9  capture population data at (1,9)

    N1PLOT(k)=N1OLD(1,9);   %Capture data at x=1,y=9 for each time increment
    N2PLOT(k)=N2OLD(1,9);
    N3PLOT(k)=N3OLD(1,9);

    figure(1)
    surf (xg,yg,N1OLD);
    %axis([0 N*DX 0 N*DX])
    zlabel('N1');
    tt=k*DT;
    title(['Prey Population (Time = ', num2str(tt,2),')'])
    F1(k)=getframe;

    figure(2);
    surf (xg,yg,N2OLD);
    zlabel('N2');
    tt=k*DT;
    title(['Predator Population (Time = ', num2str(tt,2),')'])
    F2(k)=getframe;

 
    figure(3);
    surf (xg,yg,N3OLD);
    zlabel('N3');
    title('Resource Population')
    F3(k)=getframe;

end                % end of time loop (k)Step 10

%****************************************************************************
%  Output and plotting  Step 11

%Plotting Varibles     %Build T
for i=1:M
    T(i)=DT*i;
end

%Plot surface plot of steady state densities for problem 

% figure(4);
% surf (xg,yg,N1OLD);
% zlabel('N1');
% title('Prey Population')
%
% figure(5);
% surf (xg,yg,N2OLD);
% zlabel('N2');
% title('Predator Population')
%
% figure(6);
% surf (xg,yg,N3OLD);
% zlabel('N3');
% title('Resource Population')
%
% figure(7);
% surf (xg,yg,G);
% zlabel('G');
% title('Resource Replenish Rate')

 figure(1);
 movie(F1,1,4);
 movie2avi(F1,'prey', 'fps', 4)
 
 figure(2);
 movie(F2,1,4);
 movie2avi(F2,'predator', 'fps', 4)
 
 figure(3);
 movie(F3,1,4);
 movie2avi(F3,'resources', 'fps', 4)

 
figure(4)
subplot(3,1,1);  plot(T, N1PLOT);
xlabel('time(sec)')
ylabel('Prey pop (N1)')
title('N1 Population at (j=10,k=10)')
subplot(3,1,2);  plot(T, N2PLOT);
xlabel('time(sec)')
ylabel('Predator pop (N2)')
title('N2 Population at (j=10,k=10)')
subplot(3,1,3);  plot(T, N3PLOT);
xlabel('time(sec)')
ylabel('Resource pop (N3)')
title('N3 Population at (j=10,k=10)')

end             %end of function