%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         3. practice         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nonlinear dynamics systems  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         2019.10.01.         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   Lotka-Volterra systems    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

%dynamical systems, the equation draws:
%equations: 
% dx=x*(4-x-2y) 
% dy=y*(6-x-2y)

%parameters:
xb=4; %x birth
xd=1; %x death
xi=-2; %y effect on x
yb=6; %y birth
yd=1; %y death
yi=-2; %x effect on y

figure;
%borders:
max_pontx=7;
max_ponty=7;
min_pont=0.;

%
stepsize=0.1;
x=min_pont:stepsize:max_pontx;%
y=min_pont:stepsize:max_ponty;%
%}

%
% 1.: vector field
% grid generation based on x,y:
[X,Y]=meshgrid(x,y);
% calculate the gradient of the vector field in each grid point
dx=X.*(xb-xd*X+xi*Y) ; 
dy=Y.*(yb-yd*Y+yi*X);
l=sqrt(dx.^2+dy.^2);
dx=dx./ l; %normalisation
dy=dy./l;
%plot the vector field (quiver function):
quiver(X,Y,dx, dy, 'k');
% calculate the change of x and y in every X,Y grid ponts, and draw a
% vector from it starting from the X,Y point
%}

%2. plot the lines where dx=0 or dy=0
y0=(-xb+xd*x)/xi; %xb-xd*x+xi*y=0 line's equation
x0=(-yb+yd*y)/yi;%yb-yd*y+yi*x=0 line's equation
hold on; %plot the trajectories on each other
plot(x, y0,'c', zeros(1, length(y)), y, 'c', x0, y, 'g', x, zeros(1, length(x)), 'g'); %plot the lines where dx=0 or dy=0
xlim([-0.1, max_pontx]); %plot borders
ylim([-0.1, max_ponty]);
title(['equations: dx=x*(',num2str(xb),'-',num2str(xd),'*x+',num2str(xi),'*y); dy=y*(',num2str(yb),'-',num2str(yd),'*y+',num2str(yi),'*x)']); %plot title
xlabel('amount of x'); %labels of the axes
ylabel('amount of y');
legend('direction field in given points', ['dx=0 <-  ',num2str(xb),'-',num2str(xd),'*x+',num2str(xi),'*y=0'],'dx=0 <-  x=0', ['dy=0 <-  ',num2str(yb),'-',num2str(yd),'*y+',num2str(yi),'*x=0'], 'dy=0 <-  y=0');

% Plot trajectories:
%
N=100;
sx_start=2.6;
sy_start=0.7;
sx_end=2.7;
sy_end=0.6;

sx2_start=2.65;
sy2_start=0.75;
sx2_end=2.75;
sy2_end=0.65;

sx=[linspace(sx_start,sx_end,N), linspace(sx2_start,sx2_end,N)];
sy=[linspace(sy_start,sy_end,N), linspace(sy2_start,sy2_end,N)];
%}

%
tint=[0 10]; %time interval
%calculate the differential equation for every point
for i=1:length(sx) 
    xy0=[sx(i); sy(i)]; % calculate the equation from the given x,y point (the "hold on" is on, so plot the trajectories in the same figure)
    [t, xy] = ode45(@(t, xy) [xy(1)*(xb-xd*xy(1)+xi*xy(2)); xy(2)*(yb-yd*xy(2)+yi*xy(1))],tint, xy0); %calculate the equation
    [t2, xy2] = ode45(@(t, xy) -[xy(1)*(xb-xd*xy(1)+xi*xy(2)); xy(2)*(yb-yd*xy(2)+yi*xy(1))],tint,xy0); %reverse time: f(x)=-f(x); % the end points of the trajectories are given, the dinamics runs back to the initial point during the simulation
    xy=[flipud(xy2); xy]; % the 2 outputs (normal and reverse time) are concatenated to one continous trajectory (the reverse time is turned at first)
    
    %plot(xy(:,1),xy(:,2), 'r');
    
    %
    % plot the separatrix: if a trajectory goes to one fix point, plot it with one color, if it goes to the other, plot it with the other color 
    if abs(xy(end,1)-0)<1 % the end point of the trajectory is close to (in the x axis) which fix point
        plot(xy(:,1),xy(:,2), 'r');
    elseif abs(xy(end,1)-4)<1
        plot(xy(:,1),xy(:,2), 'b');
    end
    %}  
end
%}

hold off;
