%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         4.gyakorlat         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nemlin. din. mod. a biológiában  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         2015.03.05.         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        Lotka-Volterra       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 1. Lotka-Volterra
%
% Paraméterek:
Xinit= 1;                           % x kezdeti értéke       
Yinit= 10;                          % y kezdeti értéke       
tinit=0;                            % kezdeti idõ          
Max_Time= 100;                      % szimuláció ideje       
tint =[tinit Max_Time];
alfa = 1;                           % préda természetes szaporulata   
beta = 1;                           % vadász hatása prédára (fogyasztja)
gamma = 1;                          % vadász természetes halálozása
delta = 1;                          % préda hatása a vadászra (gyarapítja)
N = 100;                            %
X_end = 15;                         % x tengely vége
Y_end = 15;                         % y tengely vége
x = 0:0.4:X_end;                    %                
y = 0:0.4:Y_end;                    % 
[X,Y] = meshgrid(x,y);              % hálót készít x,y vektorok alapján

%másik általánosabb rendszer:
% 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);


figure

% for alfa = 0.5:0.1:1
    % for beta = 0.5:0.1:1
        % for gamma = 0.5:0.1:1
           % for delta = 0.5:0.1:1
                 % for alfa2 = 0.5:0.1:1
                     % for gamma2 = 0.5:0.1:1
           
            subplot(2,2,2)
            dx= X.*((alfa-beta*Y));             % rácspontokban vektormezõ gradiensének kiszámítása 
            dy= Y.*((-gamma + delta*X));        % 
            dx= dx./sqrt(dx.^2+dy.^2);          % vektormezõ lenormálva (x)(szebb ábrához)
            dy= dy./sqrt(dx.^2+dy.^2);          %                       (y)
            % az x y rácspontokon a 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,
            quiver(X,Y,dx, dy, 'blue');
            hold on;
            
            %egyenesek kirajzolása, ahol dx=0 vagy dy=0
            y0 =alfa/beta*ones(1,length(x));  %x*(alfa-beta*y)=0 => x=0 vagy alfa-beta*y=0 -> y=alfa/beta
            x0 =gamma/delta*ones(1,length(y));  %y*(-gamma+delta*x)=0 => y=0 vagy gamma+delta*x=0 -> x=-gamma/delta  
            plot(x0, y, 'k', x, zeros(1, length(x)), 'k', x, y0,'g', zeros(1, length(y)), y, 'g'); %ábrázolom azokat az egyeneseket ahol dx=0 vagy dy=0
            
            [t, xy] = ode45(@(t, xy) [xy(1)*(alfa-beta*xy(2)); xy(2)*(-gamma+delta*xy(1))],tint,[Xinit Yinit]); %az diffegyenletet kiszámolom
            plot(xy(:,1),xy(:,2),'red');
            hold off;
            xlim([-0.1, X_end]);
            ylim([-0.1, Y_end]);
            title('egyenletek: dx=x*(alfa-beta*y); dy=y*(-gamma+delta*x)');     % ábra cím
            xlabel('x faj (préda) mennyisége');                                         % tengelyfeliratok
            ylabel('y faj (ragadozó) mennyisége');
            

            subplot (2,2,1);
            plot(t,xy(:,1),'green');
            hold on;
            plot(t,xy(:,2),'black');
            hold off;
            title(['alfa= ',num2str(alfa),'  beta= ',num2str(beta),'   gamma = ',num2str(gamma),'  delta = ',num2str(delta)]);
            xlabel('t - idõ');                                         % tengelyfeliratok
            ylabel('fajok mennyisége');
            legend('préda', 'ragadozó');
            
            
            subplot (2,2,3);
            plot3(xy(:,1),xy(:,2),t,'red');
            xlabel('préda mennyisége');
            ylabel('ragadozó mennyisége');
            zlabel('idõ');
        %}
            %pause(0.1);
            %waitforbuttonpress;
           %         end
         %       end
       %     end
     %   end
   % end
% end
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%{
% 2. Trace-det. diagram 
% a mátrixot változtatom, hogy változik a fázisportré (sé, det, trace jelölve)
x = -2:1:2;
y = -2:1:2;
[X Y] = meshgrid(x,y);

K = 500;
fix = zeros(K,length(x)*length(y));
fiy = zeros(K,length(x)*length(y));

figure;
%trace, det görbe kirajzol
subplot(1,2,2);
tr=(-8:0.1:8);
plot(tr, tr.^2./4, 'b', tr, zeros(size(tr)), 'b', zeros(size(tr)), tr, 'b');
xlabel('trace');
ylabel('det');
hold on;
%
% paraméter változtat
hatar=3;
for a = -hatar:0.2:hatar 
    A = [0 1; -4+a^2 2*a] %egyenletrendszer mátrixos alakban
   
se=eig(A);
%adott rendszerre difegyet megoldom:
 n=1;
for i = 1:length(x)
    for j = 1:length(y)
        %{
        if trace(A)>0 
           A1=-A; %idõ hátra (trajektória végpontjai a kijelölt pontok) ha elmennek a megoldások 
        else
            A1=A;
        end
        %}
        [t,fi] = ode45(@(t,y) A*y, linspace(0,5,K), [x(i) y(j)]);
        fix(:,n) = fi(:,1);
        fiy(:,n) = fi(:,2);
        n = n+1;
    end
end

%plotol trajektóriákat adott pontokból:
subplot(1,2,1);
plot(fix,fiy);
hold on;
plot(X,Y,'x');
xlim([2*min(x) 2*max(x)]);
ylim([2*min(y) 2*max(y)]);
title(['sajátértékek: ', num2str(se(1)), ', ', num2str(se(2))]);
xlabel('x');
ylabel('y');
hold off

%trace det kirajzol:
subplot(1,2,2);
plot(trace(A), det(A), 'r*');
title(['trace(A): ', num2str(trace(A)), ', det(A):', num2str(det(A))]);

waitforbuttonpress;
%pause(0.1)
end
%}





