%
%       MALTHUS
%
%       matlab code to evaluate dispersal due to convolution 
%       of a spatial dispersal probability with a population density
%       function.
%
%       p(x)    -       population density in x
%       K(x)    -       probability of individual dispersal to location
%                          x starting at location 0.
%       xl	-       maximum x coordinate (also - minimum)
%       dx	-       spatial step  delta x = 2*xl/np
%       np	-       number of points  in discretization for x
%       		   (should be a power of 2 for  efficiency of FFT)
%       c	-       center of spatial  dispersal
%       D	-       diffusion constant (per step)
%       r	-       per-capita growth constant
%       
       
np=128; xl=5; dx=2*xl/np; c=0; D=.25; 
r=.8*ones(1,np)+.5*(abs(x)<=1.5);
nsteps=20;

%
%       discretize space, define an initial population distribution
%

x=-xl+[0:np-1]*dx;
p=(abs(x+3)<=.5);
p0=p;   %  keep track of initial population for comparison
plot(x,p,'r'), axis([-xl xl 0 1.6])

%
%       define a dispersal kernel
%

K=exp(-(x-c).^2/(4*D))/sqrt(4*pi*D);

%
%       normalize K to make it have integral one
%

K=K/(dx*sum(K));

%
%       calculate the fft of K, multiplying by dx to account
%       for the additional factor of np and converting from a
%       interval length of 1 to 2*xl.  The fftshift accounts for
%       using an interval of (-xl,xl) as opposed to (0,2*xl).
%

fK=fft(fftshift(K))*dx;

%       
%       begin iterating the model 
%       

for j=1:nsteps,

%
%       calculate the fft of f; perform the convolution
%

   fp=fft(p);
   fg=fK.*fp;

%
%       fg now contains the fourier transform of the convolution;
%       invert it, multiply by post-dispersal reproduction,
%        and take a look.
%

   g=r.*abs(ifft(fg));
   hold on, plot(x,p0,'r',x,g,'b'), axis([-xl xl 0 1.6]), hold off
   pause(.1)

%
%       update p and move to the next time step
%
   p=g;

end     % of the iterations
