N=100

k=1:N;

k=fliplr(k);

hk=-1+(2.*(k-1))./(N-1);

Thetak=acos(hk);

Thetak=fliplr(Thetak);

Phik(1)=0;

for n=2:N-1;

Phik(n)=mod((Phik(n-1)+(3.6./sqrt(N)).*1/sqrt(1-hk(n).^2)),2*pi);

end

Phik(1)=0;

Phik(N)=2*pi;

R=ones(size(Thetak));

x = cos(Phik).*sin(Thetak);

y = sin(Phik).*sin(Thetak);

z = cos(Thetak);

x = x';

y = y';

z = z';

matrix = [x y z];

[V,C] = voronoin(matrix);

figure('Color','w')

plot3(x,y,z,'Marker','.','MarkerEdgeColor','b','MarkerSize',10, 'LineStyle', 'none')

axis equal

view([90 90 90])

k=N;

for k=1:length(C)

if all(C{k}~=1)

VertCell = V(C{k},:);

KVert = convhulln(VertCell);

patch('Vertices',VertCell,'Faces',KVert,'FaceColor','b','FaceAlpha',0.9)

end

end