clc;
clear variables;

%-------Parametros de entrada----------------------------------------------
Lambda1 = 0.001 ; % (Lambda) (W/m degree C)
Lambda2 = 0.06,
Cp1 = 1; % (J/kg degree C)
Cp2 = 1;

%-------Parametros geometricos---------------------------------------------
Hx = 1; % Largo de la barra
Hx1 = 0.5
Nx = 5; % Numero de nodos
dX = Hx/(Nx-2); % Delta x
Hy = 0.5; % Largo de barra
Ny = 5; % Numero de nodos en y
% Manten Hx = 1 y Nx = 5
%-------Coordenadas del sistema 1D-----------------------------------------
% Uso linspace
x = linspace (0, Hx, Nx);
y = linspace (0, Hx, Ny);

x(1) = 0; %Frontera oeste
x(Nx) = Hx; % Frontera este

% Nodos internos
for i = 2:Nx-1
    x(i) = ((i-2)*dX) + dX/2;
end
[X, Y] = meshgrid(x, y);

% Malla
disp('Coordenadas de los nodos (x_i):');
disp(x)

figure(1);
plot(X,Y, '-ko');
xlabel('x (m)');
ylabel('y (m)');
title('malla 2D');
axis equal; grid on;

% Matrices de propiedades
Lamba = zeros(Nx, Ny);
Cp = zeros(Nx, Ny);
for i = 1:Nx
for j = 1:Ny
if x(i) <= Hx1
Lambda(i,j) = Lambda2;
else
Lambda(i,j) = Lambda1;
end
end
end

Gamma = Lambda ./ Cp1
% coeficientes
a_N=zeros(Nx,Ny);
a_S=zeros(Nx,Ny);
a_W=zeros(Nx,Ny);
a_E=zeros(Nx,Ny);
a_P=zeros(Nx,Ny);
b=zeros(Nx,Ny);
% Codiciones de 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;
% Termino fuente 
for i =1:Nx
 b(i,Ny)=100;
end

for i=1:Nx
 b(i,1)=600;
end
% Calcular coeficientes para nodos internos
for i =2:Nx-1
 for j =2:Ny-1
   delta_X_WP = X(i) - X(i-1);
   delta_X_PE = X(i+1) - X(i);
   delta_y_SP = y(j) - y(j-1);
   delta_y_PN = y(j+1) - y (j);

   Gamma_P = Gamma(i,j);
   Gamma_W = Gamma(i-1,j);
   Gamma_E = Gamma(i+1j);
   Gamma_S = Gamma(ij-1);
   Gamma_N = Gamma(ij+1);

   Gamma_w = (Gamma_W* Gamma_P* delta_x_WP)/((Gamma_P* delta_x_WP/2) + (Gamma_W* delta_x_WP/2));
   Gamma_e = (Gamma_P* Gamma_E* delta_x_PE)/((Gamma_E* delta_x_PE/2) + (Gamma_P* delta_x_SP/2));
   Gamma_s = (Gamma_S* Gamma_P* delta_y_SP)/((Gamma_P* delta_y_SP/2) + (Gamma_S* delta_y_WP/2));
   Gamma_n = (Gamma_P* Gamma_N* delta_y_PN)/((Gamma_N* delta_y_PN/2) + (Gamma_P* delta_y_PN/2));

   a_W
   
end
end