function [] = spg1()
%SPG1  Reaktions-Diffusionsgleichung in 2D mit FTCS
%      s ... Substratkonzentration
%      p ... Bakterienzahl
%      g ... Pufferkonzentration

close all
clc
help spg1

N  = input(' Anzahl der raeumlichen Gitterpunkte N [40]: ');
if isempty(N)
   N = 40;
end
dt = input(' Zeitschritt dt [1.0e-2]: ');
if isempty(dt)
   dt = 1.0e-2;
end

L  = 20.0;                % x,y = -L/2 bis x,y = +L/2
dx = L/(N-1);             % Schrittweite im x-Ortsgitter
dy = L/(N-1);             % Schrittweite im y-Ortsgitter
D  = 0.1;                 % Diffusionskoeffizient

rx = D*dt/dx^2;
ry = D*dt/dy^2;

if ( (rx + ry) > 0.5 )
   disp(' WARNUNG: FTCS ist instabil')
end

x = (0:N-1)*dx - L/2;     % x(1) ,..., x(N)     : raeumliche Gitterpunkte x
y = (0:N-1)*dy - L/2;     % y(1) ,..., y(N)     : raeumliche Gitterpunkte y

s    = zeros(N, N);       % s(1,1) ,..., s(1,N) : Naeherungsloesung fuer s
                          % s(N,1) ,..., s(N,N)
sneu = zeros(N, N);

g    = ones(N, N);
gneu = zeros(N, N);

p    = ones(N, N);
pneu = zeros(N, N);

nstep = 10000;            % Anzahl der Zeitschritte

% Anfangsbedingung
t = 0.0;
s = f(x,y);               % s(x,y,0)= f(x,y)
% Graphik
figure(1)

subplot(3,1,1)
mesh( x, y, s )
view(40,30)
zlabel('s(x,y)')
title(['Substrat, t = ', num2str(t)])

subplot(3,1,2)
mesh( x, y, g )
view(40,30)
zlabel('g(x,y)')
title(['Puffer, t = ', num2str(t)])

subplot(3,1,3)
mesh( x, y, p )
view(40,30)
zlabel('p(x,y)')
title(['Bakterienzahl, t = ', num2str(t)])

% Postscript
epsfile = 'spg1a.eps'; print('-depsc', epsfile)

for istep = 1:nstep
    %
    % Explizites Verfahren fuer s und g
    %
    for i = 2:(N-1)
        for j = 2:(N-1)
            sneu(i,j) = s(i,j) + ...
                        rx * ( s(i+1,j  ) - 2*s(i,j) + s(i-1,j  ) ) + ...
                        ry * ( s(i  ,j+1) - 2*s(i,j) + s(i  ,j-1) ) - ...
                        dt * s(i,j)*g(i,j)*p(i,j) / (1.0 + s(i,j)*g(i,j));
            gneu(i,j) = g(i,j) + ...
                        rx * ( g(i+1,j  ) - 2*g(i,j) + g(i-1,j  ) ) + ...
                        ry * ( g(i  ,j+1) - 2*g(i,j) + g(i  ,j-1) ) - ...
                        dt * s(i,j)*g(i,j)*p(i,j) / (1.0 + s(i,j)*g(i,j));
        end
    end
    %
    % Euler-Schritt fuer p
    %
    for i = 1:N
        for j = 1:N
            pneu(i,j) = p(i,j) + ...
                        dt * s(i,j)*g(i,j)*p(i,j) / (1.0 + s(i,j)*g(i,j)); 
        end
    end
    %
    % Neumannsche Randbedingungen fuer s und g
    %
    sneu(:,1) = s(:,1) + 2.0 * rx * ( s(:,2)   - s(:,1) )  -  ...
                         dt * s(:,1).*g(:,1).*p(:,1)./(1.0+s(:,1).*g(:,1));
    sneu(:,N) = s(:,N) + 2.0 * rx * ( s(:,N-1) - s(:,N) )  - ...
                         dt * s(:,N).*g(:,N).*p(:,N)./(1.0+s(:,N).*g(:,N));
    sneu(1,:) = s(1,:) + 2.0 * ry * ( s(2,:)   - s(1,:) )  - ...
                         dt * s(1,:).*g(1,:).*p(1,:)./(1.0+s(1,:).*g(1,:));
    sneu(N,:) = s(N,:) + 2.0 * ry * ( s(N-1,:) - s(N,:) )  - ...
                         dt * s(N,:).*g(N,:).*p(N,:)./(1.0+s(N,:).*g(N,:));
    gneu(:,1) = g(:,1) + 2.0 * rx * ( g(:,2)   - g(:,1) )  -  ...
                         dt * s(:,1).*g(:,1).*p(:,1)./(1.0+s(:,1).*g(:,1));
    gneu(:,N) = g(:,N) + 2.0 * rx * ( g(:,N-1) - g(:,N) )  - ...
                         dt * s(:,N).*g(:,N).*p(:,N)./(1.0+s(:,N).*g(:,N));
    gneu(1,:) = g(1,:) + 2.0 * ry * ( g(2,:)   - g(1,:) )  - ...
                         dt * s(1,:).*g(1,:).*p(1,:)./(1.0+s(1,:).*g(1,:));
    gneu(N,:) = g(N,:) + 2.0 * ry * ( g(N-1,:) - g(N,:) )  - ...
                         dt * s(N,:).*g(N,:).*p(N,:)./(1.0+s(N,:).*g(N,:));
    s = sneu;
    g = gneu;
    p = pneu;
    t = t + dt;
end

% Graphik
figure(2), subplot(4,1,1), mesh( x, y, s )
view(30,60), xlabel('x'), ylabel('y'), zlabel('s(x,y)')
title(['Substrat, t = ', num2str(t)])
subplot(4,1,2), mesh( x, y, g )
view(30,60), zlabel('g(x,y)')
title(['Puffer, t = ', num2str(t)])
subplot(4,1,3), mesh( x, y, p )
view(30,60), zlabel('p(x,y)')
title(['Bakterienzahl, t = ', num2str(t)])
subplot(4,1,4), plot(x,p(N/2,:))
xlabel('x'), ylabel('p(x,y=0)')
epsfile = 'spg1b.eps'; print('-depsc', epsfile)

%--------------------------------------------------------------------------
function [s] = f(x,y)
% Anfangsbedingung
A=10.0; delta=0.1; s=zeros(length(x),length(y));
for i = 1:length(x)
    for j = 1:length(y)
        s(i,j) = A * exp( - delta * (x(i)^2 + y(j)^2) );
    end
end
%--------------------------------------------------------------------------