function MakeSpiral()
close all;
%% generate data
% data are arranged near a spiral
N = 800;
t = linspace(0,2*pi,N)';
z = cos(t);
r = 1 - 0.5*sin(t);
phi = 10*t;
x = r.*cos(phi);
y = r.*sin(phi);
sigma = 0.02; % size of the noise
X = [x,y,z] + sigma*randn(N,3);
%% display data
figure;
hold on; grid;
col = parula(N/2);
col = [col;col(end : - 1 : 1,:)];
for k = 1 : N
    plot3(X(k,1),X(k,2),X(k,3),'.','Markersize',15,'color',col(k,:));
end
view(3);
%% compute distance squared matrix
for i = 1 : N
    d(:,i) = sum((X - ones(N,1)*X(i,:)).^2,2);
end
%% compute density of points
rad = 0.5;
% density is the number of points per ball of radius rad
rad2 = rad^2;
q = zeros(N,1);
for i = 1 : N
    ind = find(d(:,i) < rad2);
    q(i) = length(ind);
end
figure;
hold on; grid;
for k = 1 : N
    plot(k,q(k),'.','Markersize',15,'color',col(k,:));
end
%% record the data
save('SpiralData.mat','X','col','d','q');
end