
 % Colleen Murphy
 % EPS 109 Independent Project
 % Fall 2015
 % 2D heat diffusion in real time with Dirichlet Boundary Conditions. The
 % initial conditions are set by loading an image, binarizing it, and
 % using the 1 values as the regions elevated in temperature.
 
clf ; clear ; close all

% Loading and binarizing the image
H = imread('cal.jpg');
I = imrotate(H,180);
G = rgb2gray(I);
level = graythresh(I);
bw = im2bw(I,level);
figure, imshow(bw)
imwrite(bw,'my_cal.jpg');

% Setting up the grid and the constants
[x,y,z] = size(bw);
N = x;
M = y;
dx = x/N;
dy = y/M;
L = 0:dx:x;
W = 0:dy:y;
eta = 0.24;
k = 1;
dt = eta/k;

% Setting the initial height of the object for heat diffusion
T = double(bw);

% Diffuse!
f = 0;
for iter = 1:2000
    iter
    Tnew = T;
    if mod(iter+9,10) == 0
        f= f+1;
        subplot(1,2,1)
        xx = 0:x/N:x-1;
        yy = 0:y/M:y-1;
        [X,Y] = meshgrid(yy,xx);
        surf(X,Y,T)
        set(gca,'xdir','reverse')
        shading flat
        colormap('default')
        zlim([0 1])
        pause(0.001)

        subplot(1,2,2)
        shading interp
        contour(X,Y,T)
        set(gca,'xdir','reverse')
        
        saveas(gcf, ['cal_',sprintf('%03d',f),'.tif'])

    end
    for i = 2:N-1
        for j = 2:M-1
            Tnew(i,j) = eta*(T(i+1,j)+ T(i-1,j) + T(i, j+1) + T(i, j-1)-4*T(i,j))+T(i,j);
        end
    end
    T = Tnew;
    
    subplot(1,2,1)
    xx = 0:x/N:x-1;
    yy = 0:y/M:y-1;
    [X,Y] = meshgrid(yy,xx);
    surf(X,Y,T)
    set(gca,'xdir','reverse')
    shading flat
    colormap('default')
    zlim([0 1])
    pause(0.001)
    
    subplot(1,2,2)
    shading interp
    contour(X,Y,T)
    set(gca,'xdir','reverse')
end
