F[x_]:=1/2 I (PolyLog[2,-I E^x]-PolyLog[2,I E^x])

G3[a_,b_]:=F[a-b]+F[b-a]-π/2 (a+b);

EnergySetup[n_,FacesCW_]:=Module[{i,rhoV,rhoF,include,LineFaces,f,Faces,NM,nonlinefaces,radii,RADII},
Faces=Sort/@FacesCW;
f=Length[Faces];
LineFaces={};rhoV=Table[Subscript[rV, i],{i,1,n-2}];rhoF=Table[Subscript[rF, i],{i,1,f-2}];
include = {};
nonlinefaces={};
For[i=1,i≤f,i++,
If[!(Faces[[i,-2]]==n-1&&Faces[[i,-1]]==n),
(* then *) AppendTo[include,i];AppendTo[nonlinefaces,FacesCW[[i]]],
(* else *) AppendTo[LineFaces,i]]];

NM=NMinimize[{(∑_(k=1)^(f-2))(∑_(i=1)^(Length[Faces[[include[[k]]]]]))If[Faces[[include[[k]],i]]<n-1,G3[(Subscript[rV, Faces[[include[[k]],i]]]),(Subscript[rF, k])],0]+(∑_(i=1)^(n-2))(Subscript[rV, i])π If[MemberQ[Flatten@Table[Faces[[LineFaces[[m]]]],{m,1,2}],i],(* then *) 1,(* else *) 2]+(∑_(j=1)^(f-2))(Subscript[rF, j])π If[MemberQ[Faces[[include[[j]]]],n-1]||MemberQ[Faces[[include[[j]]]],n],(* then *) 1,(* else *) 2],
(∑_(i=1)^(n-2))(Subscript[rV, i])+(∑_(j=1)^(f-2))(Subscript[rF, j])==0},Join[rhoV,rhoF]];
RADII=Exp[Re[Join[rhoV,rhoF]/.(NM[[2]])]];
RADII/=RADII[[1]];
radii=<||>;
For[i=1,i≤n-2,i++,
radii[i]=RADII[[i]]
];
For[i=1,i≤Length[nonlinefaces],i++,
radii[nonlinefaces[[i]]]=RADII[[n-2+i]]
];
{radii,LineFaces}
]

Boundary[n_,FacesCW_,LineFaces_,radii_,ϵ_:0]:=Module[{i,j,faces,vertices,lastvisited,d,h,H,breaker,unvisited,f,endcircle,ComputeEndcircle,ComputeEndcircle2,delete},
unvisited =FacesCW;(*Keeps track of the faces we haven't visited*)
f=Length[FacesCW];

faces=<||>; (* centers of face circles (or normals, for lines) *)
For[i=1,i≤f,i++,faces[FacesCW[[i]]]={-∞,-∞}];
vertices=<||>; (* centers of vertex circles (or normals, for lines) *)
For[i=1,i≤n-2,i++,vertices[i]={-∞,-∞}];
(* Radii is an assoc from v->rad and facelist->rad *)

(*Vertex circles going vertically upwards along first FaceLine*)
h=0; (* height along the left line *)
For[i=1,i≤Length[FacesCW[[LineFaces[[1]]]]]-2(*Number of vertex circles that lie along first LineFace*),i++,
vertices[FacesCW[[LineFaces[[1]],i]]]={0,h+radii[FacesCW[[LineFaces[[1]],i]]]};(*Set the center of this vertex circle*)
h+=2radii[[FacesCW[[LineFaces[[1]],i]]]]];(*Keeps track of our cumulative height as we climb up the first FaceLine*)

unvisited=Delete[unvisited,Position[unvisited,FacesCW[[LineFaces[[1]]]]]];

d=0;(*Horizontal distance*)

ComputeEndcircle[prev_,face_]:=
(*Returns which vertex circle we calculated last based on previous vertex circle and face*)
If[face[[Mod[Position[face,n-1][[1,1]]-2,Length[face]]+1]]==prev,
(*then*)face[[Mod[Position[face,n-1][[1,1]],Length[face]]+1]],
(*else*) face[[Mod[Position[face,n-1][[1,1]]-2,Length[face]]+1]]];

endcircle=ComputeEndcircle[n,FacesCW[[LineFaces[[1]]]]];(*Previous vertex circle is bottom line, face is first FaceLine*)

(*Face circles going horizontally across upper vertex line*)
breaker=True;
For[j=1,j≤f ,j++,
For[i=1,i≤Length[unvisited],i++,
If[MemberQ[unvisited[[i]],n-1], (*determine if this face-circle contains n-1, so it is on the upper vertex line*)
If[MemberQ[unvisited[[i]],ComputeEndcircle[endcircle,unvisited[[i]]]],
(*determine if this face-circle also contains the vertex right before n-1 to ensure we are going horizontally across in the right order*)
If[unvisited[[i]]==FacesCW[[LineFaces[[2]]]],
(* If we've reached the second faceline, then do nothing--find the other face circle containing endcircle*)
,
(* else if the face circle is not a line*)
faces[unvisited[[i]]]={d+radii[unvisited[[i]]],h}; (*Set the center of this face circle*)
d+=2radii[unvisited[[i]]]; (* update how far along the top line we are *)
endcircle=ComputeEndcircle[endcircle,unvisited[[i]]];(*Update endcircle*)
unvisited=Delete[unvisited,i];] (* remove that face from the unvisited list *)
]
]
]
];
(*Vertex circles along the second FaceLine, climbing vertically upwards*)
h=0;
For[i=1,i≤Length[FacesCW[[LineFaces[[2]]]]]-2,i++,
vertices[FacesCW[[LineFaces[[2]],i]]]={d,h+radii[FacesCW[[LineFaces[[2]],i]]]};h+=2radii[[FacesCW[[LineFaces[[2]],i]]]]];
H=h;
h=0;(*Reset height*)

ComputeEndcircle2[prev_,face_]:=
If[face[[Mod[Position[face,n][[1,1]]-2,Length[face]]+1]]==prev,
(*then*)face[[Mod[Position[face,n][[1,1]],Length[face]]+1]],
(*else*) face[[Mod[Position[face,n][[1,1]]-2,Length[face]]+1]]];

endcircle=ComputeEndcircle2[n-1,FacesCW[[LineFaces[[1]]]]];
(*Face circles along lower vertex line, moving horizontally from origin to second FaceLine*)
d=0;(*Reset distance*)
breaker=True;
For[j=1,j≤f,j++,
For[i=1,i≤Length[unvisited],i++,
If[MemberQ[unvisited[[i]],n], (* determine if this face-circle contains n *)
If[MemberQ[unvisited[[i]],endcircle(* ComputeEndcircle2[endcircle,unvisited[[i]]] *)],(* determine if this face-cirlce also contains the vertex right before n *)
If[unvisited[[i]]==FacesCW[[LineFaces[[2]]]],
(* If we've reached the second faceline, then do nothing*)
,
(* else *)
faces[unvisited[[i]]]={d+radii[unvisited[[i]]],h}; (* then we'll make this face's center be where we know it is *)
d+=2radii[unvisited[[i]]]; (* update how far along the bottom line we are *)
endcircle=ComputeEndcircle2[endcircle,unvisited[[i]]];
unvisited=Delete[unvisited,i];] (* remove that face from the unvisited list *)
]
]
]
];
(*Set the normals of the vertex lines and face lines*)
faces[FacesCW[[LineFaces[[2]]]]]={1,0};
faces[FacesCW[[LineFaces[[1]]]]]={-1,0};
vertices[n]={0,-1};
vertices[n-1]={0,1};

If[ThrowOrderError[faces,vertices,radii,d,H,n,FacesCW[[LineFaces]],0],
Print["flipping"];
Print[faces,vertices,d,H,FacesCW[[LineFaces]]];
faces[FacesCW[[LineFaces[[1]]]]]*=-1;
faces[FacesCW[[LineFaces[[2]]]]]*=-1;
For[i=1,i≤n,i++,
If[i==n-1||i==n,
,
If[Abs[vertices[i][[1]]-d]≤ϵ,
vertices[i]={0,vertices[i][[2]]},
If[Abs[vertices[i][[1]]-0]≤ϵ,
vertices[i]={d,vertices[i][[2]]}]];
]]];

(*Return the assoc's of the faces to their centers and vertices to their centers*)
{vertices,faces,d,H}
]

FindRest[{x_,y_},r_,face_,vertices_,radii_,lines_]:=(* face is all the vertices that the face intersects *)
(*{x_,y_},r: known center of opposite kind and its radius
face: all the circles of its kind that pass thru that opposite known circle
vertices: all the circles of its kind that are tangent to the one we'e looking for*)
Module[{verts=vertices,i,j,ϕ,θ,ω,face2=face,b,b2=True},

If[verts[face2[[1]]]=={-∞,-∞}||MemberQ[lines,face2[[1]]],
For[i=2,i≤Length[face2]&&b2,i++,
If[verts[face2[[i]]]!={-∞,-∞}&&!MemberQ[lines,face2[[i]]],
face2=Join[face2[[i;;]],face2[[;;i-1]]];
b2=False;
]]];

(* Print["face2 is now ",face2]; *)

b=True;
For[i=2,i≤Length[face2]&&b,i++,
If[MemberQ[lines,face2[[i]]],
AppendTo[face2,face2[[1]]];
For[j=Length[face2]-1,j>i,j--,
If[verts[face2[[j]]]=={-∞,-∞},
ϕ=ArcTan[verts[face2[[j+1]]][[1]]-x,verts[face2[[j+1]]][[2]]-y];
θ=ArcTan[radii[face2[[j+1]]]/r];(* j+1 refers to previously computed / known *)
ω=ArcTan[radii[face2[[j]]]/r];
verts[face2[[j]]]={x+Sqrt[r^2+radii[face2[[j]]]^2]Cos[ϕ-θ-ω],y+Sqrt[r^2+radii[face2[[j]]]^2]Sin[ϕ-θ-ω]}
]
];b=False;Break,
If[verts[face2[[i]]]=={-∞,-∞},
ϕ=ArcTan[verts[face2[[i-1]]][[1]]-x,verts[face2[[i-1]]][[2]]-y];
θ=ArcTan[radii[face2[[i-1]]]/r];(* i-1 refers to previously computed / known *)
ω=ArcTan[radii[face2[[i]]]/r];
verts[face2[[i]]]={x+Sqrt[r^2+radii[face2[[i]]]^2]Cos[ϕ+θ+ω],y+Sqrt[r^2+radii[face2[[i]]]^2]Sin[ϕ+θ+ω]}];
];];
Simplify[verts]
]

Cull[assoc_]:=Module[{vis={},unv={},i},For[i=1,i≤Length[Keys[assoc]],i++,If[assoc[Keys[assoc][[i]]]=={-∞,-∞},AppendTo[unv,Keys[assoc][[i]]],AppendTo[vis,Keys[assoc][[i]]]]];{vis,unv}]

Interior[n_,FacesCW_,LineFaces_,radii_,vertices_,faces_,d_,h_]:=
Module[{Vunv={},Funv={},Vvis={},Fvis={},i,j,f,verts=vertices,verts2,fcs=faces,fcs2,fi,vi,x,y,k,b},
f=Length[FacesCW];
For[i=1,i≤n,i++,
If[vertices[i]=={-∞,-∞},AppendTo[Vunv,i],AppendTo[Vvis,i]]];
For[i=1,i≤f,i++,
If[faces[FacesCW[[i]]]=={-∞,-∞},AppendTo[Funv,FacesCW[[i]]],AppendTo[Fvis,FacesCW[[i]]]]];

(* Print["initial things are ",Vunv,Vvis,Funv,Fvis]; *)

For[j=1,j≤n+Length[FacesCW]&&(Length[Vunv]+Length[Funv]>0),j++,

For[i=1,i≤Length[Fvis],i++,(* for a known face circle *)
If[Length[Intersection[Vvis,Fvis[[i]]]]*Length[Intersection[Vunv,Fvis[[i]]]]>0,
fi=Fvis[[i]];
{x,y}=fcs[fi];
verts2=FindRest[{x,y},radii[fi],fi,verts,radii,{n-1,n}];
If[ThrowOrderError[fcs,verts2,radii,d,h,n,{FacesCW[[LineFaces[[1]]]],FacesCW[[LineFaces[[2]]]]},0],verts=FindRest[{x,y},radii[fi],Reverse[fi],verts,radii,{n-1,n}],
verts=verts2];
(* 
Print[Render[radii,fcs,verts,FacesCW[[LineFaces]],n,d,h]];
Print[N[fcs],N[verts]]
*)
];
{Vvis,Vunv}=Cull[verts];
];

For[i=1,i≤Length[Vvis],i++, (* for a known vertex circle *)
If[b=False;
For[k=1,k≤Length[Fvis],k++,b=b||If[MemberQ[Fvis[[k]],Vvis[[i]]],True,False]];b,
If[b=False;For[k=1,k≤Length[Funv],k++,b=b||If[MemberQ[Funv[[k]],Vvis[[i]]],True,False]];b,
(* then *)
vi=Vvis[[i]];
{x,y}=verts[vi];
fcs2=FindRest[{x,y},radii[vi],FacesList[vi,FacesCW],fcs,radii,FacesCW[[LineFaces]]];
If[ThrowOrderError[fcs2,verts,radii,d,h,n,FacesCW[[LineFaces]],0],
(* 
Print["it thinks that this is wrong: ",N[fcs2]];
Print["it would have rendered ",Render[radii,fcs2,verts,FacesCW[[LineFaces]],n,d,h]]; 
*)
fcs=FindRest[{x,y},radii[vi],
(* Print["it thought the bad stuff was in ",fcs2]; *)
Reverse[FacesList[vi,FacesCW]],fcs,radii,FacesCW[[LineFaces]]],
fcs=fcs2];
(* 
Print[Render[radii,fcs,verts,FacesCW[[LineFaces]],n,d,h]];
Print[N[fcs],N[verts]]
*)
];
{Fvis,Funv}=Cull[fcs];
]
];
];
{verts,fcs}
]

ThrowOrderError[fcs2_,verts_,radii_,d_,h_,n_,fls_,ϵ_](* fcs2 is proposed fcs association. fls is facelines *):=
Module[{i,j,b,b2,K},
K={};
For[i=1,i≤Length[Keys[fcs2]],i++,
If[!MemberQ[fls,Keys[fcs2][[i]]],AppendTo[K,Keys[fcs2][[i]]]]
];

b=False; (* "breaker" *)
For[i=1,i≤Length[K]&&!b,i++,
If[b2=False;
-∞<fcs2[K[[i]]][[1]]<0
||-∞<fcs2[K[[i]]][[2]]<0
||fcs2[K[[i]]][[1]]>d
||fcs2[K[[i]]][[2]]>h
||(For[j=1,j≤Length[K[[i]]],j++,
If[K[[i,j]]!=n-1&&K[[i,j]]!=n&&verts[K[[i,j]]]!={-∞,-∞}&&fcs2[K[[i]]]!={-∞,-∞},
If[EuclideanDistance[verts[K[[i,j]]],fcs2[K[[i]]]]>=(1-ϵ)(radii[K[[i,j]]]+radii[K[[i]]]),b2=True;Break]]];b2),
b=True,

For[j=1,j≤Length[K]&&!b,j++,
If[i!=j&&fcs2[K[[i]]]!={-∞,-∞}&&fcs2[K[[j]]]!={-∞,-∞}&&
EuclideanDistance[fcs2[K[[i]]],fcs2[K[[j]]]]<(1+ϵ)(radii[fcs2[K[[i]]]]+radii[fcs2[K[[j]]]])
,(* then, bad *) b=True]]];
];


For[i=1,i≤Length[fls[[1]]]&&!b,i++,
If[verts[fls[[1,i]]][[1]]!=d/2+fcs2[fls[[1]]][[1]] d/2&&fls[[1,i]]<n-1,b=True]
];
For[i=1,i≤Length[fls[[2]]]&&!b,i++,
If[verts[fls[[2,i]]][[1]]!=d/2+fcs2[fls[[2]]][[1]] d/2&&fls[[2,i]]<n-1,b=True]
];
b
]

FacesList[v_,Faces_]:=
Module[{ContainV={},i,L=2,p=3,temp},
For[i=1,i≤Length[Faces],i++,
If[MemberQ[Faces[[i]],v],AppendTo[ContainV,Faces[[i]]]]];
While[L<Length[ContainV],
While[p≤Length[ContainV],
If[Length[Intersection[ContainV[[L-1]],ContainV[[p]]]]==2,
temp=ContainV[[p]];ContainV[[p]]=ContainV[[L]];ContainV[[L]]=temp;p=1+Length[ContainV],
p++]
];
L++;
p=L+1;
];
ContainV]

Packing[n_,FacesCW_,ϵ_:0]:=
Module[{radii,LineFaces,vertices,faces,d,h},
{radii,LineFaces}=EnergySetup[n,FacesCW];
{vertices,faces,d,h}=Boundary[n,FacesCW,LineFaces,radii,ϵ];
{vertices,faces}=Interior[n,FacesCW,LineFaces,radii,vertices,faces,d,h];
{vertices,faces,radii,d,h,LineFaces}];

Render[radii_,fcs_,verts_,fls_,n_,d_,h_]:=Graphics[{Blue,Thick,
Table[If[i==n-1,Line[{{-Max[Values[radii]],h},{d+Max[Values[radii]],h}}],If[i==n,Line[{{-Max[Values[radii]],0},{d+Max[Values[radii]],0}}],
Circle[verts[Keys[verts][[i]]],radii[Keys[verts][[i]]]]]],{i,1,n}],
Red,Dashed,
Table[If[Keys[fcs][[i]]==fls[[1]],Line[{{0,-Max[Values[radii]]},{0,h+Max[Values[radii]]}}],
If[Keys[fcs][[i]]==fls[[2]],Line[{{d,-Max[Values[radii]]},{d,h+Max[Values[radii]]}}],
Circle[fcs[Keys[fcs][[i]]],radii[Keys[fcs][[i]]]]]],{i,1,Length[Keys[fcs]]}]
}]

R[vhat_,d_:2]:=IdentityMatrix[d+2]+2𝒬[d].Transpose[{vhat}].{vhat}

Everything[n_,FacesCW_,ϵ_:0]:=Module[{vertices,faces,f=Length[FacesCW],radii,d,h,LineFaces,Send1,Send2,Vnye,bends={}},
Send1[x_]:=SendToZero[x,ϵ];
Send2[list_]:=Send1/@list;
{vertices,faces,radii,d,h,LineFaces}=Packing[n,FacesCW,ϵ];
Vnye=MatrixForm[RootApproximant[SuperCluster[vertices,faces,radii,FacesCW[[LineFaces]],d,h,n]]];
bends=Bend[Vnye[[1]],n]/@Range@f;
Column[{
Render[radii,faces,vertices,FacesCW[[LineFaces]],n,d,h],
Vnye,
bends,
MatrixForm[RootApproximant[
Send2/@Gram[SuperCluster[vertices,faces,radii,FacesCW[[LineFaces]],d,h,n]]
]]
}]]

InvCoords[z_,r_,d_:0]:=If[r==∞,{2d,0,z},{Norm[z]^2/r-r,1/r,z/r}]

Gram[SC_,d_:2]:=Simplify[SC.𝒬[d].Transpose[SC]];
SuperCluster[vertices_,faces_,radii_,fls_,d_,h_,n_]:=
Module[{sc={},i},
For[i=1,i≤n-2,i++,AppendTo[sc,Flatten[InvCoords[vertices[i],radii[i]]]]];
AppendTo[sc,Flatten[InvCoords[vertices[n-1],∞,h/2+vertices[n-1][[2]] h/2]]];
AppendTo[sc,Flatten[InvCoords[vertices[n],∞,h/2+vertices[n][[2]] h/2]]];
For[i=1,i≤Length[Keys[faces]],i++,
If[!MemberQ[fls,Keys[faces][[i]]],
AppendTo[sc,Flatten[InvCoords[faces[Keys[faces][[i]]],radii[Keys[faces][[i]]]]]]
]];
AppendTo[sc,Flatten[InvCoords[faces[fls[[1]]],∞,d/2+faces[fls[[1]]][[1]] d/2]]];
AppendTo[sc,Flatten[InvCoords[faces[fls[[2]]],∞,d/2+faces[fls[[2]]][[1]] d/2]]];
sc
]

SendToZero[x_,ϵ_:0]:=If[Abs[x]<ϵ,0,x];

𝒬[d_:2]:=Prepend[Prepend[Function[l,-Flatten[Prepend[l,{0,0}]]]/@IdentityMatrix[d],Flatten[Prepend[ConstantArray[0,d],{1/2,0}]]],Flatten[Prepend[ConstantArray[0,d],{0,1/2}]]]

BMat[n_]:=Table[b[i,j],{i,1,n},{j,1,n}]

Bend[V_,m_][i_]:=
Module[{sol,assignments,rest,j,l,others1,others2},
sol=Solve[BMat[m].(V[[;;m]])==(V[[;;m]]).R[V[[m+i]]]][[1]];
assignments=Association[sol];
rest=Complement[Flatten[BMat[m]],Keys[assignments]];
l=Length[rest];
others1=Solve[rest==Table[0,{j,1,l}]][[1]];
others2=<||>;
For[j=1,j≤l,j++,others2[rest[[j]]]=Subscript[a, j]];
others2=Normal[others2];
Column[{
MatrixForm[(BMat[m]/.assignments)/.others2],
MatrixForm[(BMat[m]/.assignments)/.others1]
}]
]

DrawCircle[{bhat_,b_,bx_,by_}]:=If[b!=0,Circle[{bx,by}/b,1/Abs[b]],InfiniteLine[{bx,by}bhat/2,{-by,bx}]]

DrawCircleText[fs_:500,fac_:1][{bhat_,b_,bx_,by_}]:=If[b!=0,{Text[Style[Simplify[Abs[b/fac]],FontSize->N[Abs[fs/b]]],{bx,by}/b],Circle[{bx,by}/b,1/Abs[b]]},InfiniteLine[{bx,by}bhat/2,{-by,bx}]]

Strip[cl_,co_,steps_,sizemax_:∞,pr_:{{-1,-1},{-1,-1}}]:=Module[{all,last,newlast,RZip,i,j,k,l,xmin=∞,xmax=-∞,ymin=∞,ymax=-∞,PLOTRANGE,SC=Union[cl,co]},
If[pr=={{-1,-1},{-1,-1}},
For[i=1,i≤Length[SC],i++,
If[SC[[i,2]]!=0,
If[(SC[[i,3]]-1)/Abs[SC[[i,2]]]<xmin,xmin=(SC[[i,3]]-1)/Abs[SC[[i,2]]]];
If[(SC[[i,3]]+1)/Abs[SC[[i,2]]]>xmax,xmax=(SC[[i,3]]+1)/Abs[SC[[i,2]]]];
If[(SC[[i,4]]-1)/Abs[SC[[i,2]]]<ymin,ymin=(SC[[i,4]]-1)/Abs[SC[[i,2]]]];
If[(SC[[i,4]]+1)/Abs[SC[[i,2]]]>ymax,ymax=(SC[[i,4]]+1)/Abs[SC[[i,2]]]];
]
];
PLOTRANGE={{xmin,xmax},{ymin,ymax}},
PLOTRANGE=pr
];
Graphics[{Thick,DrawCircle/@Orbit[co,cl,steps,sizemax],Blue,DrawCircle/@cl,Red,DrawCircle/@co},PlotRange->PLOTRANGE]
]

Strip2[supercl_,clindices_,pr_:{{-1,-1},{-1,-1}}]:=
Manipulate[Strip[supercl[[clindices]],supercl[[Complement[Range@Length[supercl],clindices]]],steps,∞,pr],{steps,0,20,1}]

Strip2nth[supercl_,clindices_,n_,pr_:{{-1,-1},{-1,-1}}]:=Strip[supercl[[clindices]],supercl[[Complement[Range@Length[supercl],clindices]]],n,∞,pr]

StripText[cl_,co_,steps_,sizemax_:∞,fs_:500,fac_:1,pr_:{{-1,-1},{-1,-1}}]:=Module[{all,last,newlast,RZip,i,j,k,l,xmin=∞,xmax=-∞,ymin=∞,ymax=-∞,PLOTRANGE,SC=Union[cl,co]},
If[pr=={{-1,-1},{-1,-1}},
For[i=1,i≤Length[SC],i++,
If[SC[[i,2]]!=0,
If[(SC[[i,3]]-1)/Abs[SC[[i,2]]]<xmin,xmin=(SC[[i,3]]-1)/Abs[SC[[i,2]]]];
If[(SC[[i,3]]+1)/Abs[SC[[i,2]]]>xmax,xmax=(SC[[i,3]]+1)/Abs[SC[[i,2]]]];
If[(SC[[i,4]]-1)/Abs[SC[[i,2]]]<ymin,ymin=(SC[[i,4]]-1)/Abs[SC[[i,2]]]];
If[(SC[[i,4]]+1)/Abs[SC[[i,2]]]>ymax,ymax=(SC[[i,4]]+1)/Abs[SC[[i,2]]]];
]
];
PLOTRANGE={{xmin,xmax},{ymin,ymax}},
PLOTRANGE=pr
];
Graphics[{Thick,DrawCircleText[fs,fac]/@Orbit[co,cl,steps,sizemax],Blue,DrawCircleText[fs,fac]/@cl,Red,DrawCircle/@co},PlotRange->PLOTRANGE]
]

Isolated[Gram_]:=Module[{iso={},i,j,b},
For[i=1,i≤Length[Gram],i++,
b=True;
For[j=1,j≤Length[Gram],j++,
If[i!=j&&0<Gram[[i,j]]<1,b=False]];
If[b,AppendTo[iso,i]]
];
iso
]

NoPoint[iso_,Gram_]:=Module[{i,j,k,np},
np={};
For[i=1,i≤Length[iso],i++,
For[j=1,j≤ Length[iso]-i,j++,
If[Gram[[iso[[i]],iso[[i+j]]]]==0,
AppendTo[np,{iso[[i]],iso[[i+j]]}]
]];];
np];

IsoClusters[iso_,np_]:=Module[{isocl,combos,i,j,L,sofarsogood},
combos=Subsets[iso,{1,Length[iso]}];
If[np=={},isocl=combos,
isocl={};
L=Length[combos];
For[i=1,i≤L,i++,
sofarsogood=true;
For[j=1,j≤Length[np],j++,
If[SubsetQ[combos[[i]],np[[j]]],
sofarsogood=false;]
];
If[sofarsogood==true,AppendTo[isocl,combos[[i]]]]
];];
isocl];

SendToInt[ϵ_:0][y_]:=With[{x=Simplify[y]},If[Abs[x-Floor[x]]≤ϵ ,Floor[x],
If[Abs[x-Ceiling[x]]≤ϵ ,Ceiling[x],x],x]];

Orbit[co_,cl_,steps_,sizemax_:∞,d_:2,ϵ_:0]:=Module[{i,j,k,orbit=cl,L},
For[i=1,i≤steps,i++,
L=Length[orbit];
For[j=1,j≤Length[co],j++,
For[k=1,k≤L,k++,
If[Simplify[(orbit[[k]].R[co[[j]],d])[[2]]]≤sizemax&&!MemberQ[orbit,(* SendToInt[ϵ] *)Simplify/@(orbit[[k]].R[co[[j]],d])],AppendTo[orbit,(* SendToInt[ϵ] *)Simplify/@(orbit[[k]].R[co[[j]],d])]]
]
]
];
orbit
]

Bend3[V_,sizemax_:∞,d_:2,RING_:Integers][iso_][i_]:=(* V is a supercluster matrix. iso is the indices of the cluster. 1≤i≤|cocluster| *)Module[{sol,assignments,rest,j,l,others1,others2,M=Length[V[[1]]],co,cl,orbit={},steps=0},
cl=V[[iso]];
co=V[[Complement[Range@Length[V],iso]]];
While[Length[orbit]<M,
steps++;
orbit=Orbit[co,cl,steps,sizemax,d];
];
cl=orbit[[;;M]];
sol=Solve[BMat[M].cl==cl.R[co[[i]],d],RING];
If[sol=={},Null,
sol=sol[[1]];
assignments=Association[sol];
rest=Complement[Flatten[BMat[M]],Keys[assignments]];
l=Length[rest];
If[l==0,MatrixForm[Simplify[BMat[M]/.assignments]],
others1=Solve[rest==Table[0,{j,1,l}]][[1]];
others2=<||>;
For[j=1,j≤l,j++,others2[rest[[j]]]=Subscript[a, j]];
others2=Normal[others2];
(* Column[{MatrixForm[Simplify[(BMat[M]/.assignments)/.others2]],MatrixForm[Simplify[(BMat[M]/.assignments)/.others1]]}]]*) 
If[Det[Simplify[(BMat[M]/.(assignments/.others1))/.others1]][[0]]==Symbol,Null,
MatrixForm[Simplify[(BMat[M]/.(assignments/.others1))/.others1]]]
]
]
]

Bend4[V_,d_:2,sizemax_:∞][iso_][i_]:=With[{B=Bend3[V,sizemax,d,Integers][iso][i]},
If[Length[B]>0,
{Bend3[V,sizemax,d,Reals][iso][i],"nonintegral"},
{B,"integral"}
]
]

ConvertOtherQF[qn_][x_]:=(*assumes xQx>0*)
Module[{n,Q,bhat,b,bz,ip,y=x},
n=Length[y]-2;
Q=Prepend[Function[l,Prepend[l,0]]/@IdentityMatrix[n+1],Prepend[ConstantArray[0,n+1],-qn]];
ip=y.Q.y;
If[ip>0,y/=Sqrt[ip]];
Simplify[Flatten[Prepend[y[[3;;]],{y[[1]](Sqrt[qn])+y[[2]],y[[1]](Sqrt[qn])-y[[2]]}]]]
]

SuppCoords1={
{2,(∑_(i=0)^2)v[i],(∑_(i=0)^3)v[i]},
{10,2v[0]+(∑_(i=0)^10)v[i],2v[0]+(∑_(i=0)^11)v[i]},
{14,3v[0]+v[1]+(∑_(i=0)^14)v[i],3v[0]+v[1]+(∑_(i=0)^15)v[i]},
{16,5v[0]+(∑_(i=1)^7)v[i]+(∑_(i=0)^16)v[i],3v[0]+(∑_(i=0)^17)v[i]},
{17,5v[0]+(∑_(i=1)^7)v[i]+(∑_(i=0)^17)v[i],5v[0]+(∑_(i=1)^7)v[i]+(∑_(i=0)^17)v[i]}
};
Coords1[n_]:=Module[{i,b=True,assignments={},vs},
vs=Table[-v[i]+v[i+1],{i,1,n-1}];
AppendTo[vs,-v[n]];
For[i=0,i≤n,i++,AppendTo[assignments,(v[i]->IdentityMatrix[n+1][[i+1]])]];
i=1;
While[b,
If[i>5,
(*then*) b=False,
(*else*) If[n<SuppCoords1[[i,1]],
(*then*) b=False,
(*else*) If[n==SuppCoords1[[i,1]],
(*then*) AppendTo[vs,SuppCoords1[[i,2]]];b=False,
(*else*) AppendTo[vs,SuppCoords1[[i,3]]];i++
]
]
]
];
vs/.assignments
]

SuppCoords2={
{3,(∑_(i=0)^3)v[i],(∑_(i=0)^4)v[i]},
{9,v[0]+(∑_(i=0)^9)v[i],v[0]+(∑_(i=0)^10)v[i]},
{11,2v[0]+2v[1]+(∑_(i=0)^11)v[i],2v[0]+2v[1]+(∑_(i=0)^12)v[i]},
{13,2v[0]+v[1]+v[2]+(∑_(i=0)^13)v[i],2v[0]+v[1]+v[2]+(∑_(i=0)^14)v[i]},
{13,3v[0]+2(∑_(i=0)^13)v[i],3v[0]+2(∑_(i=0)^13)v[i]}
};
Coords2[n_]:=Module[{i,b=True,assignments={},vs},
vs=Table[-v[i]+v[i+1],{i,1,n-1}];
AppendTo[vs,-v[n]];
AppendTo[vs,v[0]+2v[1]];
For[i=0,i≤n,i++,AppendTo[assignments,(v[i]->IdentityMatrix[n+1][[i+1]])]];

i=1;
While[b,
If[i>5,
(*then*) b=False,
(*else*) If[n<SuppCoords2[[i,1]],
(*then*) b=False,
(*else*) If[n==SuppCoords2[[i,1]],
(*then*) AppendTo[vs,SuppCoords2[[i,2]]];b=False,
(*else*) AppendTo[vs,SuppCoords2[[i,3]]];i++
]
]
]
];
If[n>=13,AppendTo[vs,SuppCoords2[[5,2]]]];
vs/.assignments
]

SuppCoords3={
{4,(∑_(i=0)^4)v[i],(∑_(i=0)^5)v[i]},
{9,5v[0]+3(∑_(i=1)^9)v[i],5v[0]+3(∑_(i=1)^9)v[i]},
{10,2(v[0]+v[1])+(∑_(i=2)^10)v[i],2(v[0]+v[1])+(∑_(i=2)^11)v[i]},
{12,3(∑_(i=0)^2)v[i]+(∑_(i=3)^12)v[i],3(∑_(i=0)^2)v[i]+(∑_(i=3)^13)v[i]},
{12,5v[0]+3(∑_(i=1)^8)v[i]+(∑_(i=9)^12)v[i],5v[0]+3(∑_(i=1)^8)v[i]+(∑_(i=9)^13)v[i]},
{13,2v[0]+(∑_(i=1)^13)v[i],Null},
{13,8v[0]+6(∑_(i=1)^3)v[i]+3(∑_(i=4)^13)v[i],8v[0]+6(∑_(i=1)^3)v[i]+3(∑_(i=4)^13)v[i]},{13,10v[0]+6(∑_(i=1)^7)v[i]+3(∑_(i=8)^13)v[i],10v[0]+6(∑_(i=1)^7)v[i]+3(∑_(i=8)^13)v[i]}
};
Coords3[n_]:=Module[{i,b=True,assignments={Null->Nothing},vs},
vs=Table[-v[i]+v[i+1],{i,1,n-1}];
AppendTo[vs,-v[n]];
AppendTo[vs,v[0]+3v[1]];
For[i=0,i≤n,i++,
AppendTo[assignments,(v[i]->IdentityMatrix[n+1][[i+1]])]
];
i=1;
While[b,
If[i>Length[SuppCoords3],
(*then*)b=False,
(*else*)If[n<SuppCoords3[[i,1]],
b=False,
If[n==SuppCoords3[[i,1]],
(*then*)AppendTo[vs,SuppCoords3[[i,2]]];i++,
(*else*)AppendTo[vs,SuppCoords3[[i,3]]];i++
]
]
]
];
vs/.assignments
]

SuppCoords7={
{2,v[0]+2v[1]+2v[2],v[0]+2v[1]+2v[2]+v[3]}
};
Coords7[n_]:=Module[{i,b=True,assignments={Null->Nothing},vs},
vs=Table[-v[i]+v[i+1],{i,1,n-1}];
AppendTo[vs,-v[n]];
AppendTo[vs,v[0]+3v[1]];
For[i=0,i≤n,i++,
AppendTo[assignments,(v[i]->IdentityMatrix[n+1][[i+1]])]
];
i=1;
While[b,
If[i>Length[SuppCoords7],
(*then*)b=False,
(*else*)If[n<SuppCoords7[[i,1]],
b=False,
If[n==SuppCoords7[[i,1]],
(*then*)AppendTo[vs,SuppCoords7[[i,2]]];i++,
(*else*)AppendTo[vs,SuppCoords7[[i,3]]];i++
]
]
]
];
vs/.assignments
]

DrawSphere[{bhat_,b_,bx_,by_,bz_}]:=With[{op=({{\[Piecewise], ({{.9, b>=0}, {.2, b<0}})}})},{Opacity[op],If[b==0,
Hyperplane[{bx,by,bz},{bx,by,bz}bhat/2],

Sphere[{bx,by,bz}/b,1/Abs[b]]]}]

Zip[l__]:=Module[{L=List[l],n,nth},n=Length[L[[1]]];
nth[k_][list_]:=list[[k]];
Table[nth[k]/@L,{k,1,n}]]

PiFrac[x_]:=Round[π/ArcCos[x]]

Coxeter[coords_,d_:2]:=Module[{gram,i,j,n,edges={},labels={},styles={}},
gram=Gram[coords,d];
n=Length[coords];
For[i=1,i≤n,i++,
For[j=1,j<i,j++,
If[gram[[i,j]]>0,
AppendTo[edges,{i,j}];
If[gram[[i,j]]>1,styles={styles,i<->j->Dashed},
If[gram[[i,j]]<1,labels={labels,i<->j->Framed[PiFrac[gram[[i,j]]]-2,Background->White]},labels={labels,i<->j->Framed[∞,Background->White]};styles={styles,i<->j->Thickness[0.01]}]
]
]
]
];
labels=Flatten[labels];
styles=Flatten[styles];
Graph[edges,EdgeLabels->labels,EdgeStyle->styles(*,GraphLayout->"PlanarEmbedding"*),VertexLabels->Table[k->k,{k,1,n}]]
]

McLeod=<|1->{{0,0,-1,0},{1,0,1,0},{0,0,0,-1},{1,0,0,1},{-1,1,0,0}},2->{{0,0,-1,0},{1,0,1,0},{0,0,0,-1},{2,0,0,1},{-1,1,0,0}},3->{{0,0,-1,0},{1,0,0,1},{0,0,1,-1},{-1,1,0,0}},5->{{0,0,-1,0},{1,0,1,0},{0,0,0,-1},{5,0,0,1},{-1,1,0,0},{2,2,1,1}},6->{{0,0,-1,0},{1,0,1,0},{0,0,0,-1},{6,0,0,1},{-1,1,0,0},{2,2,0,1}},
7->{{0,0,-1,0},{1,0,1,0},{0,0,1,-2},{7,0,-1,2},{-1,1,0,0},{1,1,0,1}},10->{{0,0,-1,0},{1,0,1,0},{0,0,0,-1},{10,0,0,1},{-1,1,0,0},{4,2,0,1},{8,6,3,2},{30,30,10,9},{40,40,20,11}},11->{{0,0,-1,0},{1,0,1,0},{0,0,1,-2},{11,0,-1,2},{-1,1,0,0},{2,1,0,1}},13->{{0,0,-1,0},{1,0,1,0},{0,0,0,-1},{13,0,0,1},{-1,1,0,0},{6,2,1,1},{4,3,0,1},{4,4,2,1},{52,39,13,12},{52,52,13,14}},14->{{0,0,-1,0},{1,0,1,0},{0,0,0,-1},{14,0,0,1},{-1,1,0,0},{6,2,0,1},{7,7,0,2},{4,4,2,1},{28,14,7,5}},15->{{0,0,-1,0},{1,0,1,0},{0,0,1,-2},{15,0,-1,2},{-1,1,0,0},{3,1,0,1},{15,15,4,7},{15,15,-4,8}},17->{{0,0,-1,0},{1,0,1,0},{0,0,0,-1},{17,0,0,1},{-1,1,0,0},{8,2,1,1},{4,4,1,1},{68,34,17,11},{19,8,0,3},{17,9,1,3},{136,68,17,23},{85,51,0,16},{204,102,0,35}},19->{{0,0,-1,0},{1,0,1,0},{0,0,1,-2},{19,0,-1,2},{-1,1,0,0},{4,1,0,1},{2,2,0,1}},21->{{0,0,-1,0},{1,0,1,0},{0,0,0,-1},{21,0,0,1},{-1,1,0,0},{10,2,1,1},{6,3,0,1},{6,4,2,1},{42,42,21,8},{14,14,3,3},{63,63,21,13}},30->{{0,0,-1,0},{1,0,1,0},{0,0,0,-1},{30,0,0,1},{-1,1,0,0},{14,2,0,1},{9,3,0,1},{5,5,0,1},{8,4,2,1},{6,6,3,1},{50,10,5,4}},33->{{0,0,-1,0},{1,0,1,0},{0,0,0,1},{33,0,0,1},{-1,1,0,0},{16,2,1,1},{6,6,3,1},{8,4,1,1},{11,3,1,1},{11,11,0,2},{99,33,0,10},{121,22,0,9},{90,18,3,7},{37,8,0,3},{264,66,0,23}},39->{{0,0,-1,0},{1,0,1,0},{0,0,1,-2},{39,0,-1,2},{-1,1,0,0},{9,1,0,1},{3,3,1,1},{12,3,-1,2},{26,13,-3,6},{39,13,3,7}}|>;

McLeodNorm[m_,{x_,y_,z_,w_}]:=If[Mod[m,4]==3,-2x y+2(z^2)+2z w+((m+1)/2)(w^2),-2x y+2(z^2)+2m (w^2)]

ConvertMcLeod[m_][initial_]:=Module[{coords,ip,i,j,newinitial,take},
(*5th entry of initial is (e,e)*)
newinitial=initial;
For[i=1,i≤Length[newinitial],i++,
ip=McLeodNorm[m,newinitial[[i]]];
If[ip≠2,
(*Need to normalize*)
For[j=1,j≤4,j++,
newinitial[[i,j]]=newinitial[[i,j]]/Sqrt[ip/2]];
];
If[Mod[m,4]≠3,
newinitial[[i,4]]=newinitial[[i,4]]*Sqrt[m];
,
(*Else, m%4==3, use different F*)
newinitial[[i,3]]=newinitial[[i,3]]+newinitial[[i,4]]/2;
newinitial[[i,4]]=newinitial[[i,4]]*Sqrt[m]/2;
If[newinitial[[i,2]]==0,(*If line*)
If[Sqrt[newinitial[[i,3]]^2+newinitial[[i,4]]^2]≠1,
(*if line's last coords aren't normal vector, normalize*)
For[j=1,j≤4,j++,
newinitial[[i,j]]*=1/Sqrt[newinitial[[i,3]]^2+newinitial[[i,4]]^2]];
]]]];
newinitial]

ConvertedMcLeod=Association[Function[m,m->ConvertMcLeod[m][McLeod[m]]]/@Keys[McLeod]]

BianchiMValues={1,2,3,5,6,7,10,11,13,14,15,17,19,21,30,33,39};

IsolatedBianchiClusters=\[LeftAssociation]1->{{1},{3}},2->{{1},{3}},3->{},5->{{3},{4},{3,4}},6->{{1},{3},{4},{3,4}},7->{{3},{4},{3,4}},10->{{1},{3},{4},{7},{8},{9},{1,7},{1,8},{1,9},{3,4},{3,7},{3,8},{3,9},{4,7},{4,8},{4,9},{8,9},{1,8,9},{3,4,7},{3,4,8},{3,4,9},{3,8,9},{4,8,9},{3,4,8,9}},11->{{3},{4},{3,4}},13->{{3},{4},{9},{10},{3,4},{3,9},{3,10},{4,9},{4,10},{9,10},{3,4,9},{3,4,10},{3,9,10},{4,9,10},{3,4,9,10}},14->{{1},{3},{4},{7},{8},{9},{1,8},{1,9},{3,4},{3,7},{3,8},{3,9},{4,7},{4,8},{4,9},{7,9},{3,4,7},{3,4,8},{3,4,9},{3,7,9},{4,7,9},{3,4,7,9}},15->{{3},{4},{7},{8},{3,4},{3,7},{3,8},{4,7},{4,8},{7,8},{3,4,7},{3,4,8},{3,7,8},{4,7,8},{3,4,7,8}},17->{{3},{4},{8},{11},{12},{13},{3,4},{3,8},{3,11},{3,12},{3,13},{4,8},{4,11},{4,12},{4,13},{8,11},{8,12},{8,13},{11,12},{11,13},{12,13},{3,4,8},{3,4,11},{3,4,12},{3,4,13},{3,8,11},{3,8,12},{3,8,13},{3,11,12},{3,11,13},{3,12,13},{4,8,11},{4,8,12},{4,8,13},{4,11,12},{4,11,13},{4,12,13},{8,11,12},{8,11,13},{8,12,13},{11,12,13},{3,4,8,11},{3,4,8,12},{3,4,8,13},{3,4,11,12},{3,4,11,13},{3,4,12,13},{3,8,11,12},{3,8,11,13},{3,8,12,13},{3,11,12,13},{4,8,11,12},{4,8,11,13},{4,8,12,13},{4,11,12,13},{8,11,12,13},{3,4,8,11,12},{3,4,8,11,13},{3,4,8,12,13},{3,4,11,12,13},{3,8,11,12,13},{4,8,11,12,13},{3,4,8,11,12,13}},19->{{3},{4},{3,4}},21->{{3},{4},{9},{11},{3,4},{3,9},{3,11},{4,9},{4,11},{9,11},{3,4,9},{3,4,11},{3,9,11},{4,9,11},{3,4,9,11}},30->{{1},{3},{4},{8},{9},{10},{11},{1,9},{1,10},{1,11},{3,4},{3,8},{3,9},{3,10},{3,11},{4,8},{4,9},{4,10},{4,11},{8,11},{9,11},{10,11},{1,9,11},{1,10,11},{3,4,8},{3,4,9},{3,4,10},{3,4,11},{3,8,11},{3,9,11},{3,10,11},{4,8,11},{4,9,11},{4,10,11},{3,4,8,11},{3,4,9,11},{3,4,10,11}},33->{{3},{4},{7},{10},{11},{12},{13},{15},{3,4},{3,7},{3,10},{3,11},{3,12},{3,13},{3,15},{4,7},{4,10},{4,11},{4,12},{4,13},{4,15},{7,11},{7,12},{7,13},{7,15},{10,11},{10,12},{10,13},{10,15},{11,12},{11,13},{11,15},{12,15},{13,15},{3,4,7},{3,4,10},{3,4,11},{3,4,12},{3,4,13},{3,4,15},{3,7,11},{3,7,12},{3,7,13},{3,7,15},{3,10,11},{3,10,12},{3,10,13},{3,10,15},{3,11,12},{3,11,13},{3,11,15},{3,12,15},{3,13,15},{4,7,11},{4,7,12},{4,7,13},{4,7,15},{4,10,11},{4,10,12},{4,10,13},{4,10,15},{4,11,12},{4,11,13},{4,11,15},{4,12,15},{4,13,15},{7,11,12},{7,11,13},{7,11,15},{7,12,15},{7,13,15},{10,11,12},{10,11,13},{10,11,15},{10,12,15},{10,13,15},{11,12,15},{11,13,15},{3,4,7,11},{3,4,7,12},{3,4,7,13},{3,4,7,15},{3,4,10,11},{3,4,10,12},{3,4,10,13},{3,4,10,15},{3,4,11,12},{3,4,11,13},{3,4,11,15},{3,4,12,15},{3,4,13,15},{3,7,11,12},{3,7,11,13},{3,7,11,15},{3,7,12,15},{3,7,13,15},{3,10,11,12},{3,10,11,13},{3,10,11,15},{3,10,12,15},{3,10,13,15},{3,11,12,15},{3,11,13,15},{4,7,11,12},{4,7,11,13},{4,7,11,15},{4,7,12,15},{4,7,13,15},{4,10,11,12},{4,10,11,13},{4,10,11,15},{4,10,12,15},{4,10,13,15},{4,11,12,15},{4,11,13,15},{7,11,12,15},{7,11,13,15},{10,11,12,15},{10,11,13,15},{3,4,7,11,12},{3,4,7,11,13},{3,4,7,11,15},{3,4,7,12,15},{3,4,7,13,15},{3,4,10,11,12},{3,4,10,11,13},{3,4,10,11,15},{3,4,10,12,15},{3,4,10,13,15},{3,4,11,12,15},{3,4,11,13,15},{3,7,11,12,15},{3,7,11,13,15},{3,10,11,12,15},{3,10,11,13,15},{4,7,11,12,15},{4,7,11,13,15},{4,10,11,12,15},{4,10,11,13,15},{3,4,7,11,12,15},{3,4,7,11,13,15},{3,4,10,11,12,15},{3,4,10,11,13,15}},39->{{3},{4},{9},{10},{3,4},{3,9},{3,10},{4,9},{4,10},{9,10},{3,4,9},{3,4,10},{3,9,10},{4,9,10},{3,4,9,10}}\[RightAssociation];

Cocluster[m_,iso_]:=ConvertedMcLeod[m][[Complement[Range@Length[ConvertedMcLeod[m]],iso]]]

Cluster[m_,iso_]:=ConvertedMcLeod[m][[iso]]

GetBianchiStrips[m_,steps_:10,sizemax_:100]:=Function[iso,Strip[Cluster[m,iso],Cocluster[m,iso],10,100]]/@IsolatedBianchiClusters[m]

MatrixImage[m_]:=Rasterize@Style[MatrixForm[m,TableAlignments->{Center,Center}]//TraditionalForm,30];

Double[coords_,n_,d_:2]:=Module[{G=Coxeter[coords,d],C=coords,v=coords[[n]],C1,C2},
(* n is the index of the face about which the reflection occurs, G is original coxeter, v is the coords of that face *)
C=Delete[C,n];
C1=C;
C2=Simplify[C.R[v,d]];
C=Union[C1,C2];
Column[
If[d==2,{coords,G,Graphics[DrawCircle/@coords],C1,C2,C,Coxeter[C,d],Graphics[DrawCircle/@C]},{coords,G,C1,C2,C,Coxeter[C,d]}
]
]
]