%
%       BCs
%
%       matlab code to illustrate the effect of  boundary conditions 
%       and their implementation
%       
%       boundary conditions will be implemented  at x=0, x=xl
%       by reflecting the domain
%       
%
%       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)
%       D	-       diffusion constant (per step)
%       r	-       per-capita growth constant
%       
       
np=128; xl=5; dx=2*xl/np; c=0; D=.25; 
r=1;
nsteps=15;

%
%       discretize space, define an initial population distribution
%

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

%
%       define a dispersal kernel
%

K=.5*(exp(-(x-c).^2/(4*D))/sqrt(4*pi*D)+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,
%
%       skew-reflect the data to implement wall of doom
%       pa contains the actual population before dispersal
%       pr contains the reflected populuation before dispersal
%
   pa=p;
   pr=p(np:-1:1);

%
%       calculate the fft of the pa, pr; perform the convolutions
%

   fpa=fft(pa);
   fpr=fft(pr);
   fgr=fK.*fpr;
   fga=fK.*fpa;
   fg=fga+fgr;
   
%
%       fg now contains the fourier transform of the convolution;
%       invert it, multiply by post-dispersal reproduction,
%        and take a look.
%

   gr=real(ifft(fgr));
   ga=real(ifft(fga));
   g=real(ifft(fg));
   %plot(x,pr,'y--',x,gr,'y',x,pa,'b--',x,ga,'b',x,g,'g')
   hold on, plot(x,p,'r',x,gr,'y',x,ga,'b',x,g,'g--'), axis([-xl xl -1 1]), hold off
   pause(.1)

%
%       update p, undoing the reflection and move to the next time step
%
   p=r*g.*(x>=0);

end     % of the iterations
