function B=land(A)  % program for drawing random landscapes
%   Adapted from "Matlab Guide" by Desmond and Nicholas Higham (p.95-97)
%   published by SIAM (www.siam.org)

%  A is square of dimension N=2^k+1
%  Four corners of A are used as input parameters.
%
% To call program, type:
%
%     >> rand('state',10);           %  this will make the outputs repeatable, otherwise a
%                                               %  different landscape will be produced at each run (if thats what
%                                               %   youwant).
%     >> k=2^5+1;                        %  or 2^n+1 for any positive integer n
%     >> A=zeros(k);
%     >> A([1 k],[1 k])=[1   1.25
%                         1.1  2.0];            %       this is the initial
%                                                    %       displacements of the four courners
%
%     >>  B=land(A);                         %    the fractal landscape is the
%                                                       %     matrix B
  


N=size(A,1);
f=(N+1)/2;
level=log2(N-1);
scalef=0.05*(2^(0.9*level));   % can adjust the 0.05 and the 0.9 for different effects 

                            
 B=A;

B(f,f)=mean([A(1,1),A(1,N),A(N,1),A(N,N)])+scalef*randn;  
                           %  randn is a normal random number;
                           %  you may want to use rand for a uniform (in particular, positive) 
                            %    random number, or use abs(randn) for positive displacements.
B(1,f)=mean([A(1,1),A(1,N)])+scalef*randn;
B(f,1)=mean([A(1,1),A(N,1)])+scalef*randn;
B(f,N)=mean([A(1,N),A(N,N)])+scalef*randn;
B(N,f)=mean([A(N,1),A(N,N)])+scalef*randn;

if N>3
    B(1:f,1:f)=land(B(1:f,1:f));
    B(1:f,f:N)=land(B(1:f,f:N));
    B(f:N,1:f)=land(B(f:N,1:f));
    B(f:N,f:N)=land(B(f:N,f:N));
end