%Gerardo Hernandez Duenas
function Calor2D

%Vamos a resolver la ecuacion del calor en 2D
%u_t = \nabla \cdot (kappa \grad u)

%% Parameters

Nx = 100; Ny = 100;
Lx = 10; Ly = 1;
dx = Lx/Nx; dy = Ly/Ny;

T = 10;
Movie = 0;

%% Definir la conductividad

kappa_jmhalf = zeros(Nx+2,Ny+2); kappa_jphalf = zeros(Nx+2,Ny+2); 
kappa_kmhalf = zeros(Nx+2,Ny+2); kappa_kphalf = zeros(Nx+2,Ny+2);
m = 2; rx = Lx/8; ry = Ly/8;
for k = 0:Ny+1
    k_shift = k+1;
    y_k = dy*(k-1/2);
    y_kmhalf = dy*(k-1); y_kphalf = dy*(k);
    for j = 0:Nx+1
        j_shift = j+1;
        x_j = dx*(j-1/2);
        x_jmhalf = dx*(j-1); x_jphalf = dx*(j-1);
        kappa_jmhalf(j_shift,k_shift) = 1;
        if ((x_jmhalf - Lx/4)/rx)^2 + ((y_k-Ly/4)/ry)^2 < 1
            kappa_jmhalf(j_shift,k_shift) = 0;
        end
        if ((x_jmhalf - Lx/4)/rx)^2 + ((y_k-3*Ly/4)/ry)^2 < 1
            kappa_jmhalf(j_shift,k_shift) = 0;
        end
        if ((x_jmhalf - Lx/2)/rx)^2 + ((y_k-Ly/2)/ry)^2 < 1
            kappa_jmhalf(j_shift,k_shift) = 0;
        end       
        
        kappa_jphalf(j_shift,k_shift) = 1;
        if ((x_jphalf - Lx/4)/rx)^2 + ((y_k-Ly/4)/ry)^2 < 1
            kappa_jphalf(j_shift,k_shift) = 0;
        end
        if ((x_jphalf - Lx/4)/rx)^2 + ((y_k-3*Ly/4)/ry)^2 < 1
            kappa_jphalf(j_shift,k_shift) = 0;
        end
        if ((x_jphalf - Lx/2)/rx)^2 + ((y_k-Ly/2)/ry)^2 < 1
            kappa_jphalf(j_shift,k_shift) = 0;
        end
        
        kappa_kmhalf(j_shift,k_shift) = 1;
        if ((x_j - Lx/4)/rx)^2 + ((y_kmhalf-Ly/4)/ry)^2 < 1
            kappa_kmhalf(j_shift,k_shift) = 0;
        end
        if ((x_j - Lx/4)/rx)^2 + ((y_kmhalf-3*Ly/4)/ry)^2 < 1
            kappa_kmhalf(j_shift,k_shift) = 0;
        end
        if ((x_j - Lx/2)/rx)^2 + ((y_kmhalf-Ly/2)/ry)^2 < 1
            kappa_kmhalf(j_shift,k_shift) = 0;
        end
        
       kappa_kphalf(j_shift,k_shift) = 1;
        if ((x_j - Lx/4)/rx)^2 + ((y_kphalf-Ly/4)/ry)^2 < 1
            kappa_kphalf(j_shift,k_shift) = 0;
        end
        if ((x_j - Lx/4)/rx)^2 + ((y_kphalf-3*Ly/4)/ry)^2 < 1
            kappa_kphalf(j_shift,k_shift) = 0;
        end
        if ((x_j - Lx/2)/rx)^2 + ((y_kphalf-Ly/2)/ry)^2 < 1
            kappa_kphalf(j_shift,k_shift) = 0;
        end    
    end
end

%% Definir las condiciones iniciales

u = zeros(Nx+2,Ny+2);

x_axis= zeros(Nx+2,1); y_axis= zeros(Ny+2,1);
for k = 0:Ny+1
    k_shift = k+1;
    yk = dy*(k-1/2);
    y_axis(k_shift,1) = yk;
    for j = 0:Nx+1
        j_shift = j+1;
        xj = dx*(j-1/2);
        x_axis(j_shift,1) = xj;
        u(j_shift,k_shift) = 0; %(1-cos(2*pi*(yk-Ly/2)))*(1-sin(2*pi*(xj-Ly/2)));  %1+xj+3*yk; %   
    end
end

% %[x y] = meshgrid(x_axis,y_axis);
% figure(1)
% hold on
% contour(x_axis,y_axis,transpose(kappa_kphalf),50); colorbar;
% %plot3(x,y,kappa_jmhalf,'r *'); 
% stop


%% Numerical scheme

% Condicion de estabilidad:

CFL = 0.45;

max_kappa = max(max(kappa_jmhalf));
max_kappa = max(max_kappa,max(max(kappa_jmhalf)));
max_kappa = max(max_kappa,max(max(kappa_kmhalf)));
max_kappa = max(max_kappa,max(max(kappa_kphalf)));

dt = CFL/max_kappa/(1/dx^2+1/dy^2);

u_n = u;
u_np1 = u;

nu_jmhalf = dt*kappa_jmhalf/dx^2;
nu_jphalf = dt*kappa_jphalf/dx^2;
nu_kmhalf = dt*kappa_kmhalf/dy^2;
nu_kphalf = dt*kappa_kphalf/dy^2;

[x y] = meshgrid(x_axis,y_axis);
%plot3(x,y,kappa_jmhalf); 
% contour(x_axis,y_axis,transpose(u_np1),50); colorbar;
% stop

t=0;
while t < T
    for k = 1:Ny
        k_shift = k+1;
        for j = 1:Nx
            j_shift = j+1;
            
            u_np1(j_shift,k_shift) = u_n(j_shift,k_shift);
            u_np1(j_shift,k_shift) = u_np1(j_shift,k_shift) + nu_jphalf(j_shift,k_shift)*u_n(j_shift+1,k_shift); 
            u_np1(j_shift,k_shift) = u_np1(j_shift,k_shift) - (nu_jphalf(j_shift,k_shift)+nu_jmhalf(j_shift,k_shift))*u_n(j_shift,k_shift);
            u_np1(j_shift,k_shift) = u_np1(j_shift,k_shift) + nu_jmhalf(j_shift,k_shift)*u_n(j_shift-1,k_shift);
            u_np1(j_shift,k_shift) = u_np1(j_shift,k_shift) + nu_kphalf(j_shift,k_shift)*u_n(j_shift,k_shift+1);
            u_np1(j_shift,k_shift) = u_np1(j_shift,k_shift) - (nu_kphalf(j_shift,k_shift)+nu_kmhalf(j_shift,k_shift))*u_n(j_shift,k_shift);
            u_np1(j_shift,k_shift) = u_np1(j_shift,k_shift) + nu_kmhalf(j_shift,k_shift)*u_n(j_shift,k_shift-1);
        end
    end
    
    %Condiciones de frontera
    % Frontera sur Neumann:
    u_np1(:,1) = 1-x_axis(:,1)/Lx; % k = 0
    % Frontera norte Neumann:
    u_np1(:,Ny+2) = 1-x_axis(:,1)/Lx; % k = Ny+1
    % Frontera oeste Dirichlet:
    u_np1(1,:) = 1; %sin(2*pi*3*y_axis(:,1)/Ly); % j = 0
    % Frontera este Neumann:
    u_np1(Nx+2,:) = u_np1(Nx+1,:); %j = Nx+1
    
    u_n = u_np1;
    
    t=t+dt
    
    if Movie == 1
        clf
        figure(1)
        hold on
        contour(x_axis(2:Nx+1),y_axis(2:Ny+1),transpose(u_np1(2:Nx+1,2:Ny+1)),50); colorbar;
        %plot3(x(2:Nx+1,2:Ny+1),y(2:Nx+1,2:Ny+1),u_np1(2:Nx+1,2:Ny+1)); 
        %caxis([0,0.5])
        %view(14,26)
        %plot(y_axis(2:Nx+1),u_np1(10,2:Nx+1))
        print('test','-depsc',figure(1))
        %stop
    end
end


% clf
% figure(1)
% subplot(2,1,1)
% hold on
% contour(x_axis,y_axis,transpose(kappa_kphalf),50); colorbar;
% print('Kappa','-depsc',figure(1))

clf
figure(1)
subplot(2,1,1)
hold on
contour(x_axis,y_axis,transpose(u_np1),100); colorbar;
%caxis([0,0.5])
print('SolT10','-depsc',figure(1))
