
% Close previous figure
close

% Define important granularity of mesh.  Too small gives bad cells.
% Too large freezes the computer.
k = .01;

% Create a mesh/grid of points on which to calculate level function
[xx,yy,zz] = meshgrid(-1:k:1,-1:k:1,-1:k:1);

% Define radii parameters for torus
r = .15*ones(size(xx));
R = 3*ones(size(xx));

% Genus 4 surface
%h = (x.*(x-1).^2 .*(k*x-2).^2 .*(x-3).^2 .*(x-4) + y.^2).^2 + z.^2 - r.^2;

% Setting h = 0 gives genus 2 surface
h = ((xx+1).*(xx).^2 .*(xx-1) + yy.^2).^2 + zz.^2 - r.^2;

% Sphere
h = sqrt(xx.^2 + yy.^2 + zz.^2);

% The norm of each mesh gridpoint
%hs = sqrt(xx.^2 + yy.^2 + zz.^2);

% First set all function values not close to surface far away from 0
%h(find(abs(h)>.01))  = 2;

% Find all points on the surface lying outside a ball of radius r
%h(intersect(find(h < 2),find(hs >= .51))) = 3;

% The torus function
%h = (xx.^2 + yy.^2 + zz.^2 + R.^2 - r.^2).^2 - 4*R.^2 .* (xx.^2 + yy.^2);

%# create the isosurface by thresholding at a iso-value of 10
isosurface(xx,yy,zz,h,1);

%axis equal 
% fancy effects

       %p = patch(isosurface(xx, yy, zz, h,0));
       %isonormals(xx,yy,zz,h, p)
       %set(p, 'FaceColor', 'interp', 'EdgeColor', 'none');
       daspect([1 1 1]); 
	   axis tight; 
       colormap(prism(28))
       camup([1 0 0 ]); campos([25 -55 5]) 
       camlight; lighting phong
	   %xlim([min(min(min(xx)))-.5, max(max(max(xx)))+.5]);
        

