% Parámetros del problema
L = 1000;      % Longitud del océano (en metros)
alpha = 1.4e-6; % Difusividad térmica del agua (m^2/s)
h = 10;       % Coeficiente de transferencia de calor
k = 0.6;      % Conductividad térmica del agua
u_amb = 20;   % Temperatura ambiente (grados Celsius)
u0 = 5;       % Temperatura inicial (grados Celsius)

% Dominio espacial y temporal
Nx = 50; % Número de puntos en x
Nt = 100; % Número de pasos en t
x = linspace(0, L, Nx);
t = linspace(0, 2000, Nt); % Tiempo hasta 1 día (en segundos)

% Número de términos en la serie
N_terms = 5;

% Solución de la serie de Fourier
u = zeros(Nx, Nt);
for n = 1:N_terms
    beta_n = (2*n - 1) * pi / (2 * L); % Modificación para condición de Neumann en x=L
    lambda_n = alpha * beta_n^2; % Parámetro de la ecuación
    A_n = 2 * (u_amb-u0) * sin(beta_n * L) / (beta_n * L); % 
    for i = 1:Nx
        for j = 1:Nt
            u(i, j) = u(i, j) + A_n * cos(beta_n * x(i)) * exp(-lambda_n * t(j));
        end
    end
end

% Solución
u = u_amb - u;

% Graficar evolución de la temperatura
figure;
[X, T] = meshgrid(x, t/3600); % Convertir tiempo a horas
surf(X, T, u', 'EdgeColor', 'none');
xlabel('Distancia (m)');
ylabel('Tiempo (horas)');
zlabel('Temperatura (°C)');
title('Evolución de la temperatura en el océano');
colorbar;
view(3);

% Parámetros del problema
L = 1000;      % Longitud del océano (en metros)
alpha = 1.4e-6; % Difusividad térmica del agua (m^2/s)
h = 10;       % Coeficiente de transferencia de calor
k = 0.6;      % Conductividad térmica del agua
u_amb = 20;   % Temperatura ambiente (grados Celsius)
u0 = 5;       % Temperatura inicial (grados Celsius)

% Dominio espacial y temporal
Nx = 50; % Número de puntos en x
Nt = 50; % Número de pasos en t (mucho menor para avanzar rápido)
t_max = 10000000 * 86400; % Simulación hasta 10 días en segundos
t = linspace(0, t_max, Nt); % Tiempo con pasos mucho más grandes
x = linspace(0, L, Nx);

% Número de términos en la serie
N_terms = 10;

% Solución de la serie de Fourier
u = zeros(Nx, Nt);
for n = 1:N_terms
    beta_n = (2*n - 1) * pi / (2 * L); % Modificación para condición de Neumann en x=L
    lambda_n = alpha * beta_n^2; % Parámetro de la ecuación
    A_n = 2 * (u_amb - u0) * sin(beta_n * L) / (beta_n * L);
    for i = 1:Nx
        for j = 1:Nt
            u(i, j) = u(i, j) + A_n * cos(beta_n * x(i)) * exp(-lambda_n * t(j));
        end
    end
end

% Solución
u = u_amb - u;

% Crear el objeto VideoWriter para guardar el video
v = VideoWriter('evolucion_calor_océano.mp4', 'MPEG-4'); % Nombre del archivo y formato
v.FrameRate = 5; % Velocidad de los cuadros (puedes ajustarlo según lo necesites)
open(v); % Abrir el objeto para grabar el video

% Crear la figura
figure;
colormap hot; % Mapa de calor para mejor visualización

% Crear la animación y capturar cada cuadro
for j = 1:Nt
    imagesc(x, [0, L], repmat(u(:, j)', [Nx, 1])); % Representación del perfil de temperatura
    set(gca, 'YDir', 'normal'); % Mantener el eje y en la dirección correcta
    xlabel('Distancia (m)');
    ylabel('Profundidad (m)');
    title(['Evolución del calor en el océano - t = ', num2str(t(j)/86400, '%.1f'), ' días']);
    colorbar;
    
    % Ajustar la escala de colores dinámicamente
    caxis([min(u(:)), max(u(:))]); 
    
    pause(0.2); % Control de la velocidad de la animación
    
    % Capturar el cuadro actual y agregarlo al video
    frame = getframe(gcf); % Capturar el cuadro de la figura
    writeVideo(v, frame); % Escribir el cuadro en el archivo de video
end

% Cerrar el archivo de video
close(v); % Finalizar y guardar el archivo de video