% 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.01;	% activator diffusion coefficient
r_a = 0.08;	% activator removal rate
b_a = 0.01;	% basic activator production
s_a = 0.1;
D_b = 0.3;	% inhibitor diffusion coefficient
r_b = 0.0;	% inhibitor removal rate
b_b = 0.08;	% basic inhibitor production
s_b = 0.0;
D_c = 0.002;	% second inhibitor diffusion coefficient
r_c = 0.01;	% second inhibitor removal rate
b_c = 0.0;
s_c = 1.0;

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

% FIXME: Meinhardt hack
nT = 450;
% this seems useless
KP = 1;
N = 80;
dt = 1e-0;
h = 1e-0;

% activator production rate
s = r_a * (0.95 + 0.1 * 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;
b(1,:) = 1;
c(1,:) = 1;

A_a = 1;
A_b = 0;
A_c = 1;

% SP.BAS:109
a(10,:) = A_a;
a(22,:) = A_a;
a(55,:) = A_a;
a(75,:) = A_a;
a(115,:) = A_a;
a(180,:) = A_a;
a(195,:) = A_a;
a(240,:) = A_a;
a(310,:) = A_a;
a(352,:) = A_a;
a(375,:) = A_a;
a(425,:) = A_a;
%a(455,:) = A_a;
%a(555,:) = A_a;
%a(585,:) = A_a;
%a(595,:) = A_a;
%a(630,:) = A_a;

c(10,:) = A_c;
c(22,:) = A_c;
c(55,:) = A_c;
c(75,:) = A_c;
c(115,:) = A_c;
c(180,:) = A_c;
c(195,:) = A_c;
c(240,:) = A_c;
c(310,:) = A_c;
c(352,:) = A_c;
c(375,:) = A_c;
c(425,:) = A_c;
%c(455,:) = A_c;
%c(555,:) = A_c;
%c(585,:) = A_c;
%c(595,:) = A_c;
%c(630,:) = A_c;


% 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); % v- till here is aq
    temp_a = a_star/(s_b + s_c * oc) .* (s .* ob) - r_a .* oa;
    temp_b = b_b - a_star/(s_b + s_c * oc) .* (s .* ob) - r_b * ob;
    temp_c = r_c * (oa - oc);
    na = oa + dt * (temp_a + D_a*oa*Lap);
    nb = ob + dt * (temp_b + D_b*ob*Lap);
    nc = oc + dt * (temp_c + 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)])

norm(a(end,:))