% Plot points sampled from a torus

% Radii parameters
r = 1;
R = 3

% Number of points
L = 50;

% Input angles
psi = linspace(0,2*pi,L);
theta = linspace(0,2*pi,L);
[psi,theta]=meshgrid(0:.5:2.1*pi);

%psi = 2*pi*rand(1,L);

% (x,y,z) coords defined in parametric coords
x = @(R,theta,psi) (R+r*cos(psi)).*cos(theta);
y = @(R,theta,psi) (R+r*cos(psi)).*sin(theta);
z = @(r,psi) r*sin(psi);

% 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 t = theta
 for p = psi
	X = [X x(R,t,p)];
	Y = [Y y(R,t,p)];
	
	Z = [Z z(r,p)];
 end
 end

close all
figure; hold on;
%F = plot3(X,Y,Z,'o-','MarkerSize',20)
surf(X,Y,Z);
axis equal;


U = figure; 
view(47,20);
% Uncomment below to see 'gluing the patches together
%hold on;
for j = 1:25 %00
	% Random centers of basis balls
	P = -2+3*rand(3,1);
	S = [];
	r = 2+1.1*rand(1);

	P1 = P(1)*ones(size(X));
	P2 = P(2)*ones(size(X));
	P3 = P(3)*ones(size(X));
	A  = sqrt((X-P1).^2 + (Y-P2).^2 + (Z-P3).^2);
	[NI,NJ]  = find(A < r);
	k = length(NI);
	% you are creating an 700 by 700 jatrix in low d cases
	for i = 1:k
	xs(i) = X(NI(i),NJ(i));
	ys(i) = Y(NI(i),NJ(i));
	zs(i) = Z(NI(i),NJ(i));
	i
	end

	figure(U);
	scratch
	pause(1)
	
	
end




