N = 8;
sideLength = 2^N + 1;
terrain = zeros(sideLength, sideLength);
rng(6,'twister');

for i = 1:N
    terrain = doDiamond(terrain, N, i);
    terrain = doSquare(terrain, N, i);
end

%surfc(terrain);
set(gca,'color',[.5 .5 .5])
%colormap(blue)
h=surfl(terrain)
set(h,'edgecolor','none')
%contourf(peaks)


function terrain = doDiamond(terrain, N, iterNum)
    delta = 2^(N-iterNum);
    numBox = 2^(iterNum-1);
    boxSize = 2^(N)/(2^(iterNum-1));
    
    %{
    1 x x x 2
    x x x x x
    x x c x x
    x x x x x
    3 x x x 4
    %}
    
    for box_row = 1:numBox
        for box_column = 1:numBox
    
            center = [(box_row - 1)*(boxSize)+delta+1, ((box_column - 1)*boxSize)+delta+1];
            corner1 = [center(1) - delta, center(2) - delta];
            corner2 = [center(1) + delta, center(2) - delta];
            corner3 = [center(1) - delta, center(2) + delta];
            corner4 = [center(1) + delta, center(2) + delta];
            
            terrain(center(1),center(2)) = (terrain(corner1(1),corner1(2)) + ...
                terrain(corner2(1),corner2(2)) + ...
                terrain(corner3(1),corner3(2)) + ...
                terrain(corner4(1),corner4(2)))/4; 
            terrain(center(1),center(2)) = terrain(center(1),center(2)) + ...
                ((rand() - 0.5)*2)*2^(-iterNum);
        end
    end
end

function terrain = doSquare(terrain, N, iterNum)
    delta = 2^(N-iterNum);
    numBox = 2^(iterNum-1);
    boxSize = 2^(N)/(2^(iterNum-1));
    
    %{
    x x 1 x x
    x x x x x
    2 x c x 3
    x x x x x
    x x 4 x x
    %}
    
    for box_row = 1:numBox
        for box_column = 1:numBox
            center = [(box_row - 1)*(boxSize)+delta+1, (box_column - 1)*(boxSize)+delta+1];
            sq = zeros(4, 2);
            sq(1,1) = center(1);
            sq(1,2) = center(2) - delta;
            sq(2,1) = center(1) - delta;
            sq(2,2) = center(2);
            sq(3,1) = center(1) + delta;
            sq(3,2) = center(2);
            sq(4,1) = center(1);
            sq(4,2) = center(2) + delta;
           
            for sq_num = 1:4
                if(box_row ~= 1 && sq_num == 1)
                    continue;
                end
                if(box_column ~= 1 && sq_num == 2)
                    continue;
                end

                corner1 = [sq(sq_num, 1), sq(sq_num, 2) - delta];
                corner2 = [sq(sq_num, 1) - delta, sq(sq_num, 2)];
                corner3 = [sq(sq_num, 1) + delta, sq(sq_num, 2)];
                corner4 = [sq(sq_num, 1), sq(sq_num, 2) + delta];
               
                if(sq(sq_num, 1) - delta < 1)
                    %corner2(1) = length(terrain) - delta;  
                    terrain(sq(sq_num, 1),sq(sq_num,2)) = (terrain(corner1(1),corner1(2)) + ...
                    terrain(corner3(1),corner3(2)) + ...
                    terrain(corner4(1),corner4(2)))/3; 
                %end
                elseif(sq(sq_num, 2) - delta < 1)
                    %corner1(2) = length(terrain) - delta;
                    terrain(sq(sq_num, 1),sq(sq_num,2)) = (terrain(corner4(1),corner4(2)) + ...
                    terrain(corner2(1),corner2(2)) + ...
                    terrain(corner3(1),corner3(2)))/3; 
                %end
                elseif(sq(sq_num, 1) + delta > length(terrain))
                    %corner3(1) = 1 + delta;
                    terrain(sq(sq_num, 1),sq(sq_num,2)) = (terrain(corner1(1),corner1(2)) + ...
                    terrain(corner2(1),corner2(2)) + ...
                    terrain(corner4(1),corner4(2)))/3; 
                %end
                elseif(sq(sq_num, 2) + delta > length(terrain))
                   % corner4(2) = 1 + delta;  
                   terrain(sq(sq_num, 1),sq(sq_num,2)) = (terrain(corner1(1),corner1(2)) + ...
                    terrain(corner2(1),corner2(2)) + ...
                    terrain(corner3(1),corner3(2)))/3; 
                else 
                    terrain(sq(sq_num, 1),sq(sq_num,2)) = (terrain(corner1(1),corner1(2)) + ...
                    terrain(corner2(1),corner2(2)) + ...
                    terrain(corner3(1),corner3(2)) + ...
                    terrain(corner4(1),corner4(2)))/4; 
                end
                
                
                terrain(sq(sq_num, 1),sq(sq_num,2)) = terrain(sq(sq_num, 1),sq(sq_num,2)) + ...
                ((rand() - 0.5)*2)*2^(-iterNum);
                
            end
            
        end
    end
end