function fractal_landscape(k)
% k should be a power of 2, the landscape generated will be 2^k+1 by 2^k+1

% INITIALIZATION %
n = 2^k+1;
landscape = zeros(n);
boolmap = zeros(n);

n = n-1;
sizeLand = length(landscape);
randRange = [-1,1];
f = 0;

%Plot initial flat landscape and save to file
[X,Y] = meshgrid(1:length(landscape(1,:)), 1:length(landscape(:,1)));
%mesh(X,Y,landscape);
surf(X,Y,landscape, 'EdgeColor', 'none');
title('Fractal Landscape');
f = f+1;
%saveas(gcf, ['movie',sprintf('%03d',f),'.png'])
pause(.1)

while n > 1
    step = n/2;
    % DIAMOND STEP %   
    for i=step+1:n:sizeLand-1
        for j=step+1:n:sizeLand-1
            val = rand() * (randRange(2)-randRange(1)) + randRange(1);
            if boolmap(i,j) == 0
                landscape(i,j) = val + (landscape(i-step,j-step) + landscape(i-step,j+step) + landscape(i+step,j-step) + landscape(i+step,j+step))/4;
                boolmap(i,j) = 1;
            end
        end
    end
    
    %Plot intermediate landscape between diamond and square steps
    [X,Y] = meshgrid(1:length(landscape(1,:)), 1:length(landscape(:,1)));
    %mesh(X,Y,landscape);
    surf(X,Y,landscape, 'EdgeColor', 'none');
    title('Fractal Landscape');
    f = f+1;
    %saveas(gcf, ['movie',sprintf('%03d',f),'.png'])
    pause(.1)

    % SQUARE STEP %
    jinit = 1+step;
    j = jinit;
    for i=1:step:sizeLand
        while j <= sizeLand
            val = 0;
            ctr = 4;
            
            if i > 1
                val = val + landscape(i-step,j);
            else
                ctr = ctr-1;
            end             
            
            if i < sizeLand
                val = val + landscape(i+step,j);
            else
                ctr = ctr-1;
            end

            if j > 1
                val = val + landscape(i,j-step);
            else
                ctr = ctr-1;
            end

            if j < sizeLand
                val = val + landscape(i,j+step);
            else
                ctr = ctr-1;
            end
            
            val = val/ctr + rand() * (randRange(2)-randRange(1)) + randRange(1);
                
            if boolmap(i,j) == 0
                landscape(i,j) = val;
                boolmap(i,j) = 1;
            end
            
            j = j+step;
        end
        
        if jinit == 1+step
            jinit = 1;
        else
            jinit = 1+step;
        end
        j = jinit;
    end
    
    randRange = randRange/2;
    n = n / 2;
   
    %Plot result after square step
    [X,Y] = meshgrid(1:length(landscape(1,:)), 1:length(landscape(:,1)));
    %mesh(X,Y,landscape);
    surf(X,Y,landscape, 'EdgeColor', 'none');
    title('Fractal Landscape');
    f = f+1;
    %saveas(gcf, ['movie',sprintf('%03d',f),'.png'])
    pause(.1)
end

%Plot final result
%[X,Y] = meshgrid(1:length(landscape(1,:)), 1:length(landscape(:,1)));
%mesh(X,Y,landscape);
%surf(X,Y,landscape, 'EdgeColor', 'none');

