% Phytophthora model
% Wopke van der Werf, Walter Rossing, Geert-Jan Kessel, Hans Heesterbeek & Jim Powell
% 25 April 2001

%
% parameters
%

% space
nrows = 64;
nplants = 128;
rows = [-nrows/2:1:nrows/2-1];
plants = [-nplants/2:1:nplants/2-1];
row_distance = 0.75;				% meters
plant_distance = 0.30;			% meters

X = rows * row_distance;
Y = plants * plant_distance;
[x,y] = meshgrid(X,Y);

% crop; lai
laimax = 5;							% m^2
laigrowthrate = 0.10; 			% m^2/m^2/d

% disease parameters; sporulation intensity (spores/lesion area); interception_probability;
% infection_probability;

sporulation_intensity = 10000;	% spores/cm^2
infection_probability = 0.01;
interception_probability = 0.1;
psusceptible = .3333;

% dispersal parameters
a = 1; 								% 1/time
b = 1;								% -
c = 1.2;								% distance/time
										% c/a: mean dispersal distance

% time
ngens = 90;

%
%	set up mixture of susceptible and resistant potato plants
%

plants=rand(nplants,nrows);
resistent_plants=(plants>=psusceptible);
susceptible_plants=(plants<psusceptible);

%	Setting up a movie

M=moviein(ngens); 

%
%	Initial conditions
% 

lat1  = zeros(nplants,nrows);
lat2  = zeros(nplants,nrows);
lat3  = zeros(nplants,nrows);
lat4  = zeros(nplants,nrows);
lat5  = zeros(nplants,nrows);
latent = lat1 + lat2 + lat3 + lat4 + lat5;
nlinf = zeros(nplants,nrows);
lai   = zeros(nplants,nrows);
landed_sporangia = zeros(nplants,nrows);

% lat6=zeros(nplants,nrows);
% lat7=zeros(nplants,nrows);

% initial_inoculum=rand(nplants,nrows);
% source = initial_inoculum .* susceptible_plants;
% inf = (source>=.5);

inf   = zeros(nplants,nrows);
inf(nplants/2-1:nplants/2+1, nrows/2-1:nrows/2+1) = 10;

%
%	Define dispersal kernel for sporangia	
% 

not_normalized_kernel = 0.5 * b * a/c * (a/c * sqrt(x.^2 + y.^2)).^(b-1) .* exp(-a/c * sqrt(x.^2 + y.^2).^b);
integral = trapz(trapz(not_normalized_kernel * plant_distance * row_distance));
kernel = not_normalized_kernel/integral;
fkernel = fft2(kernel) * row_distance * plant_distance;

%
% Time loop
%
for j=1:ngens

%
%	Advance the life cycles using adjusted N-B model
%

lai = laimax * (1-exp(-laigrowthrate * ngens/laimax));

produced_sporangia = sporulation_intensity * inf;
nlinf = nlinf + inf;
inf   = lat5;
lat5  = lat4;
lat4  = lat3;
lat3  = lat2;
lat2  = lat1;
% lat1  = susceptible_plants .* infection_probability .* ((lai-(latent+nlinf+inf)/10000)./lai) .* landed_sporangia;
% lat1  = susceptible_plants .* infection_probability .* exp(-(latent+nlinf+inf)/10000)./lai)) .* landed_sporangia;
% avail = 10000 * lai .* (1-p);
avail = lai * 10000 - latent - inf - nlinf;
lat1  = avail .* susceptible_plants .* (1-exp(- infection_probability .* ...
   landed_sporangia .* avail ./ (lai .* 10000) ./ avail));

lesions = inf + nlinf;
latent = lat1 + lat2 + lat3 + lat4 + lat5;
p = .0001 * (latent + lesions)./lai;

%
% 	Do the dispersal step
%

fproduced_sporangia = fft2(produced_sporangia);

%	Now do the convolutions.
%	The fftshift serves to center the probability functions.
%

landed_sporangia = interception_probability * real(fftshift(ifft2(fkernel .* fproduced_sporangia)));

%
%	Now plot the solution
%
figure(1);pcolor(x,y,p), shading flat, caxis([0 1]), colormap(hot), axis square, colorbar;
% figure(2);pcolor(x,y,inf), shading flat, colormap(hot), axis square, colorbar;
% figure(3);pcolor(x,y,nlinf), shading flat, colormap(hot), axis square, colorbar;

%subplot(1,2,1)
%pcolor(X,Y,h0),shading flat, colormap hot, title('Host'),colorbar
%subplot(1,2,2)
%pcolor(X,Y,p0),shading flat, colormap jet, title('Parasite'),colorbar

M(:,j) = getframe; %	If you want movies

end
%	
%	Here ends the iteration
%

%surf(X,Y,h0-p0), shading interp
%  
%  To play the movie,  use this command:
%  movie(M,2,8)
