% Turing morfogenizis:
% dU/dt= 4*U+2*V+1*dU^2/dx^2
% dV/dt= -8*V-26*U+17*dV^2/dx^2
% 2 gradiens metszéspontja: pi/2
%feltételek mintázat képzõdéshez:
% a két diff nem lehet egyforma gyors
% nem lehet szimetrikus az [a, b; c, d] mátrix 
% legalább 2D-s rendszer kell (u,v)

dt = 0.00002; %idõbeli lépésköz %D*2*dt/dx^2<1
dx = pi/100; %térbeli lépésköz

x = 0:dx:pi; %viszgált pontok az egyenesen [0 pi]
nx = length(x);
nt = 200000; %idõtartomány vége

%diffúziós együtthatók: diffúzió sebesslge
%a diffúzió destabilizál
Du=1;
Dv=17;%1; %6; % elhal a rsz: 0-ra áll be %17; %stabil mintázat lesz 

% kölcsönhatások, hatások együtthatói:
cuu=4;   %u hat u-ra
cvu=2;   %v hat u-ra
cvv=-8;  %v hat v-re
cuv=-26; %u hat v-re

u = -1*ones(nx,1); %a minta kezdeti értékei
u(round(1:round(nx/4)))=1;
v = zeros(nx,1); %a minta kezdeti értékei

%
% skálázott színek mintázathoz, colormap-et definiál
pic=zeros(1, nx, 3); %3D-s rgb colormap
pic(:,:,1)=(u'+1)/2; % piros (0 1 közt 255 érték) %(0 1 közé skálázva az értékek), mert elszáll a  rendszer
pic(:,:,3)=(v'+1)/2; % kék (0 1 közt 255 érték)

%a végén megjelenített összesítõ szín kép
bigpic=zeros(nt+1,nx,3);
bigpic(1,:,1)=(u'+1)/2;
bigpic(1,:,3)=(v'+1)/2;
%}

%ábrázolás:
figure;
subplot(2,1,1);
%kezdeti érték rajzol:
plot(x, u, 'r', x, v, 'b'); 
grid(gca, 'minor');
   % ylim([-50 50]);
    xlim([0 pi]);
    
subplot(2,1,2)
image(pic); %a 3D-s rgb mátrix színes képként kirajzolása

    waitforbuttonpress;
    
for t = 1:nt %minden idõpillanatban
    u2 = u; %u2 = értékek az elõzõ idõplillanatban; u = új értékek
    v2 = v; %v2 = értékek az elõzõ idõplillanatban; v = új értékek
    
    for i = 1:nx %minden pontban kiszámolom az aktuális értékét
    
    % Neumann perem    
    if i==1
    u_prev=u2(1);
    v_prev=v2(1);
    else
    u_prev=u2(i-1);
    v_prev=v2(i-1);
    end
    
    if i==nx
    u_next=u2(end);
    v_next=v2(end);
    else
    u_next=u2(i+1);
    v_next=v2(i+1);    
    end
    
                % Turing egyenletek: 
                u(i) = u2(i) +cuu*dt*u2(i) +cvu*dt*v2(i) +Du*dt/dx^2*(u_next-2*u2(i)+u_prev); %dU/dt=a*U+b*V+D*dU^2/dx^2 % U függõ+ V függõ +diff tag
                v(i) = v2(i) +cvv*dt*v2(i) +cuv*dt*u2(i) +Dv*dt/dx^2*(v_next-2*v2(i)+v_prev); %dV/dt=c*V+d*U+D*dV^2/dx^2 % V függõ+ U függõ +diff tag
    end
    
    %ábrázol:
    if mod(t, 100)==0 %minden x. lépésben
    subplot(2,1,1)
    plot(x, u, 'r', x, v, 'b');
    grid(gca, 'minor');
    %ylim([-50 50]);
    xlim([0 pi]);
    title(['lépés: ', num2str(t)]);
    
    %
    subplot(2,1,2)
    pic(:,:,1)=min(1, max(0, (u'+1)/2)); % (0 1) közé skálázva értékek, mert elszáll a  rendszer
    pic(:,:,3)=min(1, max(0, (v'+1)/2));
    image(pic);
    %}
    
    %waitforbuttonpress;
     pause(0.01);
    end
     
    bigpic(t+1,:,1)=min(1, max(0, (u'+1)/2)); % (0 1) közé skálázva értékek, mert elszáll a  rendszer
    bigpic(t+1,:,3)=min(1, max(0, (v'+1)/2));
end

% a mintázat alakulása az idõben:
figure;
image(bigpic);
ylabel ('idõ');

