% Trial#28022007
% Chaitanya Athale, EMBL Heidelberg
% Numerical code for solving the 2 reaction system of reaction and diffusion reduced to 1D.
%REF.: Miura and Maini (2004) Review: Periodic pattern formation in
%reaction-diffusion systems: An introduction for numerical simulation
function sol = turingDiff
close all; clear all;
%-- Model parameters
P=parPDE( 0 );
cellS=P(1);
numC=P(2);
%-- solver parameters
geom = 2; % 0:slab, 1:cylinder, 2:sphere
rspan = linspace(0.001,cellS,numC);   % r is in micrometer
tspan = linspace(0,500,50);           % t is in seconds
%--pde solver
sol=pdepe( geom, @kinppasePDE, @kinppaseIC, @kinppaseBC, rspan, tspan);
A=sol(:,:,1);
I=sol(:,:,2);
figure(1),
subplot(1,2,1),surf(rspan,tspan, A), xlabel('Radial dist.'),ylabel('Time'),zlabel('A'),title('Activ'), colorbar;
subplot(1,2,2),surf(rspan,tspan, I), xlabel('Radial dist.'),ylabel('Time'),zlabel('I'), title('Inhib'), colorbar;
%end profiles
%make a movie
for s=1:1:10
    figure(2), plot(rspan,A(1,:),'x-b'),legend('Activator');
    M(s)=getframe;
    figure(3),plot(rspan,I(1,:),'o-r'),legend('Inhibitor');
    N(s)=getframe;    
end
    
fn =11;
for n=1:1:size(A,1)
    figure(2), plot(rspan,A(n,:),'x-b'),legend('Activator'),text(cellS/2, 0, ['Time = ', sprintf('%02i',n) ]);
    M(fn)=getframe;
    figure(3),plot(rspan,I(n,:),'o-r'),legend('Inhibitor'),text(cellS/2, 0, ['Time = ', sprintf('%02i',n) ]);
    N(fn)=getframe;
    fn = fn+1;
end
movie2avi(M,'turingAct.avi','COMPRESSION','none');
movie2avi(N,'turingInh.avi','COMPRESSION','none');

function [c,f,s] = kinppasePDE( r, t, u, dudr )
% the return values
c = zeros(2,1);
f = zeros(2,1);
s = zeros(2,1);
p=parPDE( 0 );
%-- set the parameters
Da=p(3); %diff of A
Di=p(4); %diff of I
aa=p(5); %autoactiv of A
ai=p(6); %activ of I by A
ia=p(7); %inhib of A by I %7
ii=p(8); %inhib of I by I %8
%%---describe the varibles
A = u(1);
I = u(2);
%-- solution
At=aa*A - ia*I - A^3;
It=ai*A - ii*I;
%--
c = ones(2,1);
f = [ Da*dudr(1);Di*dudr(2)];
s = [At; It];
%%%-- initial conditions
function u0 = kinppaseIC(r)
P = parPDE( 0 );
%--
u0 = zeros(2,1);
u0(1) = rand;%/nCell;
u0(2) = rand;
%-- Boundary conditions
function [pl,ql,pr,qr] = kinppaseBC(xl,ul,xr,ur,t)
%%values
pl = zeros(2,1);
pr = zeros(2,1);
%%function (flux)
ql = ones(2,1);
qr = ones(2,1);
%-- Function to call and change parameters
function P = parPDE( n )
%-- use n to increment any given parameter
cellSize = 1;       %1
nCell =100;         %2
Da=0.0002;          %3
Di=0.01;            %4
aa=0.6; %autoacr A  %5
ai=1.5;%act I by A %6
ia=1;  %inhib of A by I %7
ii=2;  %inhib of I by I %8
P=[cellSize; nCell; Da; Di; aa; ai; ia; ii];