
%set up size of plot
psize = 50; %x and y dimension
ax = [0 psize 0 psize -10 10]; %axis for later

%Initial conditions for Dam_1
Land = zeros(50,50);
Land(1:25,1:25)= peaks(25);
Land(26:50,1:25) = peaks(25);
Land(:,26:50) = 10;
Water = Land-0.3;
Land(3:47,29:47) = 5;
Both = max(Land,Water); %This is the higher of Land and Water 

%make a colormap for land and water. brown and blue basically
colormap([205/255,133/255,63/255;0/255,191/255,255/255])

%plot Land and Water
L = surf(Land); %L for Land
axis(ax)        %use ax from earlier for axis, so axis doesnt auto change
axis off        %hide axis
hold on
W = surf(Water);%W for Water
axis(ax)
axis off
hold off
alpha(W,0.5)    %This makes the water transparent!

set(W,'FaceColor','interp','EdgeColor','interp') %Changes grid color

%This allows the water and land to be different colors
C1 = ones(psize,psize); %for 1st row of colormap
C2 = C1 + 1;            %second row
set(L,'CData',C1); %as far as I understand, we're making the Land use the
set(W,'CData',C2); %first row of the colormap, and the water the second row

drawnow
pause 

%Now that the initial conditions are made, we start the main loop
for loop = 1:1000
    
    %for the gif:
%     loop
%     if loop <= 550        
%         saveas(gcf, ['Movie2_',num2str(loop, '%03d'),'.png'])       
%     end
% 
%     if loop > 550 && mod(loop, 50) == 0
%         saveas(gcf, ['Movie2_',num2str((550 + (loop - 550)/50), '%03d'),'.png'])
%     end
    
    if loop == 20 %after 20 loops, blow a hole in the dam!
        Land(24:27,26:29) = 7;
    end
     
    %We want the boundaries to be like walls now 
    %Water(1,:) = Land(1,:);
    %Water(psize,:) = Land(psize,:);
    %Water(:,1) = Land(:,1);
    %Water(:,psize) = Land(:,psize);
    
    Both = max(Land,Water);
    
    %move the water! Basically check 4 adjacent spots and find which is
    %lowest. Move some water in that direction. I removed half of
    %difference from current spot and added it to the spot where the water
    %is moving to, so the elevations end up equaling
    for i = 2:psize-1
        iL = i-1; %for use below, basically left of current spot
        iR = i+1; %right of current spot
        
        for j = 2:psize-1
            jD = j-1; %down from current spot
            jU = j+1; %up from current spot
            
            %if the water at this spot is above land and theres an adjacent
            %spot that has lower elivation than this one. Otherwise, this
            %water is in a pool of water and doesn't need to move.
            if Water(i,j) > Land(i,j) && (Both(iL,j) < Both(i,j) || Both(iR,j) < Both(i,j) || Both(i,jU) < Both(i,j) || Both(i,jD) < Both(i,j))
                
                adjacent = [Both(iL,j),Both(iR,j),Both(i,jU),Both(i,jD)];
                if Both(iL,j) == min(adjacent) %if left is lowest spot
                    slope =  Both(i,j) - Both(iL,j);
                    Water(i,j) = Water(i,j) - slope/2;
                    Water(iL,j) = Water(iL,j) + slope/2;
                end
                if Both(iR,j) == min(adjacent) %if right is lowest spot
                    slope =  Both(i,j) - Both(iR,j);
                    Water(i,j) = Water(i,j) - slope/2;
                    Water(iR,j) = Water(iR,j) + slope/2;
                end
                if Both(i,jU) == min(adjacent) %if up is lowest spot
                    slope =  Both(i,j) - Both(i,jU);
                    Water(i,j) = Water(i,j) - slope/2;
                    Water(i,jU) = Water(i,jU) + slope/2;
                end
                if Both(i,jD) == min(adjacent) %if down is lowest spot
                    slope =  Both(i,j) - Both(i,jD);
                    Water(i,j) = Water(i,j) - slope/2;
                    Water(i,jD) = Water(i,jD) + slope/2;
                end
                %notice that if two directions are equally the minimum, the
                %lack of 'else' here means it will go both directions
                
                %if the water goes below Land level, bring it back up
                if Water(i,j) < Land(i,j) 
                    Water(i,j) = Land(i,j);
                end
                if Water(i,j) == Land(i,j) %put it back down a little
                    Water(i,j) = Land(i,j)-0.3;%so you cant see it
                end
                Both = max(Land,Water);   
            end  
        end
    end
    
    %plot Land and Water
    L = surf(Land); %L for Land
    axis(ax)        %use ax from earlier for axis, so axis doesnt change
    axis off        %hide axis
    hold on
    W = surf(Water);%W for Water
    axis(ax)
    axis off
    hold off
    alpha(W,0.5)    %This makes the water transparent!

    set(W,'FaceColor','interp','EdgeColor','interp') %Changes grid color

    % This part allows the Water and Land to be different colors
    set(L,'CData',C1); %without this the colors sort of mix
    set(W,'CData',C2);

    drawnow
    
    pause (0.1)

end