n1 = 11;
n2 = 11;
n3 = 11;
data = zeros(n1, n2, n3);
s = size(data);
c = (s+1)/2;

for i1=1:n1,
    for i2=1:n2,
	for i3=1:n3,
	    vd = [i1, i2, i3] - c;
	    d(i1,i2,i3) = sqrt(dot(vd,vd));
	    if (i3 >= n3/2),
	        vd(3) = 0;
		d(i1,i2,i3) = min(d(i1,i2,i3), 2.0*sqrt(dot(vd,vd)));
	    end
	end
    end
end

thresh = 3.0;
fv = isosurface(d, thresh);
p = patch(fv);
isonormals(d, p);
daspect([1 1 1]);
view(3);
camlight;
lighting phong;
axis([0 10 0 10 0 10 0 255]);
set(p,'FaceColor','red','EdgeColor','none');
set(p,'AmbientStrength',0.2);
set(p,'SpecularColorReflectance',0.5,'SpecularExponent',25);

