

N=101;                              %%Size of grid (resolution)
L=10;                               %%length of axis
height=1;                           %%original height scale
p=.75;                               %%parameter: when p is high, new 
                                    %%pyramids have a bigger effect

q=5;                                %%positive integer, how many new 
                                    %%pyramids added before one is 
                                    %%subtracted
                                    
T=zeros(N,N);                       %%Mountain grid


for iter=1:(20*q +1)
     
H= (p*height)/(-2)^(mod(iter,q));   %%Height of new pyramid
A=floor(rand()*N)+1;                %%coordinates of
B=floor(rand()*N)+1;                %%pyramid peak (A,B)

      

%%%%%%%This section creates the pyramid%%%%%%%%%%%
T1=zeros(N,N);                 
T1(A,:) = H;
    
for i=1:A-1
    T1(i,:) = (i-1)*(H/(A-1));
end
for i=A+1:N
    T1(i,:) = H - (i-A)*(H/(N-A));
end
T2=zeros(N,N);    
T2(:,B) = H; 
for i=1:B-1
    T2(:,i) = (i-1)*(H/(B-1));
end
for i=B+1:N
    T2(:,i) = H - (i-B)*(H/(N-B));
end

%%%%%%This creates the new mountain%%%%%%%%%%%%%%%

Tnew=min(T1,T2);
T=Tnew+T;                       

%%%%%%%%%%Finds height of mountain%%%%%%%%%%%%%%%%


for i=1:N                       %%This part finds 
    for j=1:N                   %%the height of the
    if T(i,j)<0                 %%mountain in this
        T(i,j)=abs(T(i,j));     %%iteration, which
    end                         %%affects the height
    if T(i,j)>height            %%of the next
        height = T(i,j);        %%pyramid. This 
    end                         %%makes sure that
    end                         %%the mountain will
end                             %%have more dramatic
                              

%%%%%%%This section plots%%%%%%%%%%%%%%%%%%%%%%%%%%
xx = 0:L/(N-1):L;
yy = 0:L/(N-1):L;
[X,Y] = meshgrid(yy,xx);
surf(X,Y,T);
view([+34.5 14]);
pause(0.2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end
