

clear
close all
ii = 1;
for i = 1:1


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

% Create a random point at which to center the sphere
a  = .2;
x0 = -a+2*a*rand(1);
y0 = -a+2*a*rand(1);
z0 = -a+2*a*rand(1);

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

% 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;

% The sphere at (x0,y0,z0)
hs = sqrt((xx-x0).^2 + (yy-y0).^2 + (zz-z0).^2);

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

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

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

% Close previous figure
%close all
%set(0,'DefaultFigureWindowStyle','docked')

% Allowing a hold on, and iterating the procedure we see the surface
% 'glued back together again'
%hold on;

%# create the isosurface by thresholding at a iso-value of 10
disp('isosurface'); tic;


isosurface(xx,yy,zz,hs,1);
%isosurface(xx,yy,zz,hs,1);
axis off; hold on;
toc;
%axis equal 
% fancy effects

       %p = patch(isosurface(xx, yy, zz, h,0.25));
       %isonormals(xx,yy,zz,h, p)
       %set(p, 'FaceColor', 'y', 'EdgeColor', 'none');
       daspect([1 1 1]); 
       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]);
       %view([118,18])


for i=1:20
  camorbit(1,0,'data',[0 1 1])
  M(ii) = getframe(1);
  java.lang.Runtime.getRuntime.freeMemory
  ii = ii+1;
end

end
movie2avi(M,'surface_rotate.avi');
%movie(M)


