%%
%  Author: Peter Polcz (ppolcz@gmail.com) 
%  Created on 2023. October 25. (2023a)
%

r = 1.5;
a = 5;
b = 18;
D = 2;
c = 4;

K_Hopf = b*(a*D+c)/(a*D-c);
K = K_Hopf * (0.8 + 0.4);

syms t x y real
f = [
    r*x*(1 - x/K) - y*(a*x)/(b+x)
    -c*y + D*y*(a*x)/(b+x)
    ];
x = [x;y];

LimX = K_Hopf*1.5;
LimY = K_Hopf;

Eq1 = [
    b*c/(a*D-c)
    b*D*r*(-K*c-b*c+K*a*D)/K/(a*D-c)^2
    ];

%%% 
% Generate anonymous functions
% 
f_fh = matlabFunction(f,'vars',{t,x});
J_fh = matlabFunction(jacobian(f,x),'vars',{x});
f1_fh = matlabFunction(f(1),'vars',x);
f2_fh = matlabFunction(f(2),'vars',x);

%%%
% Compute equilibria
%
sol = solve(f,x);
Eq = double([sol.x,sol.y])';

%%%
%
Resolution = 31;
[x1,x2] = meshgrid(linspace(0,LimX,Resolution),linspace(0,LimY,Resolution));

f1 = f1_fh(x1,x2);
f2 = f2_fh(x1,x2);
norm_f = sqrt(f1.^2 + f2.^2);
f1 = f1 ./ norm_f;
f2 = f2 ./ norm_f;


fig = figure(312);
delete(fig.Children);
Tl = tiledlayout(1,2,'TileSpacing','compact','Padding','compact');
ax = nexttile;
hold on, grid on, box on;

quiver(x1,x2,f1,f2,'DisplayName','Vector field')

for i = 1:size(Eq,2)
    plot(Eq(1,i),Eq(2,i),'.','MarkerSize',50)
end

term_event = @(t,x) hp_ode_terminal_event_ball(t,x,[0;0],norm([LimX,LimY]));
odeopts = odeset('Events',term_event);

Color_1 = [0 0.4470 0.7410];
Color_2 = [0.8500 0.3250 0.0980];
Color_3 = [0.9290 0.6940 0.1250];
Color_4 = [0.4940 0.1840 0.5560];
Color_5 = [0.4660 0.6740 0.1880];
Color_6 = [0.3010 0.7450 0.9330];
Color_7 = [0.6350 0.0780 0.1840];

xlim([0,LimX]);
ylim([0,LimY]);

[~,x] = ode45(f_fh,[0 1000],[K;0.1],odeopts);
Pl = plot(x(1,1),x(1,2),'.','Color',Color_3,'MarkerSize',25);
plot(x(:,1),x(:,2),'Color',Pl.Color,'LineWidth',1.5);

[~,x] = ode45(f_fh,[0 1000],Eq1+[1;1]*K/10,odeopts);
Pl = plot(x(1,1),x(1,2),'.','Color',Color_7,'MarkerSize',25);
plot(x(:,1),x(:,2),'Color',Pl.Color,'LineWidth',1.5);



ax = nexttile;
hold on, grid on, box on;

Tr = linspace(-3,3,101);
Pl = plot(Tr,Tr.^2/4);
for i = 1:size(Eq,2)
    J = J_fh(Eq(:,i));
    d = det(J);
    t = trace(J);

    plot(t,d,'.','MarkerSize',50)
end
XLim = ax.XLim;
xgrid = linspace(XLim(1),XLim(2),10000);
YLim = ax.YLim;
plot(XLim,[0 0],'Color',Pl.Color);
plot([0 0],YLim,'Color',Pl.Color);
plot(xgrid,xgrid.^2/4,'Color',Pl.Color);
ax.XLim = XLim;
ax.YLim = YLim;

