nx = 11;
ny = 11;
nz = 11;
xmin =  0.0;
xmax = 10.0;
ymin =  0.0;
ymax = 10.0;
zmin =  0.0;
zmax = 10.0;
pmin = [xmin ymin zmin];
pmax = [xmax ymax zmax];
prange = pmax - pmin;
pcenter = (pmin + pmax)/2.0;

x = zeros(nx, ny, nz);
y = zeros(nx, ny, nz);
z = zeros(nx, ny, nz);
for ix=1:nx,
    for iy=1:ny,
	for iz=1:nz,
	    x(ix, iy, iz) = ix;
	    y(ix, iy, iz) = iy;
	    z(ix, iy, iz) = iz;
	end
    end
end

data = zeros(nx, ny, nz);
for ix=1:nx,
    for iy=1:ny,
	for iz=1:nz,
	    p = [x(ix, iy, iz) y(ix, iy, iz) z(ix, iy, iz)];
	    d2 = p - pcenter;
	    data(ix,iy,iz) = sqrt(dot(d2,d2));
            if (iz >= nz/2),
 	        d2(3) = 0;
 		data(ix,iy,iz) = min(data(ix,iy,iz), 2.0*sqrt(dot(d2,d2)));
 	    end
	end
    end
end

[xgrad, ygrad, zgrad] = gradient(data);

thresh = 3.0;
fv = isosurface(data, thresh);
p = patch(fv);
n = isonormals(data, p);
daspect([1 1 1]);
view(3);
camlight;
lighting gouraud;
axis([0 10 0 10 0 10 0 255]);
set(p,'FaceColor','red','EdgeColor','none');

ncone = 100;
cx = zeros(ncone, 1);
cy = zeros(ncone, 1);
cz = zeros(ncone, 1);

for i=1:100
    cx(i) = prange(1)*rand() + pmin(1);
    cy(i) = prange(2)*rand() + pmin(2);
    cz(i) = prange(3)*rand() + pmin(3);
end

cones = coneplot(x, y, z, xgrad, ygrad, zgrad, cx, cy, cz);
