clear
clf

N = 100;
initDistriFileName = 'face.obj';

% Discover what L, min x/y, max x/y are from the file
minX = Inf;
minY = Inf;
maxX = -Inf;
maxY = -Inf;
initDistriFileID = fopen(initDistriFileName);
nextLine = fgetl(initDistriFileID);
while (isa(nextLine, 'char') || (strcmp(nextLine, '') == 1))
    lineTokens = strsplit(nextLine);
    if (lineTokens{1} == 'v')
        x = str2double(lineTokens{2});
        y = str2double(lineTokens{3});
        if (x > maxX)
            maxX = x;
        end
        if (x < minX)
            minX = x;
        end
        if (y > maxY)
            maxY = y;
        end
        if (y < minY)
            minY = y;
        end
    end
    nextLine = fgetl(initDistriFileID);
end
fclose(initDistriFileID);
L = max(maxX-minX, maxY-minY);
dx = L/N;
dy = L/N;
xx = minX:dx:minX+L;
yy = minY:dy:minY+L;
[X,Y] = meshgrid(yy,xx);

% Actually construct initial distribution from file info
T = zeros(N+1);
T(:,:) = -Inf;
initDistriFileID = fopen(initDistriFileName);
nextLine = fgetl(initDistriFileID);
while (isa(nextLine, 'char') || (strcmp(nextLine, '') == 1))
    lineTokens = strsplit(nextLine);
    if (lineTokens{1} == 'v')
        x = str2double(lineTokens{2});
        y = str2double(lineTokens{3});
        z = str2double(lineTokens{4});
        i = round((x-minX)/dx)+1;
        j = round((y-minY)/dy)+1;
        if (z > T(i,j))  % only use highest z-value for initial distribution, to ensure continuity of model if it's a closed surface
            T(i,j) = z;
        end
    end
    nextLine = fgetl(initDistriFileID);
end
fclose(initDistriFileID);

% Find initial z bounds to assist in graphing
minZ = Inf;
maxZ = -Inf;
for i = 1:size(T,1)
    for j = 1:size(T,2)
        if (T(i,j) ~= -Inf)  % This was filled
            if (T(i,j) > maxZ)
                maxZ = T(i,j);
            end
            if (T(i,j) < minZ)
                minZ = T(i,j);
            end
        end
    end
end

% Set all unfilled parts to have minimum z value
for i = 1:size(T,1)
    for j = 1:size(T,2)
        if (T(i,j) == -Inf)  % This was not filled
            %numNeighbors = 0;
            %totalSumZ = 0;
            %if (i > 1)
            %    totalSumZ = totalSumZ + max(T(i-1,j),minZ);
            %    numNeighbors = numNeighbors + 1;
            %end
            %if (i < size(T,1))
            %    totalSumZ = totalSumZ + max(T(i+1,j),minZ);
            %    numNeighbors = numNeighbors + 1;
            %end
            %if (j > 1)
            %    totalSumZ = totalSumZ + max(T(i,j-1),minZ);
            %    numNeighbors = numNeighbors + 1;
            %end
            %if (j < size(T,2))
            %    totalSumZ = totalSumZ + max(T(i,j+1),minZ);
            %    numNeighbors = numNeighbors + 1;
            %end
            %T(i,j) = totalSumZ / numNeighbors;
            T(i,j) = minZ;
        end
    end
end

k = 1;
dt = 0.02;
numIterations = 2000;
eta = (k*dt)/(dx^2);
f = 0;
for iter=1:numIterations
    % Dirichlet boundary conditions
    T(:,1) = minZ;
    T(:,N+1) = minZ;
    T(1,:) = minZ;
    T(N+1,:) = minZ;
    % insulating boundary conditions
    %T(:,1) = T(:,2);
    %T(:,N+1) = T(:,N);
    %T(1,:) = T(2,:);
    %T(N+1,:) = T(N,:);
    % +periodic boundary conditions
    %T(:,N+1) = T(:,1);
    %T(N+1,:) = T(1,:);
    
    Tnew=T; % just to copy value at the boundaries at 1 and N+1
    for i=2:N
        for j=2:N
            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;
    if mod(iter+9,10) == 0
        f=f+1;
        surf(X,Y,T)
        %shading interp
        zlim([minZ, maxZ])
        scrsz = get(0,'ScreenSize');
        set(gcf,'Position',[1 1 scrsz(3) scrsz(4)])
        saveas(gcf, ['movie1_',sprintf('%03d',f),'.png'])
        pause(0.01)
    end
end