function X = geodesicSphere(maxgen)
% Make a geodesic sphere
%
% X = geodesicSphere(maxgen)
%
% Inputs: number of subdivisions
% Outputs: Nx3 matrix of surface points

% Icosahedron start
phi=(1+sqrt(5))/2;
X=[0 1 phi;
    0 -1 phi;
    0 1 -phi;
    0 -1 -phi;
    1 phi 0;
    -1 phi 0;
    1 -phi 0;
    -1 phi 0;
    -1 -phi 0;
    phi 0 1;
    -phi 0 1;
    phi 0 -1;
    -phi 0 -1];

for i=1:size(X,1)
    X(i,:)=X(i,:)/norm(X(i,:));
end

for gen=1:maxgen
    tri=convhulln(X);
    lin=[tri(:,1:2); tri(:,2:3); [tri(:,3) tri(:,1)]];
    used=sparse(size(X,1),size(X,1));
    X2=[];
    for i=1:size(lin,1)
        if (used(lin(i,1),lin(i,2))==0)
            NX=(X(lin(i,1),:)+X(lin(i,2),:))/2;
            NX=NX/norm(NX);
            X2 = [X2; NX];
            used(lin(i,1),lin(i,2))=1;
        end
    end
    X=[X ; X2];
end