clear all 
close all 
clf 

L=20;
N=100;

Z=zeros(N+1,N+1);

for i=1:N+1
    for j=1:N+1
    %courtyard 
    Z(5:10,20:60)=20;
    Z(10:45,20:25)=20;
    Z(40:45,20:60)=20;
    Z(10:45,55:60)=20;
    %other part 
    Z(45:50,20:60)=40;
    Z(45:85,20:25)=40;
    Z(80:85,20:60)=40;
    Z(45:85,55:60)=40;
    dr=sqrt((i-65)^2+(j-40)^2);
    h=40;
    r=8;
    if dr<r
        Z(i,j)=40+h*(cos(dr/r))^2;
    end
    %inner square
    Z(73:75,32:48)=60;
    Z(57:73,30:32)=60;
    Z(55:57,32:48)=60;
    Z(57:73,48:50)=60;
    dr=sqrt((i-75)^2+(j-48)^2);
    h=45;
    r=2;
    if dr<r
        Z(i,j)=35+h*(cos(dr/r))^2;
    end 
    dr=sqrt((i-55)^2+(j-48)^2);
    h=45;
    r=2;
    if dr<r
        Z(i,j)=35+h*(cos(dr/r))^2;
    end 
    dr=sqrt((i-75)^2+(j-30)^2);
    h=40;
    r=2;
    if dr<r
        Z(i,j)=35+h*(cos(dr/r))^2;
    end 
    dr=sqrt((i-55)^2+(j-32)^2);
    h=45;
    r=2;
    if dr<r
        Z(i,j)=40+h*(cos(dr/r))^2;
    end
    %domes 
    dr=sqrt((i-53)^2+(j-28)^2);
    h=35;
    r=3;
    if dr<r
        Z(i,j)=25+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-53)^2+(j-40)^2);
    h=35;
    r=3;
    if dr<r
        Z(i,j)=25+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-53)^2+(j-52)^2);
    h=35;
    r=3;
    if dr<r
        Z(i,j)=25+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-65)^2+(j-28)^2);
    h=35;
    r=3;
    if dr<r
        Z(i,j)=25+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-77)^2+(j-28)^2);
    h=35;
    r=3;
    if dr<r
        Z(i,j)=25+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-77)^2+(j-40)^2);
    h=35;
    r=3;
    if dr<r
        Z(i,j)=25+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-77)^2+(j-53)^2);
    h=35;
    r=3;
    if dr<r
        Z(i,j)=25+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-65)^2+(j-53)^2);
    h=35;
    r=3;
    if dr<r
        Z(i,j)=25+h*(cos(dr/r))^2;
    end
    %minarets 
    dr=sqrt((i-45)^2+(j-18)^2);
    h=55;
    r=3;
    if dr<r
        Z(i,j)=45+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-45)^2+(j-62)^2);
    h=55;
    r=3;
    if dr<r
        Z(i,j)=45+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-85)^2+(j-20)^2);
    h=55;
    r=3;
    if dr<r
        Z(i,j)=45+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-85)^2+(j-62)^2);
    h=55;
    r=3;
    if dr<r
        Z(i,j)=45+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-10)^2+(j-20)^2);
    h=55;
    r=3;
    if dr<r
        Z(i,j)=45+h*(cos(dr/r))^2;
    end
    dr=sqrt((i-10)^2+(j-60)^2);
    h=55;
    r=3;
    if dr<r
        Z(i,j)=45+h*(cos(dr/r))^2;
    end
    end 
end 

xx=0:L/N:L;
yy=0:L/N:L;
[X,Y]=meshgrid(yy,xx);
surf(X,Y,Z);
axis([0 L 0 L 0 100]);