%
%  BC2D
%  
%  a script using reflections to implement boundary conditions in 2 dims
%  reflecting (solid) boundary in y
%  lethal (absorbing) boundary in x
%  
%  dispersal is implemented using Allen's `ripple' function
%  

     xl=10; yl=10; np=64; dx=2*xl/np; dy=2*yl/np; % parameters for the
                                                   % x,y coordinates
                                                   
     x=linspace(0,xl-dx,np); x2=linspace(-xl,xl-dx,2*np);
     y=linspace(0,yl-dy,np); y2=linspace(-yl,yl-dy,2*np);

     [X,Y]=meshgrid(x,y);                          % grid of interest
     [X2,Y2]=meshgrid(x2,y2);                      % reflection grid

     t=2; D=.1; v=2;                                % set parameters

     K2D=exp(-(sqrt(X2.^2+Y2.^2)-v*t).^2./(4*D*t));   % ripple dispersal
     K2D=K2D/(dx*dy*trapz(trapz(K2D)));            % normalize

     fK=fft2(K2D)*dx*dy;                           % FFT for convolution
     p0=(((X-xl/2).^2+(Y-yl/2).^2)<=4);
     py=[p0; p0(np:-1:1,:)];                       % reflect in y
     px=[-py(:, np:-1:1) py];                      % reflect in x
     
     fp2=fft2(px);
     p2=real(fftshift(ifft2(fp2.*fK)));
     pt=p2(np+1:2*np,np+1:2*np);
     pcolor(X,Y,pt), shading interp, colormap hot, axis square