% Plot points sampled from a torus

% Radii parameters
r = 1;
R = 3

% Number of points
L = 100;

% Input angles
psi = linspace(0,2*pi,L);
theta = linspace(0,2*pi,L);
%psi = 2*pi*rand(1,L);

% (x,y,z) coords defined in parametric coords
%h =@(x,y,z,r,k) (x.*(k*x-1).^2 .*(x-2) + y.^2).^2 + z.^2 - r.^2;
h =@(x,y,z,r,k) (x.*(x-1).^2 .*(k*x-2).^2 .*(x-3).^2 .*(x-4) + y.^2).^2 + z.^2 - r.^2;

% Define empty (X,Y,Z)
X = [];
Y = [];
Z = [];

% Twisted Band / Pinched Annulus when set p in Z to psi
% Make a pinched annulus % Y = ones(size(X));
for k = .1 %:.1:2
figure; hold on;

S = [];
M = 5000;
E = [];
for i = 1:4000
D = -5+10*rand(3,M);

D(4,:) = h(D(1,:),D(2,:),D(3,:),.01*ones(size(D(3,:))),k);

K = find(abs(D(4,:)) < .01);

D(1,K) = D(1,K) - ones(size(D(1,K)));

E = [E D(:,K)];

end
view(47,82);
plot3(E(1,:),E(2,:),E(3,:),'o')
pause(.1)
end


% Uncomment below to see 'gluing the patches together
U = figure;
%hold on;
S = [];
E(5,:) = 100;

for j = 1:1
	view(47-10+20*rand(1),20);
	% Random centers of basis balls
	P = -1+3*rand(3,1);
	r = .5*2.5*rand(1);
	k = 1;


	E(5,:) = sqrt( 	 (E(1,:) - P(1)*ones(1,length(M))).^2 + ...
			  	     (E(2,:) - P(2)*ones(1,length(M))).^2 + ...
   					 (E(3,:) - P(3)*ones(1,length(M))).^2 ); 

	N = find(E(5,:) < r); 

	
	
		figure(U)
		plot3(E(1,N),E(2,N),E(3,N),'or');
		%plot3(0,0,0)
		pause(1)
end

