%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         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= 1;                          % y kezdeti értéke       
tinit=0;                            % kezdeti idõ          
Max_Time= 30;                      % szimuláció ideje       
tint =[tinit Max_Time];
alfa = 8;                           % x természetes szaporulata   
alfa2=5;                            % y szaporodik
beta =-2;                           % y hatása x 
beta2=-1;                           % x hatása a y
gamma = 1;                          % x természetes halálozása
gamma2=1;                           % y természetes halála
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);

%{
%4.a egy egyenesrõl indítok trajektóriákat:
N=100; %hány pontot számolok ki egy egyenes mentén
sx5k=0;
sy5k=1;
sx5v=1;
sy5v=0;
%sx5=linspace(sxk,sxv,N);
%sy5=linspace(syk,syv,N);

%4.c az ábra 4 határoló egyeneseirõl indítok trajektóriákat:
sxk=0.1;
syk=0.1;
sxv=0.1;
syv=Y_end;

sx2k=0.1;
sy2k=0.1;
sx2v=X_end;
sy2v=0.1;

sx3k=X_end;
sy3k=0.1;
sx3v=X_end;
sy3v=Y_end;

sx4k=0.1;
sy4k=Y_end;
sx4v=X_end;
sy4v=Y_end;

sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N), linspace(sx3k,sx3v,N), linspace(sx4k,sx4v,N), linspace(sx5k,sx5v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N), linspace(sy3k,sy3v,N), linspace(sy4k,sy4v,N), linspace(sy5k,sy5v,N)];
%}

%7. 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é
% a nyeregponthoz közeli egyenesek
N=100;
sxk=1.9;
syk=3;
sxv=2;
syv=2.9;

sx2k=2;
sy2k=3.1;
sx2v=2.1;
sy2v=3;

sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N)];

figure

 %for alfa = 01:0.5:15
     %for beta2 = -2:0.1:2
        % 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
           
      
            dx= X.*((alfa-gamma*X+beta*Y));             % rácspontokban vektormezõ gradiensének kiszámítása 
            dy= Y.*((alfa2-gamma2*Y +beta2*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
            y1 =(-alfa+gamma*x)/beta;
            y2 =(alfa2+beta2*x)/gamma2;  
            plot(x, y1, 'k', x, zeros(1, length(x)), 'g', x, y2,'g', zeros(1, length(y)), y, 'k'); %ábrázolom azokat az egyeneseket ahol dx=0 vagy dy=0
            
            for i=1:length(sx)
            Xinit=sx(i);
            Yinit=sy(i);
                
            [t, xy] = ode45(@(t, xy) [xy(1)*(alfa-gamma*xy(1)+beta*xy(2)); xy(2)*(alfa2-gamma2*xy(2)+beta2*xy(1))],tint,[Xinit Yinit]); %az diffegyenletet kiszámolom
            plot(xy(:,1),xy(:,2),'red');
            
             %
            %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 plusz a másik irányba is kirajzolom az adott színnel:
            if abs(xy(end,1)-0)<1 %melyik egyensúlyi helyzethez van adott közel (az x tengelyen) a trajektória végpontja
            [t, xyv] = ode45(@(t, xy) -[xy(1)*(alfa-gamma*xy(1)+beta*xy(2)); xy(2)*(alfa2-gamma2*xy(2)+beta2*xy(1))],tint,[Xinit Yinit]); %az idõ visszafelé: f(x)=-f(x)
                plot(xy(:,1),xy(:,2), 'm', xyv(:,1),xyv(:,2), 'm');
            elseif abs(xy(end,1)-8)<1
            [t, xyv] = ode45(@(t, xy) -[xy(1)*(alfa-gamma*xy(1)+beta*xy(2)); xy(2)*(alfa2-gamma2*xy(2)+beta2*xy(1))],tint,[Xinit Yinit]); %az idõ visszafelé: f(x)=-f(x)
                plot(xy(:,1),xy(:,2), 'c', xyv(:,1),xyv(:,2), 'c');
            end
            %}
            %{
            pontossag=1;
            if (xy(end,1)-0)<pontossag
                plot(xy(:,1),xy(:,2),'m');
            elseif (xy(end,1)-8)<pontossag
                plot(xy(:,1),xy(:,2),'c');
            end
            %}
            %hold off;
            xlim([-0.1, X_end]);
            ylim([-0.1, Y_end]);
           title(['beta2=', num2str(beta2)]);
            xlabel('x faj mennyisége');                                         % tengelyfeliratok
            ylabel('y faj mennyisége');
            
            end
            
        %}
            %pause(0.1);
           % waitforbuttonpress;
           %         end
         %       end
       %     end
     %   end
   % end
 %end
%}







