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