%% 253 demo: isosurfaces
% 2015 W1: Colin Macdonald, public domain

x1 = linspace(-2, 2, 64);
y1 = x1;
z1 = x1;
[x,y,z] = meshgrid(x1,y1,z1);

figure(1); clf;
title('isosurface ("3D contour") example')
xlabel('x'); ylabel('y'); zlabel('z')

%% you can try different functions here
f = x.^2 + y.^2 + z.^2;
%f = x.^2 + y.^3 + sin(z);

% crude way to cut a slice out
f(x < 0 & x > y) = nan;

% first, plot the k = 1 isosurface
k = 1

%isosurface(x,y,z,f, k);  % simple version

				% but to control color, we use:
tic
H = patch(isosurface(x, y, z, f, k));
toc
tic
isonormals(x, y, z, f, H)
toc
tic
set(H, 'FaceColor', 'red', 'EdgeColor', 'none');
toc

view(3)
grid on
% comment this out if using octave
camlight left
lighting phong
axis equal
xlim([-1.2, 1.2])
ylim([-1.2, 1.2])
zlim([-1.2, 1.2])


disp('press enter to add another isosurface')
pause

% now plot the k=2 surface
k = 2
H = patch(isosurface(x, y, z, f, k));
isonormals(x,y,z,f, H)
set(H, 'FaceColor', 'blue', 'EdgeColor', 'none');
set(H, 'FaceAlpha', 0.6);


disp('press enter to add another isosurface')
pause

% now plot the k=2 surface
k = 1/3
H = patch(isosurface(x, y, z, f, k));
isonormals(x,y,z,f, H)
set(H, 'FaceColor', 'green', 'EdgeColor', 'none');
set(H, 'FaceAlpha', 0.6);
%axis normal
%axis equal
lighting phong
