clf
clear

xmin = 0;
xmax = 60;
N=100;
x = xmin:xmax/N:xmax;
y = x;
f=0;
T = zeros(N+1,N+1);

k = 0.0088;
dx = xmax/N;
dy = dx;
eta = 0.24; %Reasonable constant for 2D
dt = (eta/k)*(dx)^2;

[X,Y] = meshgrid(x,y);

%Set up dike and sill
for i=1:N+1
    for j=1:N+1
        if (x(i)<=15) && (x(i)>=10)
            if (y(j)>=0) && (y(j)<=60)
                T(j,i) = 1000;
            end
        end
        if (x(i)<=60) && (x(i)>=15)
            if (y(j)>=20) && (y(j)<=25)
                T(j,i) = 1000;
            end
        end
    end
end

for iter=1:400
    Tnew=T;
    for k=2:N
        for l=2:N
            Tnew(k,l) = eta*(T(k+1,l)+T(k-1,l)+T(k,l+1)+T(k,l-1)-4*T(k,l))+T(k,l); 
        end
    end
    T=Tnew;
    
    %Set up insulating boundary conditions
    for m=2:N
        for n=2:N
            T(1,n) = T(2,n);
            T(N+1,n) = T(N,n);
            T(m,1) = T(m,2);
            T(m,N+1) = T(m,N);
            T(1,1) = T(m,1);
            T(N+1,1) = T(N,1);
        end
    end
    
    %Plot distribution and save images to make into gif
    if mod(iter+39,10) == 0
        f=f+1;
        surf(X,Y,T)
        title('Temperature Distribution')
        colormap jet
        axis([0 xmax 0 xmax 0 1000])
        caxis([0 1000])
        colorbar;
        zlabel('Temperature')
        xlabel('Distance')
        ylabel('Distance')
        saveas(gcf, ['final_',sprintf('%03d',f),'.png'])
        pause(0.01)
        
        contourf(X,Y,T)
        colormap jet
        caxis([0 1000])
        colorbar
        saveas(gcf, ['topview_',sprintf('%03d',f),'.png'])
        pause(0.01)
    end
end
