function nb2d
% 2D Spatial Nicholson-Bailey model
%
% The simple Nicholson-Bailey model can be extended to incorporate space.
% To this end we consider a spatial grid on which the Nicholson-Bailey
% dynamics take place. At each site (i, j) in the grid the dynamics are
% (more or less) given by the simple Nicholson-Bailey model, but in
% addition we allow that hosts and parasitoids disperse to all immediately
% neighbouring sites.

N = 100;         % spatial grid size
end_time = 1000; % number of time steps
dh = 0.1;        % host dispersal rate
dp = 0.1;        % parasitoid dispersal rate
R = 2;           % number of offspring of an unparasitized host surviving
                 % to the next year.
a = 1;           % a proportionality constant
c = 1;           % number of parasitoids that hatch from an infected host

% Ht and Pt represent the densities of the host and parasitoid population
% at year t.
Ht = zeros(N, N);
Pt = zeros(N, N);

% Auxiliary matrices for H* and P*
Hsum = zeros(N, N);
Psum = zeros(N, N);

% Left and right index, used for periodic boundary conditions
left = zeros(N);
left(1) = N;
for i=2:N
    left(i) = i - 1;
end

right = zeros(N);
right(N) = 1;
for i=1:N-1
    right(i) = i + 1;
end

% Initial conditions
for i=1:N
    for j=1:N
        Ht(i,j) = 1;
    end
end
Pt(N/2, N/2) = 1;
    
for t=1:end_time

    % Calculate sum over all 8 neighbouring fields
    for i=1:N
        for j=1:N            
            Hsum(i,j) = Ht(left(i),  left(j)) + ...
                        Ht(i,        left(j)) + ...
                        Ht(right(i), left(j)) + ...
                        Ht(left(i),  j) + ...
                        Ht(right(i), j) + ...
                        Ht(left(i),  right(j)) + ...
                        Ht(i,        right(j)) + ...
                        Ht(right(i), right(j));
            Psum(i,j) = Pt(left(i),  left(j)) + ...
                        Pt(i,        left(j)) + ...
                        Pt(right(i), left(j)) + ...
                        Pt(left(i),  j) + ...
                        Pt(right(i), j) + ...
                        Pt(left(i),  right(j)) + ...
                        Pt(i,        right(j)) + ...
                        Pt(right(i), right(j));
        end
    end
    
    % Calculate auxiliary variables H* and P*
    Hx = (1 - dh)*Ht + dh/8 * Hsum;
    Px = (1 - dp)*Pt + dp/8 * Psum;
    
    % Calculate new host and parasitoid populations
    expp = exp(-a*Px);
    Ht = R * Hx .* expp;
    Pt = c * Hx .* (1 - expp);
end

% Display host and parasitoid populations
figure('Units', 'normalized', 'Position', [0.1 0.1 0.5 0.25]);

subplot(1, 2, 1);
pcolor(log(Ht));
colormap(gray(100));
shading interp;
title('Logarithm of host population');

subplot(1, 2, 2);
pcolor(log(Pt));
colormap(gray(100));
shading interp;
title('Logarithm of parasitoid population');
