%3. gyakorlat: Lotka-Volterra rendszer
%dinamikai rendszer, az egyenlet rajzol:
%egyenletek: dx=x*(8-x-2y); dy=y*(5-x-y);
figure;
%ábrázolási határok:
max_pontx=9;%20;%
max_ponty=6;%20;%
min_pont=0.;
%
lepes=0.1;% lépésköz
x=min_pont:lepes:max_pontx;%
y=min_pont:lepes:max_ponty;%
%}
%{
% [0,0] környékén vizsgál:
lepes=0.01;
x=min_pont:lepes:0.5;
y=min_pont:lepes:0.5;
%}

%
% rácsot képez x,y alapján:
[X,Y]=meshgrid(x,y);

dx=X.*(8-X-2*Y);
dy=Y.*(5-X-Y);
dx=dx./ sqrt(dx.^2+dy.^2); %vektormezõ lenormálva
dy=dy./sqrt(dx.^2+dy.^2);
quiver(X,Y,dx, dy, 'k'); %vektormezõ kirajzolás
% az x y rácspontokon az x és az y irányba változást számolja ki és az
% eredményül kapott vektort ábrázolja az adott pontból kiindulva,
%}

y0=4-x/2; % 8-x-2*y=0 egyenes egyenlete
x0=5-y; % 5-x-y=0 egyenes egyenlete
hold on;
plot(x, y0,'b', zeros(1, length(y)), y, 'b', x0, y, 'r', x, zeros(1, length(x)), 'r'); %ábrázolom azokat az egyeneseket ahol dx=0 vagy dy=0
xlim([-0.1, max_pontx]); %ábrázolás határa
ylim([-0.1, max_ponty]);
title('egyenletek: dx=x*(8-x-2*y); dy=y*(5-x-y)'); %ábra cím
xlabel('x faj mennyisége'); %tengelyfeliratok
ylabel('y faj mennyisége');
legend('dx=0, mert x=0', 'dx=0, mert 8-x-2*y=0', 'dy=0, mert y=0', 'dy=0, mert 5-x-y=0'); %jelmagyarázat
%legend('vektormezõ gradiense adott pontokban', 'dx=0, mert x=0', 'dx=0, mert 8-x-2*y=0', 'dy=0, mert y=0', 'dy=0, mert 5-x-y=0');

% rács pontjaiból indulok:
%sx=reshape(X, 1, length(x)*length(y)); %kezdeti x pontok sorba rendezve
%sy=reshape(Y, 1, length(x)*length(y)); %kezdeti y pontok sorba rendezve


N=100; %10; %hány pontot számolok ki
%
%intervallum határok
sxk=0.1;%1;%0;%
syk=0.1;%3;%1;%
sxv=0.1;%2;%1;%
syv=max_ponty;%2;%0;%

%az intervallum, ahonnan indítom a trajektóriákat
%sx=linspace(sxk,sxv,N);
%sy=linspace(syk,syv,N);

sx2k=0.1;%2;%0;%
sy2k=0.1;%4;%6;%
sx2v=max_pontx;%3;%9;%
sy2v=0.1;%3;%0;%

%2 pont készletet fûzök egymáshoz
%sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N)];
%sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N)];

sx3k=max_pontx;
sy3k=0.1;
sx3v=max_pontx;
sy3v=max_ponty;

sx4k=0.1;
sy4k=max_ponty;
sx4v=max_pontx;
sy4v=max_ponty;

% az ábra határain 4 ponthalmaz
sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N), linspace(sx3k,sx3v,N), linspace(sx4k,sx4v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N), linspace(sy3k,sy3v,N), linspace(sy4k,sy4v,N)];
%}

%{
%idõben vissza: f(x)=-f(x)
%a trajektória végpontjait adom meg, a futás során visszafele halad a folyamat a kiindulási hely felé
sx5k=2.3;%
sy5k=2.7;%
sx5v=2.6;%
sy5v=2.85;%

%sx=linspace(sx5k,sx5v,N);
%sy=linspace(sy5k,sy5v,N);

sx6k=1.6;%
sy6k=3.1;%
sx6v=1.9;%
sy6v=3.25;%

sx=[linspace(sx5k,sx5v,N), linspace(sx6k,sx6v,N)];
sy=[linspace(sy5k,sy5v,N), linspace(sy6k,sy6v,N)];
%}

%{
%a kiiindulási pontok egy kör mentén:
%a körök középpontjai az egyensúlyi pontok
    mg1=[0;0];
    mg2=[8;0];
    mg3=[0;5];
    mg4=[2;3];
   
    Nl=100; %pontok száma körvonalon
    phi=linspace(0,2*pi,Nl); %szögek
    r=0.05; %a kör sugara
    
     kor=repmat(mg4,1,length(phi))+[r*sin(phi); r*cos(phi)]; %a kör kp-ját sokszorosítom  a repmattal és módosítom az irány és sugár alapján
     %A populáció mérete nem lehet negatív, ezért ezeket ki kell zárnunk
     kor_jo=kor(:,kor(1,:)>=0 & kor(2,:)>=0);
     sx=kor_jo(1,:);
     sy=kor_jo(2,:);
%}

%
tint=[0 10]; %futási idõsáv
for i=1:length(sx) %minden kiindulási pontra kiszámolom a diffegyenletet
    xy0=[sx(i); sy(i)]; %mindig az adott x,y pontból számolom ki a diffegyenletet (a hold on be van kapcsolva így mindent egyszerre rajzolok ki)
    [t, xy] = ode45(@(t, xy) [xy(1)*(8.5-xy(1)-2*xy(2)); xy(2)*(5.5-xy(1)-xy(2))],tint,xy0); %az egyenletet kiszámolom
    %[t, xy] = ode45(@(t, xy) -[xy(1)*(8-xy(1)-2*xy(2)); xy(2)*(5-xy(1)-xy(2))],tint,xy0); %az idõ visszafelé: f(x)=-f(x)
    %plot(xy(:,1),xy(:,2), 'g');
    %
    %a szeparatrix ábrázolása: ha az egyik stabil pontba fut be a mg akkor az egyik, ha a másikba, akkor a másik színnel ábrázolom a trajektóriát
    if abs(xy(end,1)-0)<1 %melyik egyensúlyi helyzethez van adott közel (az x tengelyen) a trajektória végpontja
        plot(xy(:,1),xy(:,2), 'g');
    elseif abs(xy(end,1)-8)<1
        plot(xy(:,1),xy(:,2), 'c');
    end
    %}  
end
%}

hold off;
