% Coupled phase oscillators on the torus
clear; close all; clc;

%% Parameters
omega1 = 1;
omega2 = 1.3;
k1 = .4;
k2 = .4;

Tmax = 120;
dt = 0.01;
t = 0:dt:Tmax;
N = length(t);

% Initial conditions
theta1 = zeros(1,N);
theta2 = zeros(1,N);
theta1(1) = 0.3;
theta2(1) = 2.0;

%% Time integration (Euler)
for n = 1:N-1
    theta1(n+1) = theta1(n) + dt*( ...
        omega1 + k1*sin(theta2(n) - theta1(n)) );

    theta2(n+1) = theta2(n) + dt*( ...
        omega2 + k2*sin(theta1(n) - theta2(n)) );
end

% Reduce mod 2*pi for all plots
theta1_mod = mod(theta1, 2*pi);
theta2_mod = mod(theta2, 2*pi);

%% Break lines at jumps for unit square plot
jump_threshold = pi;  % if jump > pi, insert NaN
theta1_plot = theta1_mod;
theta2_plot = theta2_mod;

for n = 2:N
    if abs(theta1_plot(n) - theta1_plot(n-1)) > jump_threshold || ...
       abs(theta2_plot(n) - theta2_plot(n-1)) > jump_threshold
        theta1_plot(n-1) = NaN;
        theta2_plot(n-1) = NaN;
    end
end

%% Torus embedding parameters
R = 2;
r = 0.6;

x = (R + r*cos(theta2_mod)) .* cos(theta1_mod);
y = (R + r*cos(theta2_mod)) .* sin(theta1_mod);
z = r * sin(theta2_mod);

%% Torus surface
[u,v] = meshgrid(linspace(0,2*pi,60));
Xt = (R + r*cos(v)) .* cos(u);
Yt = (R + r*cos(v)) .* sin(u);
Zt = r * sin(v);

%% Figure setup
figure('units','normalized','outerposition',[0 0 1 1]);

% --- (1) Unit square representation ---
subplot(1,3,1); hold on;

% Growing trajectory on unit square with jumps handled
hSqPath = plot(theta1_plot(1), theta2_plot(1), ...
    'Color',[0.2 0.2 0.8], 'LineWidth',1.5);

% Moving point
hSq = plot(theta1_plot(1), theta2_plot(1), ...
    'ro', 'MarkerFaceColor','r');

xlim([0 2*pi]); ylim([0 2*pi]);
axis square
xlabel('\theta_1'); ylabel('\theta_2');

% --- (2) Embedded torus ---
subplot(1,3,2); hold on;

surf(Xt, Yt, Zt, ...
    'FaceAlpha',0.18, ...
    'EdgeColor','none', ...
    'FaceColor',[0.6 0.7 1]);

hPath = plot3(x(1), y(1), z(1), ...
    'Color',[0.2 0.2 0.8], ...
    'LineWidth',1.5);

hTorus = plot3(x(1), y(1), z(1), ...
    'ro','MarkerFaceColor','r');

axis equal
axis tight
view(0,75)
camlight headlight
lighting gouraud

% --- (3) Two phases on the circle ---
subplot(1,3,3); hold on;
th = linspace(0,2*pi,300);
plot(cos(th), sin(th), 'k','LineWidth',1);

hC1 = plot(cos(theta1_mod(1)), sin(theta1_mod(1)), ...
    'ro','MarkerFaceColor','r','MarkerSize',8);
hC2 = plot(cos(theta2_mod(1)), sin(theta2_mod(1)), ...
    'bo','MarkerFaceColor','b','MarkerSize',8);
axis square

%% Animation
step = 10;
for n = 2:step:N
    % Unit square: growing path + moving point
    set(hSqPath, 'XData', theta1_plot(1:n), 'YData', theta2_plot(1:n));
    set(hSq, 'XData', theta1_plot(n), 'YData', theta2_plot(n));

    % Torus
    set(hPath, 'XData', x(1:n), 'YData', y(1:n), 'ZData', z(1:n));
    set(hTorus, 'XData', x(n), 'YData', y(n), 'ZData', z(n));

    % Circle points
    set(hC1, 'XData', cos(theta1_mod(n)), 'YData', sin(theta1_mod(n)));
    set(hC2, 'XData', cos(theta2_mod(n)), 'YData', sin(theta2_mod(n)));

    drawnow
end