% MODEL FOR EARTHWORM DISPERSAL -- by Christien H. Ettema, 26 April 2001 --
% --------------------------------------------------------------------------------
% This model simulates the dispersal of earthworms in an arable field. Dispersal is 
% possible in two life stages (cocoons and juveniles), and takes place by active, 
% random movement (juveniles) or by passive dispersal, in soil stuck to tractor wheels. 
%
% The life cycle is as follows: 
% In fall (t=n), adult earthworms reproduce (cocoons) and subsequently die in winter. 
%    The cocoons do not have legs or wings and are generally too heavy to be blown 
%    around by wind. However, some passive dispersal takes place because soil with 
%    cocoons gets stuck to tractor wheels and is deposited elsewhere. 
% In spring, the dispersed cocoons hatch, and juveniles emerge. 
%    The juveniles disperse both actively and passively, by random walks and tractor 
%    wheels, respectively. It is assumed that the probability of "tractor dispersal" 
%    is similar for juveniles and cocoons, but mortality associated with 
%    tractor dispersal (squashing) is accounted only for juveniles. 
% In the next fall, juveniles mature to adults (t=n+1). It is assumed that mating 
%    and reproduction is a local affair (in other words, adults don't disperse). 


np=128;
xl=50; 
dx=2*xl/np;
x=linspace(-xl,xl-dx,np);
[X,Y]=meshgrid(x,x);
ngens=10;

% Set passive dispersal parameters and calculate passive dispersal kernel Kp.
%
% For Kp, it is assumed that the tractor only drives parallel to the x-axis, and that
%   dispersal in the y-direction is only because of the width of the tractor wheels W. 
%   With every turn of the tractor wheels, soil clots may fall off, creating a 
%   probability that cocoons reach that specific location. 
%
% wc = wheelcircumference of tractor
% Dp = mean distance of passive dispersal by tractor, same for cocoons as for juveniles
% W = width of wheel
% Kp = passive dispersal kernel
%

wc=6;
Dp=10;
W=1;
Kp=0.5*(1+cos((2*pi.*X)./wc)).*exp(-abs(X)./Dp).*exp(-Y.^2./W);
mass=dx^2*trapz(trapz(Kp));												% Normalize Kp
Kp=Kp/mass;
figure(1), plot(x,Kp(np/2+1,:));
fKp=dx^2*fft2(Kp);															% Take FT of Kp

%
% Set active dispersal parameters and calculate active dispersal kernel Ka
% Da = mean distance of active dispersal (juveniles only)
%

Da=2.4;																			
Ka=(1/(4*pi*Da^2))*exp((-X.^2-Y.^2)./(4*Da^2));
mass=dx^2*trapz(trapz(Ka));												% Normalize Ka
Ka=Ka/mass;
fKa=dx^2*fft2(Ka);															% Take FT of Ka

%
% Define start population of earthworm adults 'A0' 
%

A=((X.^2+Y.^2)<=1);
A0=A;

%
% Start loop
%

for i=1:ngens
   
%   
% Calculate dispersed cocoons 'dc'
% f = fecundity: number of cocoons per adult earthworm
% q = fraction of cocoons and juveniles picked up by the tractor wheel
% dc = dispersed cocoons
%

f=90*ones(np); 
q=0.25;  																		
fA=fft2(A.*f.*q);										% Take FT of A*f*q
fdc=fA.*fKp;                   					% Multiply FT's of A and Kp
dc=real(fftshift(ifft2(fdc)));					% Invert to get dispersed cocoons
dc=dc+((1-q).*f.*A);									% Add fraction of cocoons that "stayed"

%
% Post-dispersal cocoon hatching
% h = hatching success
% juv = total hatched juveniles
%

h=0.7*ones(np);  
juv=dc.*h; 

%
% Calculate passively dispersed juveniles djp, using the same 'q' and dispersal 
%   kernel 'Kp' as for cocoons:
%

fjuv=fft2(q.*juv);									% Take FT of q*juv	
fdjp=fjuv.*fKp;										% Multiply FT's of juv and Kp
djp=real(fftshift(ifft2(fdjp)));					% Invert to get passively dispersed 
															%                         juveniles 
%
% Calculate actively dispersed juveniles dja, assuming random walk dispersal:
%

fjuv=fft2((1-q).*juv);								% Take FT of (1-q)*juv
fdja=fjuv.*fKa;										% Multiply FT's of juv and Ka
dja=real(fftshift(ifft2(fdja)));					% Invert to get actively dispersed 
															% 								 juveniles
%
% Add djp and dja for total dispersed juveniles, adding mortality for tractor-
%    dispersed juveniles
% 
% mdjp = mortality due to tractor "squashing" of passively dispersed juveniles
% dj = dispersed juveniles
%

mdjp=0.9;
dj=((1-mdjp).*djp)+dja;

%
% Growth of juveniles to adult, until a maximum carrying capacity Kmax
% mat = maturation rate
%

Kmax=300;  mat=0.25;
Anew=mat.*dj;
Anew=min(Anew,Kmax);
figure(2), pcolor(X,Y,Anew), shading flat, hold on, contour(X,Y,A0,'m'), colorbar;
pause(.5)
A=Anew;																			% Update A

end

















   






