%
% The formation of spiral waves based on Murray 1989
% and the centered finite difference scheme with periodic 
% boundary conditions. 
%
% Mathematical Biology, Jae-Hun Jung March, 2007
%

clear all, close all
alpha = 0.1; 
%alpha=10;
beta = 0.01;% beta=1;
Du = 0.05; %Du=0.0
gamma =0.6; %gamma = 1;
Dv = gamma*Du;

N = 50; M = 50;

x = linspace(-1,1,N);
y = linspace(-1,1,M);

h = 2/(N-1); g = 2/(M-1); k = min(h,g)/100;

u1 = rand(M,N);    w1 = u1;
v =rand(M+2,N+2);  z = v;

u1=zeros(M,N);
for ix = 1:N
    for iy = 1:M
        w1(iy,ix) = sin(10*pi*y(iy))*sin(10*pi*x(ix));
    end
end

v(2:M+1,2:N+1) = u1;z(2:M+1,2:N+1) = w1;
v(2:M+1,1)   = u1(:,N-1);z(2:M+1,1)   = w1(:,N-1);
v(2:M+1,N+2) = u1(:,2);z(2:M+1,N+2) = w1(:,2);
v(1,2:N+1)   = u1(M-1,:);z(1,2:N+1)   = w1(M-1,:);
v(M+2,2:N+1) = u1(2,:); z(M+2,2:N+1) = w1(2,:); 

for i = 1:1000
    time = i*k; 
      tmpv = v(2:M+1,2:N+1);
      tmpw = z(2:M+1,2:N+1);
      u1(:,:) = tmpv + k*Du*(v(2:M+1,3:N+2)-2*v(2:M+1,2:N+1)+v(2:M+1,1:N))/h^2 ...
               +k*Du*(v(3:M+2,2:N+1)-2*v(2:M+1,2:N+1)+v(1:M,2:N+1))/g^2 ... 
               +k*(alpha*(1-tmpv)-tmpv.*tmpw.^2);
      w1(:,:) = tmpw+k*Dv*(z(2:M+1,3:N+2)-2*z(2:M+1,2:N+1)+z(2:M+1,1:N))/h^2 ...
               +k*Dv*(z(3:M+2,2:N+1)-2*z(2:M+1,2:N+1)+z(1:M,2:N+1))/g^2 ... 
               +k*(tmpv.*tmpw.^2-(alpha+beta)*tmpw./(1+tmpw+tmpv));
 
    v(2:M+1,2:N+1) = u1;      z(2:M+1,2:N+1) = w1; 
    v(2:M+1,1)   = u1(:,N-1); z(2:M+1,1)   = w1(:,N-1);
    v(2:M+1,N+2) = u1(:,2);   z(2:M+1,N+2) = w1(:,2);
    v(1,2:N+1)   = u1(M-1,:); z(1,2:N+1)   = w1(M-1,:);
    v(M+2,2:N+1) = u1(2,:);   z(M+2,2:N+1) = w1(2,:);
    
    imagesc(u1)
    A(i) = getframe
    pause(0.001)
end