% 1D time dependent activator-inhibitor system on periodic domain

L = 4;		% length of periodic domain
N = 256;	% space sample points
LifeT = 1;	% life time of the sea shell
nT = 256;	% time sample points

% system parameters
D_a = 0.005;	% activator diffusion coefficient
r_a = 0.0;	% activator removal rate
b_a = 0.0;	% basic activator production
D_b = 0.0;	% inhibitor diffusion coefficient
r_b = 0.0;	% inhibitor removal rate
b_b = 0.01;	% basic inhibitor production
D_c = 0.3;	% second inhibitor diffusion coefficient
r_c = 0.0;	% second inhibitor removal rate
b_c = 0.02;
s_a = 0.05;

% help variables
h = L/N;	% grid size
dt = LifeT/nT;	% time step

% FIXME: Meinhardt hack
nT = 144;
KP = 400;
N = 40;
dt = 1.0;
h = 1.0;

% activator production rate
s = r_a * (0.99 + 0.02 * rand(1,N));

% Build Laplace Operator (Finite Differences)
e = ones(N,1) / h^2;
Lap = spdiags([e -2*e e], -1:1, N, N);
Lap(N,1) = e(1);
Lap(1,N) = e(1);

% set initial values
a = zeros(nT, N);
b = zeros(nT, N);
c = zeros(nT, N);
a(1,:) = 1.0e-3 * rand();
b(1,:) = 1.0e-3 * rand();
c(1,:) = 1.0e-3 * rand();

% solve problem (Fixed Step Euler Forward)
for n=1:nT-1
  oa = a(n,:);
  ob = b(n,:);
  oc = c(n,:);
  for k=1:KP
    a_star = (oa.^2 + b_a)./(1 + s_a * oa.^2);
    aqb = s .* a_star .* ob;
    aqc = s .* a_star .* oc;
    na = oa + dt .* (aqb .+ aqc .- r_a .* oa + D_a*oa*Lap);
    nb = ob + dt .* (b_b .- aqb .- r_b .* ob + D_b*ob*Lap);
    nc = oc + dt .* (b_c .- aqc .- r_c .* oc + D_c*oc*Lap);
    oa = na;
    ob = nb;
    oc = nc;
  end
  a(n+1,:) = na;
  b(n+1,:) = nb;
  c(n+1,:) = nc;
end

subplot(1,3,1), imagesc(a);
subplot(1,3,2), imagesc(b);
subplot(1,3,3), imagesc(c);

function y = makecolor(x)
% makecolor: return a vector containing 64 steps of the actual color percentage
    if (x == 0)
        y = zeros(64,1);
    else
        y = [0:x/63:x]';
    end
end

colormap([makecolor(0.1), makecolor(0.7), makecolor(1)])