function [x, y, z] = RandomCylinder(ravg, zmin, zmax, ni, nj)
  % -----------------------------------------------------------------------
  % RandomCylinder - define a cylinder of randomly varying radius, with its
  %                  axis coinciding with the z axis.
  %
  % Args: ravg - average radius
  %       zmin - first point on axis of cylinder
  %       zmax - last point on axis of cylinder
  %       ni - number of gridlines around the cylinder
  %       nj - number of gridlines along the axis
  %
  % Rtns: x - ni x nj matrix of x coordinates of grid points
  %       y - ni x nj matrix of y coordinates of grid points
  %       z - ni x nj matrix of z coordinates of grid points
  % -----------------------------------------------------------------------

  % let r be some random function...
  r = 2 * ravg * rand(nj, 1);
  % smooth it a bit...
  for k = 1:20; r = ([0;0;r] + [0;r;0] + [r;0;0])/3; r = r(2:nj+1); end;

  % r is a column vector, theta is a row vector, outer product is a matrix
  theta = 0:2*pi/(ni-1):2*pi;
  x = r * cos(theta);
  y = r * sin(theta);

  % z is a column vector, so we repeat it ni times to make the matrix
  z = repmat((zmin + (zmax-zmin)*(0:nj-1)/(nj-1))', 1, ni);
end
