function [] = icecream(N)

   h = 2/N;
   [x,y,z] = meshgrid(-1:h:1,-1:h:1,-1:h:1);
   phi1 = sqrt(x.^2+y.^2+z.^2)-0.5;
   phi2 = z;
   phi = min(phi1,phi2);
   phi3 = -z+sqrt(x.^2+y.^2)-0.5;
   phi = max(phi,phi3);
   clf
   isosurface(phi,0)

end
