clear all; 
close all;
a = 2;
R = 1;
Nr = 25;
r = linspace(a-R,a+R,Nr);
Ntheta = 25;
theta = linspace(0,2*pi,Ntheta);
% outer products can be used to create matrices for x and y
% the product of N-by-1 (column) vector r' and 1-by-M (row) vector cos(theta) creates a matrix N-by-M
x = r'*cos(theta); 
y = r'*sin(theta);
[n,m] = size(x);

figure(1);
view([0,0,1]); hold on;
mesh(x,y,zeros(n,m)); hold off;

Urad = 0.1*sqrt(R^2-(r-a).^2);
Umax = max(Urad);
u = Urad'*ones(1,m);

figure(2); 
V = linspace(0,Umax,5);
h = contour(x,y,u,V); 
clabel(h); 

figure(3);
axis([-(R+a),(R+a),-(R+a),(R+a),0,4*Umax]); hold on;
xlabel('x'); ylabel('y'); zlabel('u'); title('The torus');
view(-20,20);
surf(x,y,u+2*Umax); 
surf(x,y,-u+2*Umax); 
mesh(x,y,zeros(n,m)); hold off;