G[1]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,-1}};
G[2]={{-(3/2),0,-1,(1/2)},{0,-1,0,0},{1,0,1,-1},{-(1/2),0,-1,-(1/2)}};
G[3]={{-(3/2),-((Sqrt[3])/2),(1/2),(1/2)},{((Sqrt[3])/2),(1/2),-((Sqrt[3])/2),-((Sqrt[3])/2)},{-(1/2),-((Sqrt[3])/2),-(1/2),(1/2)},{-(1/2),-((Sqrt[3])/2),(1/2),-(1/2)}};
G[4]={{-(3/2),((Sqrt[3])/2),(1/2),(1/2)},{-((Sqrt[3])/2),(1/2),((Sqrt[3])/2),((Sqrt[3])/2)},{-(1/2),((Sqrt[3])/2),-(1/2),(1/2)},{-(1/2),((Sqrt[3])/2),(1/2),-(1/2)}};

length=2; (*This parameter sets the number of reflections with the group elements*)
(matrices=
Nest[(
Union[#~Join~Flatten[Outer[Dot,{G[1],G[2],G[3],G[4]},#,1],1]])&,{IdentityMatrix[4]},length])//Length//Timing

ball0v[1]={1,Sqrt[(1)/2] Sin[θ]Cos[ϕ],Sqrt[(1)/2]Sin[θ]Sin[ϕ],Cos[θ] (1)/2+(1)/2};
ball0v[2]={1,Sqrt[(1-3/5)/2] Sin[θ]Cos[ϕ],(1-3/5)/2Sin[θ]Sin[ϕ]+(1+3/5)/2,Sqrt[(1-3/5)/2]Cos[θ]};
ball0v[3]={1,-((Cos[ϕ] Sin[θ])/(2 (Sqrt[5])))+(1/2) (Sqrt[3]) ((4/5)+(1/5) Sin[θ] Sin[ϕ]),-(1/2) (Sqrt[(3/5)]) Cos[ϕ] Sin[θ]+(1/2) (-(4/5)-(1/5) Sin[θ] Sin[ϕ]),((Cos[θ])/(Sqrt[5]))};
ball0v[4]={1,-((Cos[ϕ] Sin[θ])/(2 (Sqrt[5])))-(1/2) (Sqrt[3]) ((4/5)+(1/5) Sin[θ] Sin[ϕ]),(1/2) (Sqrt[(3/5)]) Cos[ϕ] Sin[θ]+(1/2) (-(4/5)-(1/5) Sin[θ] Sin[ϕ]),((Cos[θ])/(Sqrt[5]))};

Balls=Drop[Flatten[Table[(matrices.G[j].ball0v[j])[[i]]/(matrices.G[j].ball0v[j])[[i,1]],{j,1,4},{i,1,Length[matrices]}],1],None,1];

Balls=N[Balls];

Length[Balls]

graphics=ParametricPlot3D[Balls,{ϕ,  0, 2Pi},{θ,0, 2 Pi},PlotRange->1,Mesh->None,Boxed->False,PlotStyle->Yellow,Axes->None]

Export["KSz1.pdf",graphics]