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

fid = fopen('ex_cone.dat', 'wt');
fprintf(fid, '%12.3f', data);
fclose(fid);


