clc;
clear variables;

%-------Parametros de entrada----------------------------------------------
Lambda1 = 0; % λ (W/m°C) 
Cp1 = 130; % (J/kg°C)
Gamma = Lambda/Cp; % Γ

%-------Parametros geometricos---------------------------------------------
Hx = 1; % Largo de la barra
Nx = 7; % Numero de nodos
dx = Hx/(Nx-2); % Δx
dy=Hy/(Nx-2);
Hm=0.5;
%-------Coordenadas del sistema 1D-----------------------------------------
x(1) = 0; x(Nx) = Hx; % Fronteras
y(1)=0;y(Ny)=Hy;

% Nodos internos
for i = 2:Nx-1
    x(i) = ((i-2)*dx + dx/2);
end
for j=2:Ny-1
y(j)=(cj-2)*dy+dy/2);
end
% Malla
[X,Y]= meshgrid(x,y);
disp(x)
% propiedades de materiales
Lambda=zeros(Nx,Ny);
Cp=Zeros(Nx,Ny);
for i=1;Nx
for j=1:Ny
if x(i)<= Hm
lambda(i,j)=Lambda1;
Cp(i,j)=Cp1;
else
Lambda(i,j)=Lambda2;
Cp(i,j)=Cp2;
end
end
end
%calcular Gamma
Gamma=Lambda./Cp;

%-------Coeficientes-------------------------------------------------------
a_N = Zeros (Nx,Ny);% coeficiente norte
a_S= zeros (Nx,Ny);%  coeficiente sur
a_W = zeros(Nx, y); % Coeficiente oeste
a_E = zeros(Nx, y); % Coeficiente este
a_P = zeros(Nx, y); % Coeficiente central (Centroide)
b = zeros(Nx, y);   % Término fuente

% Coeficientes para nodos internos
for i = 2:Nx-1
    delta_x_WP = x(i) - x(i-1); % Distancia entre nodo W y P, δWP
    delta_x_PE = x(i+1) - x(i); % Distancia entre nodo P y E, δPE
    
    a_W(i) = Gamma / delta_x_WP;
    a_E(i) = Gamma / delta_x_PE;
    a_P(i) = a_W(i) + a_E(i);
end

%calcular coeficientes para nodos internos 
for i=2Nx-1
 for j=2Ny-1
delta_x_WP=x(i)-x(i-1);
delta_x_PE=x(i-1)-x(i);
delta_y_SP=y(i)-y(i-1);
delta_y_PN=y(i+1)-y(i);

Gamma_P=Gamma(i,j);

% Coeficientes para nodos frontera
a_P(1,:) = 1; 
a_E(1,:) = 0; 
a_S(1,:) = 0; 
a_N(1,:) = 0;
a_W(1;:) = 0;

a_P(Nx,:)=1;
a_E(Nx,:)=0;
a_S(Nx,:)=0;
a_W(Nx,:)=0;
a_N(Nx,:)=0;

a_p(:,1)=1;
a_E(:,1)=0;
a_S(:,1)=0;
a_W(:,1)=0;
a_N(:,1)=0;

a_p(:,Ny)=1;
a_E(:,Ny)=0;
a_S(:,Ny)=0;
a_W(:,Ny)=0;
a_N(:,Ny)=0;

b(1) = 0; % Condición de frontera en x(1)

a_P(Nx) = 1; 
a_W(Nx) = 0; 
a_E(Nx) = 0; 

b(Nx) = 100; % Condición de frontera en x(Nx)

% Coeficientes
disp('Coeficientes a_W:');
disp(a_W);
disp('Coeficientes a_E:');
disp(a_E);
disp('Coeficientes a_P:');
disp(a_P);
disp('Términos fuente b:');
disp(b);

%-------Inicialización de temperaturas-------------------------------------
T = zeros(Nx, 1);   
T(1) = 0;            
T(Nx) = 100;         

% Asignar valor inicial nodos internos
for i = 2:Nx-1
    T(i) = 50;
end

% Mostrar temperatura inicial
disp('Temperatura inicial T(x):');
disp(T);

%-------Método de Jacobi--------------------------------------------------
max_iter = 1000;        % Número máximo de iteraciones
tolerancia = 1e-20;    % Tolerancia para convergencia

for iter = 1:max_iter
    % Guardar temperaturas de la iteración anterior
    T_old = T;
    
    % Actualizar temperaturas en nodos internos
    for i = 2:Nx-1
        T(i) = (a_W(i) * T_old(i-1) + a_E(i) * T_old(i+1) + b(i)) / a_P(i);
    end
    
    disp(['Iteración: ', num2str(iter)]);
    disp(T);
    
    % Verificar convergencia
    if max(abs(T - T_old)) < tolerancia
        disp(['Convergencia alcanzada en la iteración: ', num2str(iter)]);
        break;
    end
end

% Mostrar resultados
% disp('Distribución de temperatura T(x):');
% disp(T);

%-------Gráfica de resultados---------------------------------------------
figure;
plot(x, T, 'b-o', 'LineWidth', 2, 'DisplayName', 'Solución numérica');
xlabel('Posición (x)');
ylabel('Temperatura (°C)');
title('Distribución de temperatura en la barra');
legend;
grid on;