%% TASK 5: 3D plot and video 
% Nodes and solution from Task 5
r = nodes(1,:)';  % r-coordinate
z = nodes(2,:)';  % z-coordinate
u_sol = u;        % solution vector

% Create a fine phi grid for rotation
phi = linspace(0, 2*pi, 50);

% Create 3D coordinates
[Phi, R] = meshgrid(phi, r);
Z = repmat(z, 1, length(phi));
X = R .* cos(Phi);
Y = R .* sin(Phi);

% Interpolate u onto this grid (simple nearest neighbor)
U = repmat(u_sol, 1, length(phi));

% Plot
figure
surf(X, Y, Z, U, 'EdgeColor','none');
axis equal
xlabel('x'); ylabel('y'); zlabel('z');
title('Axisymmetric 3D solution');
colorbar
view(30,30)

% Assume you already have the 3D plot (X,Y,Z,U)
figure
h = surf(X,Y,Z,U, 'EdgeColor','none');
axis equal
colorbar
xlabel('x'); ylabel('y'); zlabel('z');
title('Axisymmetric 3D solution');

v = VideoWriter('mug_rotation.avi'); % create video file
open(v);

for angle = 0:2:360
    view(angle, 30); % rotate camera around z-axis
    drawnow;
    frame = getframe(gcf); % capture current frame
    writeVideo(v, frame);
end
close(v);


%% TASK 8: 9 snapshot of the solution over time: 

% --- Snapshot times (seconds)
snapshot_times = [0 50 100 150 200 250 300 350 400];
snapshot_steps = round(snapshot_times / tau) + 1;

Nsnap = length(snapshot_steps);
snapshots = zeros(Nnodes, Nsnap);

% Initial condition (t = 0)
snapshots(:,1) = u0;
snapshot_counter = 2;

% --- Time stepping parameters
T_end = snapshot_times(end);
Nt = ceil(T_end / tau);

% --- System matrix
A = (1/tau)*M + K;

% LU factorization (constant matrix)
[L,U,P,Q] = lu(A);

% --- Initialize solution
u = u0;

% --- Time stepping loop
for k = 1:Nt
    % Right-hand side
    b = (1/tau)*M*u + F;

    % Solve system
    u_next = Q * (U \ (L \ (P * b)));
    u = u_next;

    % Store snapshots at requested times
    if snapshot_counter <= Nsnap
        if k+1 == snapshot_steps(snapshot_counter)
            snapshots(:,snapshot_counter) = u;
            snapshot_counter = snapshot_counter + 1;
        end
    end
end

%% --- 3x3 grid plot (overview figure)
figure(20)
for i = 1:Nsnap
    subplot(3,3,i)
    pdeplot(model, 'XYData', snapshots(:,i), 'Mesh','off');
    axis equal
    title(['T = ', num2str(snapshot_times(i)), ' s'])
end
colorbar
sgtitle('Transient temperature evolution')

% Save combined figure
saveas(gcf, 'Task8_Transient_3x3.png');