% Linear flow on the torus: dtheta1/dt = 4, dtheta2/dt = 3
clear; close all; clc;

% Frequencies
omega1 = 3*sqrt(2);
omega2 = 3;

% Time discretization
Tmax = 60;
N = 30000;
t = linspace(0, Tmax, N);

% Initial conditions
theta1_0 = 0;
theta2_0 = 0;

% Solutions
theta1 = mod(theta1_0 + omega1*t, 2*pi);
theta2 = mod(theta2_0 + omega2*t, 2*pi);

% Torus parameters
R = 2;     % Major radius
r = 0.7;   % Minor radius

% Trajectory embedding
x = (R + r*cos(theta2)) .* cos(theta1);
y = (R + r*cos(theta2)) .* sin(theta1);
z = r * sin(theta2);

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

figure('units','normalized','outerposition',[0 0 1 1]);
hold on;

surf(Xt, Yt, Zt, ...
    'FaceColor', [0.6 0.6 0.9], ...
    'FaceAlpha', 0.25, ...
    'EdgeColor', 'none');

% Lighting
camlight headlight
lighting gouraud

% Plot full trajectory (static)
plot3(x, y, z, '.', ...
    'MarkerSize', 2, ...
    'Color', [0 0 1]);

% Initialize moving red dot
hDot = plot3(x(1), y(1), z(1), 'ro', ...
    'MarkerFaceColor', 'r', ...
    'MarkerSize', 8);

axis equal off
view(3)

% Animation loop
step = 20;   % Increase for faster animation
for k = 1:step:N
    set(hDot, 'XData', x(k), 'YData', y(k), 'ZData', z(k));
    drawnow
end