nx = 11;
ny = 11;
nz = 11;

data = zeros(nx, ny, nz);
fid = fopen('ex_cone.dat', 'r');
for iz=1:nz
    data(1:nx,1:ny,iz) = fscanf(fid, '%g', [nx, ny]);
end
fclose(fid);

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

thresh = 3.0;
fv = isosurface(data, thresh);
p = patch(fv);
n = isonormals(d, 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);
