(*
    (c) Alex Kontorovich 2017

    Companion to ``Letter to Baragar on Crystallographic Sphere Packings''
*)

n=3; (* this is the dimension of R^n, packed by (n-1) spheres *)

Q=-IdentityMatrix[n+2];Q[[1,1]]=Q[[2,2]]=0;Q[[1,2]]=Q[[2,1]]=1/2; (* This is the matrix of the form Q *)

(* Code to display a sphere (or plane) given its inverse coordinates -- Thanks to Kei Nakamura for sharing this *)
MR=100;
showSph[w_]:=If[w[[2]]≠0,{Sphere[Take[w,-3]/w[[2]],Abs[1/w[[2]]]]},Module[{n=Take[w,-3],h=w[[1]]/2},If[n[[1]]==0,Module[{p1={1,0,0},p2={0,n[[3]],-n[[2]]}},{Polygon[{h n+1.75 MR p1,h n+1.75 MR p2,h n-1.75 MR p1,h n-1.75 MR p2}]}],If[n[[2]]==0,Module[{p1={0,1,0},p2={-n[[3]],0,n[[1]]}},{Polygon[{h n+1.75 MR p1,h n+1.75 MR p2,h n-1.75 MR p1,h n-1.75 MR p2}]}],If[n[[3]]==0,Module[{p1={0,0,1},p2={n[[2]],-n[[1]],0}},{Polygon[{h n+1.75 MR p1,h n+1.75 MR p2,h n-1.75 MR p1,h n-1.75 MR p2}]}],Module[{p1={n[[2]],-n[[1]],0},p2={0,n[[3]],-n[[2]]}},{Polygon[{h n+1.75 MR p1,h n+1.75 MR p2,h n-1.75 MR p1,h n-1.75 MR p2}]}]]]]]]

(* Now we determine the walls of the Coxeter diagram *)
v1={0,0,0,0,-1}; (* chosen arbitrarily to be the xy-plane*)
Show[Graphics3D[{Opacity[.2],showSph[v1]}],PlotRange->(range={{-2,2},{-2,2},{-2,2}})]

(* Next walls (4), (5), and (6) *)
v4={0,0,0,1,0};
v5={0,0,((Sqrt[3])/2),-(1/2),0};
v6={2,0,-1,0,0};
Show[
Graphics3D[{Opacity[.2],showSph[v1]}],
Graphics3D[{Opacity[.2],showSph[v4]}],
Graphics3D[{Opacity[.2],showSph[v5]}],
Graphics3D[{Opacity[.2],showSph[v6]}],PlotRange->range]

(* general inversive coordinates *)
vec=Table[Subscript[bx,j],{j,1,n}];vec=Join[{btil,b},vec];
(* Now solve for v3 knowing its products with 1, 4, 5, and 6 *)
v3=vec/.Solve[
vec.Q.v1==0 &&
vec.Q.v4==1/2 &&
vec.Q.v5==0 &&
vec.Q.v6==0 &&
vec.Q.vec==-1
,vec
][[1]]
v3//TeXForm
Show[
Graphics3D[{Opacity[.2],showSph[v1]}],
Graphics3D[{Opacity[.2],showSph[v3]}],
Graphics3D[{Opacity[.2],showSph[v4]}],
Graphics3D[{Opacity[.2],showSph[v5]}],
Graphics3D[{Opacity[.2],showSph[v6]}],PlotRange->2range]

(* Now find v2 using v3 *)
v2=
{btil/.Solve[{btil,0,0,0,1}.Q.v3==1/Sqrt[2],btil][[1]]
,0,0,0,1}

(* collect inversive coordinates together *)
V={v1,v2,v3,v4,v5,v6};
V//MatrixForm

Show[
Graphics3D[{Opacity[.4],showSph/@V}],
PlotRange->2.5range,ViewPoint -> 
   {-2.855269116857883, 
    -0.6284834705971442, 
    1.7036568896065218}, 
  ViewVertical -> 
   {-0.14306063186958326, 
    0.12800329019222398, 
    0.9814014536921377}]

(* From the Structure Theorem, choose the cluster to be v1, cocluster is the rest *)
cluster={v1};
cocluster={v2,v3,v4,v5,v6};

(* This is the formula for reflection through v *)
vecToMat[v_]:=IdentityMatrix[n+2]+2(Transpose[{(Q.v)}].{v});

(* Compute reflection matrices *)
Rmats=vecToMat/@cocluster;
(* (* For exporting to TeX *)
For[i=1,i≤Length[Rmats],i++,
Print[TeXForm[Rmats[[i]]]]]
*)
MatrixForm/@Rmats

(* Code for generating sphere packing *)
morecircs[from_,level_]:=If[level<levelmax, (* this controls the level *)
Module[{j},For[j=1,j≤Length[Rmats],j++,If[
from≠j,(* the maps are all involutions *)
w=Simplify[w.Rmats[[j]]];(* Apply jth action *)
If[(w[[2]] <sizemax  && !MemberQ[pic,w]), (* if the new bend is not too big, and we haven't seen this sphere already *)
AppendTo[pic,w]; (* add it to the list *)

morecircs[j,level+1]; (* Iterate *)
];
w=w.Rmats[[j]]; (* undo the jth involution *)
]
]
]
];

sizemax=70;
levelmax=50;

(* we take orbits one at a time -- though in this case the cluster only has one sphere, so there is only one orbit *)
AllPics=Table[{},{k,1,Length[cluster]}];
(* initialize pics *)

(* loop over cluster *)
For[k=1,k≤Length[cluster],k++,

w=cluster[[k]]; (* set w to be the cluster vector and initialize "pic" *)
pic={w};

Print[Timing[
morecircs[0,1]; (* do the loop *)
]];

AllPics[[k]]=Union[pic];
Print["Done with k=",k," out of ",Length[cluster],". Produced ",Length[AllPics[[k]]], " spheres."];
];

(* choose which spheres to show *)
show=Select[AllPics[[1]],#[[2]]<100&];
Length[show]

(* show spheres with cocluster *)
Show[Graphics3D[{Blue,Thick,Opacity[.5],showSph/@(show)}],Graphics3D[{Red,Thick,Opacity[0.2],showSph/@cocluster}],PlotRange->3range,ViewPoint -> 
   {-2.855269116857883, 
    -0.6284834705971442, 
    1.7036568896065218}, 
  ViewVertical -> 
   {-0.14306063186958326, 
    0.12800329019222398, 
    0.9814014536921377}]

(* What are the bends? Turns out we need to rescale this image by Sqrt[3/2]; then they are all integers*)
Transpose[show][[2]]Sqrt[3/2]
%//Union