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

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],G[5],G[6]},#,1],1]])&,{IdentityMatrix[4]},length])//Length//Timing

s1:={-1/3};s2:={17/19};s3={1/3};s4={7/11};
i:=1;j:=1;

ball0v[1]=Prepend[{ Sqrt[(1-s1[[i]])/2] Sin[θ]Cos[ϕ],Sqrt[(1-s1[[i]])/2]Sin[θ]Sin[ϕ],Cos[θ] (1-s1[[i]])/2+(1+s1[[i]])/2},1];

ball0v[2]=Prepend[{ Sqrt[(1-s2[[j]])/2] Sin[θ]Cos[ϕ],Sqrt[(1-s2[[j]])/2]Sin[θ]Sin[ϕ],Cos[θ] (1-s2[[j]])/2+(1+s2[[j]])/2}.RotationMatrix[π, {0,((-(Sqrt[2]))/3),(2/3)}],1];
ball0v[3]=Prepend[{ Sqrt[(1-s2[[j]])/2] Sin[θ]Cos[ϕ],Sqrt[(1-s2[[j]])/2]Sin[θ]Sin[ϕ],Cos[θ] (1-s2[[j]])/2+(1+s2[[j]])/2}.RotationMatrix[π, {0,((-(Sqrt[2]))/3),(2/3)}].RotationMatrix[2π/3, {0,0,1}],1];
ball0v[4]=Prepend[{ Sqrt[(1-s2[[i]])/2] Sin[θ]Cos[ϕ],Sqrt[(1-s2[[i]])/2]Sin[θ]Sin[ϕ],Cos[θ] (1-s2[[i]])/2+(1+s2[[i]])/2}.RotationMatrix[π, {0,((-(Sqrt[2]))/3),(2/3)}].RotationMatrix[-2π/3, {0,0,1}],1];
ball0v[5]=Prepend[{ Sqrt[(1-s4[[j]])/2] Sin[θ]Cos[ϕ],Sqrt[(1-s4[[j]])/2]Sin[θ]Sin[ϕ],Cos[θ] (1-s4[[j]])/2+(1+s4[[j]])/2}.RotationMatrix[π, {0,((-(Sqrt[2]))/3),((-2)/3)}].RotationMatrix[π, {0,1,0}],1];
ball0v[6]=Prepend[{ Sqrt[(1-s4[[j]])/2] Sin[θ]Cos[ϕ],Sqrt[(1-s4[[j]])/2]Sin[θ]Sin[ϕ],Cos[θ] (1-s4[[j]])/2+(1+s4[[j]])/2}.RotationMatrix[π, {0,((-(Sqrt[2]))/3),((-2)/3)}].RotationMatrix[π, {0,1,0}].RotationMatrix[2π/3, {0,0,1}],1];
ball0v[7]=Prepend[{ Sqrt[(1-s4[[j]])/2] Sin[θ]Cos[ϕ],Sqrt[(1-s4[[j]])/2]Sin[θ]Sin[ϕ],Cos[θ] (1-s4[[j]])/2+(1+s4[[j]])/2}.RotationMatrix[π, {0,((-(Sqrt[2]))/3),((-2)/3)}].RotationMatrix[π, {0,1,0}].RotationMatrix[-2π/3, {0,0,1}],1];
ball0v[8]=Prepend[{ Sqrt[(1-s3[[j]])/2] Sin[θ]Cos[ϕ],Sqrt[(1-s3[[j]])/2]Sin[θ]Sin[ϕ],Cos[θ] (1-s3[[j]])/2+(1+s3[[j]])/2}.RotationMatrix[π, {0,1,0}],1];

Balls=Drop[Flatten[Table[(matrices.ball0v[j])[[i]]/(matrices.ball0v[j])[[i,1]],{j,1,8},{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["Cube0.pdf",graphics]