%% selected_ch11_14_figures_revised.m
% Figures for the report on Hirsch--Smale--Devaney Chapter 11--14.
% Run this file directly. It creates a folder named revised_report_figures
% in the current working directory and exports all figures as PNG files.
%
% Topics:
% 1. SIR and SIRS infectious disease models
% 2. van der Pol equation and Hopf bifurcation
% 3. Kepler problem and conic sections
% 4. Lorenz attractor and sensitive dependence

clear; close all; clc;
outDir = fullfile(pwd, 'revised_report_figures');
if ~exist(outDir, 'dir')
    mkdir(outDir);
end

%% Basic plotting settings
set(groot, 'defaultAxesFontSize', 10);
set(groot, 'defaultLineLineWidth', 1.8);
set(groot, 'defaultAxesLineWidth', 0.8);

%% ---------------- 1. SIR model ----------------
beta = 1.0;
nu = 0.35;
threshold = nu / beta;

sir = @(t,u) [-beta*u(1)*u(2); ...
               beta*u(1)*u(2) - nu*u(2); ...
               nu*u(2)];

% Figure 1: SIR time series
[t,u] = ode45(sir, [0 45], [0.92 0.06 0.02]);
S = u(:,1); I = u(:,2); R = u(:,3);
[~,idx] = max(I);
figure('Position',[100 100 760 430]);
plot(t,S,'DisplayName','S(t): susceptible'); hold on;
plot(t,I,'DisplayName','I(t): infected');
plot(t,R,'DisplayName','R(t): recovered');
xline(t(idx),'--','infection peak','LabelOrientation','horizontal');
yline(threshold,':','threshold \nu/\beta');
xlabel('time'); ylabel('population proportion');
title('SIR time series'); legend('Location','best'); grid on;
exportgraphics(gcf, fullfile(outDir,'sir_time_series_v2.png'), 'Resolution', 220);

% Figure 2: SIR phase portrait
figure('Position',[100 100 660 540]);
[Sg,Ig] = meshgrid(linspace(0.02,1.05,22), linspace(0.005,0.8,20));
U = -beta*Sg.*Ig;
V = beta*Sg.*Ig - nu*Ig;
N = sqrt(U.^2 + V.^2); N(N==0)=1;
quiver(Sg,Ig,U./N,V./N,0.45,'Color',[0.55 0.55 0.55]); hold on;
ics = [0.92 0.06 0.02; 0.75 0.20 0.05; 0.55 0.35 0.10; 0.30 0.45 0.25; 0.90 0.25 0.00];
for k = 1:size(ics,1)
    [~,uu] = ode45(sir, [0 65], ics(k,:));
    plot(uu(:,1), uu(:,2));
end
xline(threshold,'--','S=\nu/\beta','LineWidth',1.5);
Sline = linspace(0.05,1.1,300);
C = 0.92 + 0.06 - (nu/beta)*log(0.92);
Iline = C - Sline + (nu/beta)*log(Sline);
plot(Sline(Iline>0), Iline(Iline>0),'k-','LineWidth',1.0,'DisplayName','one level curve');
xlim([0 1.05]); ylim([0 0.72]); xlabel('S'); ylabel('I');
title('SIR phase portrait and threshold line'); grid on;
exportgraphics(gcf, fullfile(outDir,'sir_phase_v2.png'), 'Resolution', 220);

% Figure 3: outbreak size under different initial susceptible levels
figure('Position',[100 100 760 430]); hold on;
caseData = [0.25 0.05; 0.45 0.05; 0.92 0.05];
caseName = {'S_0 < \nu/\beta', 'S_0 slightly above', 'S_0 large'};
for k = 1:size(caseData,1)
    s0 = caseData(k,1); i0 = caseData(k,2); r0 = max(0,1-s0-i0);
    [tc,uc] = ode45(sir, [0 45], [s0 i0 r0]);
    plot(tc, uc(:,2), 'DisplayName', caseName{k});
end
xlabel('time'); ylabel('I(t)'); title('Different initial susceptible levels give different outbreaks');
legend('Location','best'); grid on;
exportgraphics(gcf, fullfile(outDir,'sir_threshold_cases_v2.png'), 'Resolution', 220);

%% SIRS model
betaS = 1; nuS = 1; muS = 1; tau = 2;
sirs = @(t,u) [-betaS*u(1)*u(2) + muS*(tau-u(1)-u(2)); ...
                betaS*u(1)*u(2) - nuS*u(2)];
figure('Position',[100 100 660 540]); hold on;
[Sg,Ig] = meshgrid(linspace(0,2,25), linspace(0,2,25));
mask = Sg + Ig <= tau + 1e-12;
U = -betaS*Sg.*Ig + muS*(tau-Sg-Ig);
V = betaS*Sg.*Ig - nuS*Ig;
N = sqrt(U.^2+V.^2); N(N==0)=1;
quiver(Sg(mask),Ig(mask),U(mask)./N(mask),V(mask)./N(mask),0.42,'Color',[0.55 0.55 0.55]);
ics = [1.8 0.05; 1.4 0.5; 0.4 1.0; 0.1 1.4; 1.0 0.1];
for k=1:size(ics,1)
    [~,uu] = ode45(sirs, [0 20], ics(k,:));
    plot(uu(:,1),uu(:,2));
end
Sline = linspace(0,2,300);
plot(Sline,tau-Sline,'k-','DisplayName','S+I=\tau');
xline(nuS/betaS,'--','S=\nu/\beta','LineWidth',1.4);
Inull = muS*(tau-Sline)./(betaS*Sline+muS);
plot(Sline,Inull,'LineWidth',2.0,'DisplayName','S-nullcline');
plot(tau,0,'ko','MarkerFaceColor','k');
plot(nuS/betaS, muS*(tau-nuS/betaS)/(nuS+muS), 'ko','MarkerFaceColor','k');
xlim([0 2.05]); ylim([0 2.05]); xlabel('S'); ylabel('I');
title('SIRS phase portrait in the invariant triangle'); legend('Location','northeast'); grid on;
exportgraphics(gcf, fullfile(outDir,'sirs_phase_v2.png'), 'Resolution', 220);

%% ---------------- 2. van der Pol and Hopf ----------------
vdp = @(mu) @(t,u) [u(2)-u(1)^3+mu*u(1); -u(1)];

% Figure 5: nullclines and direction field
figure('Position',[100 100 660 540]);
[xg,yg] = meshgrid(linspace(-2.2,2.2,25), linspace(-2.7,2.7,25));
U = yg - xg.^3 + xg; V = -xg;
N = sqrt(U.^2+V.^2); N(N==0)=1;
quiver(xg,yg,U./N,V./N,0.45,'Color',[0.55 0.55 0.55]); hold on;
xlineData = linspace(-2.1,2.1,400);
plot(xlineData, xlineData.^3-xlineData,'LineWidth',2.2,'DisplayName','x-nullcline: y=x^3-x');
xline(0,'--','y-nullcline: x=0','LineWidth',1.5);
plot(0,0,'ko','MarkerFaceColor','k');
xlim([-2.2 2.2]); ylim([-2.7 2.7]); xlabel('x'); ylabel('y');
title('van der Pol nullclines and direction field'); legend('Location','northwest'); grid on;
exportgraphics(gcf, fullfile(outDir,'vdp_nullclines_v2.png'), 'Resolution', 220);

% Figure 6: attracting limit cycle
figure('Position',[100 100 660 540]); hold on;
ics = [0.1 0.1; 1.8 1.8; -2.2 1.2; 2.0 -1.7; -0.5 2.2];
for k=1:size(ics,1)
    [~,uu] = ode45(vdp(1), [0 45], ics(k,:));
    plot(uu(:,1),uu(:,2));
end
[tLC,uLC] = ode45(vdp(1), [0 80], [2.5 0]);
startIdx = max(1, length(tLC)-1300);
plot(uLC(startIdx:end,1), uLC(startIdx:end,2), 'k-', 'LineWidth',2.1, 'DisplayName','attracting cycle');
xlabel('x'); ylabel('y'); title('Attracting limit cycle of the van der Pol equation');
xlim([-2.4 2.4]); ylim([-3.0 3.0]); legend('Location','best'); grid on;
exportgraphics(gcf, fullfile(outDir,'vdp_limit_cycle_v2.png'), 'Resolution', 220);

% Figure 7: energy balance
figure('Position',[100 100 660 430]);
x = linspace(-2,2,500); Wdot = x.^2.*(1-x.^2);
plot(x,Wdot,'LineWidth',2.2); hold on; yline(0,'k-'); xline(-1,'--'); xline(1,'--');
xlabel('x'); ylabel('dW/dt = x^2(1-x^2)'); title('Energy balance in the van der Pol equation'); grid on;
exportgraphics(gcf, fullfile(outDir,'vdp_energy_balance_v2.png'), 'Resolution', 220);

% Figure 8: Hopf phase portraits
figure('Position',[100 100 850 720]);
muList = [-0.5 -0.1 0.05 0.2];
for j=1:4
    subplot(2,2,j); hold on;
    mu = muList(j);
    ics = [0.2 0.1; 1.6 0.6; -1.6 1.0; 0.5 -1.8];
    for k=1:size(ics,1)
        [~,uu] = ode45(vdp(mu), [0 55], ics(k,:));
        plot(uu(:,1), uu(:,2));
    end
    plot(0,0,'ko','MarkerFaceColor','k');
    title(['\mu = ', num2str(mu)]); xlabel('x'); ylabel('y');
    xlim([-2 2]); ylim([-2.4 2.4]); grid on;
end
exportgraphics(gcf, fullfile(outDir,'hopf_phase_portraits_v2.png'), 'Resolution', 220);

% Figure 9: Hopf bifurcation diagram, numerical amplitude
muGrid = linspace(-0.5,0.8,70); amp = zeros(size(muGrid));
for j=1:length(muGrid)
    mu = muGrid(j);
    if mu > 0
        [~,uu] = ode45(vdp(mu), [0 180], [1.5 0.2]);
        rad = sqrt(uu(:,1).^2 + uu(:,2).^2);
        tail = max(1,length(rad)-800):length(rad);
        amp(j) = max(rad(tail));
    end
end
figure('Position',[100 100 660 430]);
plot(muGrid, amp, 'LineWidth',2.1); hold on; xline(0,'--','bifurcation value');
xlabel('parameter \mu'); ylabel('approximate radius');
title('Hopf bifurcation: birth of a small attracting cycle'); grid on;
exportgraphics(gcf, fullfile(outDir,'hopf_bifurcation_diagram_v2.png'), 'Resolution', 220);

%% ---------------- 3. Kepler problem ----------------
% Figure 10: equal-area law
e = 0.55; p = 1.0;
th = linspace(0,2*pi,1000); rr = p./(1+e*cos(th));
xx = rr.*cos(th); yy = rr.*sin(th);
figure('Position',[100 100 660 540]); hold on;
plot(xx,yy,'LineWidth',2.2,'DisplayName','elliptic orbit');
plot(0,0,'o','MarkerSize',9,'MarkerFaceColor',[1 0.65 0.1],'MarkerEdgeColor','k','DisplayName','force center/focus');
sectors = [0.35 0.72; 2.55 2.92];
for j=1:2
    tt = linspace(sectors(j,1), sectors(j,2), 100); rr = p./(1+e*cos(tt));
    fill([0 rr.*cos(tt) 0], [0 rr.*sin(tt) 0], [0.4 0.7 1.0], 'FaceAlpha',0.28, 'EdgeColor','none');
    plot([0 rr(1)*cos(tt(1))],[0 rr(1)*sin(tt(1))],'Color',[0.35 0.35 0.35]);
    plot([0 rr(end)*cos(tt(end))],[0 rr(end)*sin(tt(end))],'Color',[0.35 0.35 0.35]);
end
axis equal; xlabel('x'); ylabel('y'); title('Equal-area law from angular momentum conservation'); legend('Location','northeast'); grid on;
exportgraphics(gcf, fullfile(outDir,'kepler_equal_area_v2.png'), 'Resolution', 220);

% Figure 11: effective potential
ell = 1; rgrid = linspace(0.18,6,1000);
Veff = ell^2./(2*rgrid.^2) - 1./rgrid;
figure('Position',[100 100 720 440]);
plot(rgrid,Veff,'LineWidth',2.3,'DisplayName','V_{eff}(r) = ell^2/(2r^2)-1/r'); hold on;
yline(-0.28,'--','E<0: bounded'); yline(0,'--','E=0: parabolic boundary'); yline(0.16,'--','E>0: unbounded');
xlim([0.18 6]); ylim([-0.55 0.5]); xlabel('r'); ylabel('energy');
title('Effective potential for the Newtonian central force problem'); legend('Location','best'); grid on;
exportgraphics(gcf, fullfile(outDir,'kepler_effective_potential_v2.png'), 'Resolution', 220);

% Figure 12: conic sections
figure('Position',[100 100 720 540]); hold on;
eList = [0.55 1.0 1.35]; labels = {'E<0: ellipse','E=0: parabola','E>0: hyperbola'};
for j=1:3
    e = eList(j); p = 1;
    if e < 1
        th = linspace(0,2*pi,1000);
    elseif abs(e-1) < 1e-12
        th = linspace(-2.25,2.25,800);
    else
        lim = acos(-1/e)-0.08; th = linspace(-lim,lim,900);
    end
    rr = p./(1+e*cos(th));
    plot(rr.*cos(th), rr.*sin(th), 'DisplayName', labels{j}, 'LineWidth',2.2);
end
plot(0,0,'o','MarkerSize',8,'MarkerFaceColor',[1 0.65 0.1],'MarkerEdgeColor','k','DisplayName','focus');
axis equal; xlim([-3.8 3.8]); ylim([-2.8 2.8]); xlabel('x'); ylabel('y');
title('Kepler orbits as conic sections'); legend('Location','best'); grid on;
exportgraphics(gcf, fullfile(outDir,'kepler_conics_v2.png'), 'Resolution', 220);

% Figure 13: logical chain
figure('Position',[100 100 760 280]); axis off;
boxText = {'Newton law','central force','angular momentum','W=1/r','conic orbit'};
xpos = [0.08 0.28 0.50 0.70 0.90]; ypos = 0.55;
for j=1:length(boxText)
    annotation('textbox',[xpos(j)-0.07 ypos-0.08 0.14 0.16], 'String',boxText{j}, 'HorizontalAlignment','center', 'VerticalAlignment','middle', 'FitBoxToText','on', 'BackgroundColor','white');
end
for j=1:length(boxText)-1
    annotation('arrow',[xpos(j)+0.07 xpos(j+1)-0.07],[ypos ypos]);
end
title('Logical chain of the Kepler derivation');
exportgraphics(gcf, fullfile(outDir,'kepler_logic_chain_v2.png'), 'Resolution', 220);

%% ---------------- 4. Lorenz system ----------------
sigma = 10; b = 8/3; rho = 28;
lorenz = @(t,u) [sigma*(u(2)-u(1)); rho*u(1)-u(2)-u(1)*u(3); u(1)*u(2)-b*u(3)];

% Figure 14: Lorenz attractor
[tL,uL] = ode45(lorenz, [0 55], [0 2 0]);
figure('Position',[100 100 720 540]);
plot3(uL(800:end,1), uL(800:end,2), uL(800:end,3), 'LineWidth',0.7);
xlabel('x'); ylabel('y'); zlabel('z'); title('Lorenz attractor, sigma=10, b=8/3, r=28'); grid on; view(-58,23);
exportgraphics(gcf, fullfile(outDir,'lorenz_attractor_v2.png'), 'Resolution', 220);

% Figure 15: coordinate time series
figure('Position',[100 100 760 430]);
idx = (tL>=10 & tL<=30);
plot(tL(idx),uL(idx,1),'DisplayName','x(t)'); hold on;
plot(tL(idx),uL(idx,2),'DisplayName','y(t)');
plot(tL(idx),uL(idx,3),'DisplayName','z(t)');
xlabel('time'); ylabel('coordinates'); title('Lorenz coordinate time series'); legend('Location','best'); grid on;
exportgraphics(gcf, fullfile(outDir,'lorenz_time_series_v2.png'), 'Resolution', 220);

% Figure 16 and 17: sensitive dependence
[tA,uA] = ode45(lorenz, [0 35], [0 2.00 0]);
[tB,uB] = ode45(lorenz, [0 35], [0 2.01 0]);
tCommon = linspace(0,35,5000);
uAi = interp1(tA,uA,tCommon); uBi = interp1(tB,uB,tCommon);
dist = sqrt(sum((uAi-uBi).^2,2));
figure('Position',[100 100 760 430]);
semilogy(tCommon, dist, 'LineWidth',2.1); yline(1,'--');
xlabel('time'); ylabel('distance between two solutions'); title('Sensitive dependence in the Lorenz system'); grid on;
exportgraphics(gcf, fullfile(outDir,'lorenz_sensitivity_distance_v2.png'), 'Resolution', 220);
figure('Position',[100 100 760 430]);
plot(tCommon,uAi(:,1),'DisplayName','x(t), initial y=2.00'); hold on;
plot(tCommon,uBi(:,1),'--','DisplayName','x(t), initial y=2.01');
xlim([0 28]); xlabel('time'); ylabel('x-coordinate');
title('Two close Lorenz initial conditions eventually separate'); legend('Location','best'); grid on;
exportgraphics(gcf, fullfile(outDir,'lorenz_sensitivity_traces_v2.png'), 'Resolution', 220);

% Figure 18: numerical Poincare section x = 0, positive crossing
% Use a longer integration for a denser section plot.
[tSec,uSec] = ode45(lorenz, [0 180], [0 2 0]);
ys = []; zs = [];
for k=1:length(tSec)-1
    dxdt = sigma*(uSec(k,2)-uSec(k,1));
    if uSec(k,1) < 0 && uSec(k+1,1) >= 0 && dxdt > 0
        a = (0-uSec(k,1))/(uSec(k+1,1)-uSec(k,1));
        ys(end+1) = uSec(k,2) + a*(uSec(k+1,2)-uSec(k,2)); %#ok<SAGROW>
        zs(end+1) = uSec(k,3) + a*(uSec(k+1,3)-uSec(k,3)); %#ok<SAGROW>
    end
end
figure('Position',[100 100 660 540]);
scatter(ys,zs,14,'filled','MarkerFaceAlpha',0.55);
xlabel('y on section x=0'); ylabel('z on section x=0');
title('Numerical Poincare section: x = 0, positive crossing'); grid on;
exportgraphics(gcf, fullfile(outDir,'lorenz_poincare_section_v2.png'), 'Resolution', 220);

% Figure 19: Lorenz equilibria versus r
rGrid = linspace(0,35,500); amp = sqrt(max(b*(rGrid-1),0));
figure('Position',[100 100 720 430]);
plot(rGrid, zeros(size(rGrid)), 'DisplayName','origin equilibrium branch'); hold on;
plot(rGrid, amp, 'DisplayName','Q_+: x=y=+sqrt(b(r-1))');
plot(rGrid, -amp, 'DisplayName','Q_-: x=y=-sqrt(b(r-1))');
xline(1,'--','r=1'); xlabel('r'); ylabel('x-coordinate of equilibria');
title('Lorenz equilibria and the pitchfork-like bifurcation at r=1'); legend('Location','best'); grid on;
exportgraphics(gcf, fullfile(outDir,'lorenz_equilibria_bifurcation_v2.png'), 'Resolution', 220);

fprintf('Done. Figures saved in: %s\n', outDir);