function Lotka_Volterra_Phase_plane(alpha)
% This function shows the Prey-Predator dynamics
% as described by the level curves of the conserved quantity H, 
% α *x + y − log(x^α*y) = H
% α is a constant
% x represents the quantity of preys
% y represents the quantity of predator
%
%
% For more details, the reader can refer to 
% Mathematical Biology: I. An Introduction,
% J.D. Murray
%

if nargin == 0
    alpha = 1;
end

%% Define step size and axis range
stepsize_c = 0.01;
stepsize_q = 0.5;
max_u_size = 8;
max_v_size = 8;

%% Generate the grid points and compute the Lotka Volterra function
[Xc,Yc]=meshgrid(0:stepsize_c:max_u_size,0:stepsize_c:max_v_size);
Zc = log(alpha*Xc+Yc-log(Xc.^alpha.*Yc));

%% Create graphic
figure;
box on
hold on
grid on
axis equal
axis([0 max_u_size 0 max_v_size]);
contour(Xc,Yc,Zc,0.5:0.1:4);

% [C,Hc] = contour(Xc,Yc,Zc,0.5:0.1:4);
% RR = get(Hc,'Children');
% for i=1:numel(RR)
%     XD = get(RR(i),'XData');
%     YD = get(RR(i),'YData');
%     [XD,YD]
% end

% [PXc,PYc] = gradient(Zc,stepsize_c,stepsize_c);
%
% selectX = ismember(0:stepsize_c:max_u_size,...
%                    stepsize_q:stepsize_q:max_u_size-stepsize_q);
% selectY = ismember(0:stepsize_c:max_v_size,...
%                    stepsize_q:stepsize_q:max_v_size-stepsize_q);
% [Xq,Yq]=meshgrid(stepsize_q:stepsize_q:max_u_size-stepsize_q,...
%                  stepsize_q:stepsize_q:max_v_size-stepsize_q);
% quiver(Xq,Yq,PXc(selectX,selectY),PYc(selectX,selectY));

%% Create label and title
xlabel('x')
ylabel('y    ','Rotation',0);
title(['\alpha = ' num2str(alpha)],'Interpreter','Tex');
filename = [mfilename '_' sprintf('%03d',floor(100*alpha))];
saveas(gcf,filename,'png');
plot2svg([filename '.svg']);