function [x, y, z] = cylinder( p, len, r)
% constuct the mesh for a cylinder
% p is the position of the center of the face
% l the length
% r the radius
axial = 20;
rotational = 50;
radial = 20;
ll = linspace(0, len, axial);
th1 = linspace(0, 2*pi, rotational);
[l, th] = meshgrid(ll, th1);
x1 = p(1) + l;
y1 = p(2) + r*cos(th);
z1 = p(3) + r*sin(th);
r1 = linspace(0, r, radial);
[r2, th] = meshgrid( r1, th1 );
x2 = p(1) * ones(size(r2));
y2 = p(2) + r2.*cos(th);
z2 = p(3) + r2.*sin(th);
x3 = (p(1) + len)*ones(size(r2)); 
[r3, th] = meshgrid(r1(end:-1:1), th1);
y3 = p(2) + r3.*cos(th);
z3 = p(3) + r3.*sin(th);
x = [x2 x1 x3];
y = [y2 y1 y3];
z = [z2 z1 z3];