(*
	(c) Alex Kontorovich 2012
*)

(* This code is for drawing various pictures in Apollonian-3 and Apollonian-5 gaskets *)

(* First draw an Apollonian-3 gasket: *)

(* Parameters *)
sizemax=300;  (* Largest curvature to draw *)
levelmax=40;  (* Largest level to go down *)
labels=True; (* Whether to label the circle by its curvature *)
fscale=700; (* Size of the text inside a circle *)
debug=False;

b0={-1,4,4,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)
b0={-1,2,2,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)
b0={-2,4,5,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)
b0={-4,7,10,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)



 (* Check that this is an integer *)
Sqrt[(b0[[1]]b0[[2]]+b0[[1]]b0[[3]]+b0[[2]]b0[[3]])/2];


(* Code *)

R2ToC[v_]:=v[[1]]+I v[[2]];
CToR2[z_]:={Re[z],Im[z]};

(* Initializations *)
(* Find circle center and rad so that reflection of "standard" pic is the given one; use last solution found *)
{X0,Y0,K}={x0,y0,k}/.Flatten[Take[Solve[
-1/b0[[1]]==k^2/(2(y0-1))
&&
1/b0[[2]]==k^2/(2(y0+1))
&&
1/b0[[3]]==k^2/(x0^2+y0^2-1)
&&k>0&&x0>0
,
{x0,y0,k}],1]];

(* To flip a point *)
circInvP[cir_,p_]:=Simplify[{cir[[1]],cir[[2]]}+cir[[3]]^2(p-{cir[[1]],cir[[2]]})/Norm[(p-{cir[[1]],cir[[2]]})]^2];
(* To flip a circle, given {x,y,rad} *)
circInvCir[cirFix_,cirToMove_]:=Module[{near,far,nearTil,farTil,unitV,ans},
unitV=(Drop[cirToMove,-1]-Drop[cirFix,-1])/Norm[(Drop[cirToMove,-1]-Drop[cirFix,-1])];

near=Drop[cirToMove,-1]-cirToMove[[3]] unitV;
far=Drop[cirToMove,-1]+cirToMove[[3]] unitV;

nearTil=circInvP[cirFix,near];
farTil=circInvP[cirFix,far];

ans=Flatten[{(nearTil+farTil)/2,Norm[nearTil-farTil]/2}]//FullSimplify;

ans
];

(* Get "standard" sextuplet configuration *)
xx0=xx/.Flatten[Take[Solve[ xx^2+(1/2)^2==(3/2)^2&&xx>0,xx],1]];

stdPic=Show[Graphics[
{Dotted, Red,
Line[{{-2,-1},{5,-1}}]
,
Line[{{-2,1},{5,1}}]
,
Circle[{0,0},1  ]
,
Circle[{xx0,1/2},1/2  ]
,
Circle[{xx0,-1/2},1/2  ]
,
Circle[{2xx0,0},1  ]
,
Point[{0,0}]
}
]
];

(* Fill out b0 = curvatures  *)
b0[[4]]=1/circInvCir[{X0,Y0,K},{xx0,-1/2,1/2}][[3]];
b0[[5]]=1/circInvCir[{X0,Y0,K},{xx0,1/2,1/2}][[3]];
b0[[6]]=1/circInvCir[{X0,Y0,K},{2 xx0,0,1}][[3]];

(* Fill out w0 = curv centers *)
w0={
(* flip line y=1 *)
b0[[1]]
R2ToC[
{X0,Y0}
-
1/2 ({X0,Y0}-{X0,1}) K^2/ Norm[{X0,Y0}-{X0,1}]^2
]

,
(* flip line y=-1 *)
b0[[2]]
R2ToC[
{X0,Y0}
-
1/2 ({X0,Y0}-{X0,-1}) K^2/ Norm[{X0,Y0}-{X0,-1}]^2
]
,
b0[[3]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{0,0,1}],-1]
]
,
b0[[4]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{xx0,-1/2,1/2}],-1]
]
,
b0[[5]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{xx0,1/2,1/2}],-1]
]
,
b0[[6]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{2xx0,0,1}],-1]
]
}//FullSimplify;

pic0=(*Show[Graphics[*)
{
Circle[CToR2[w0[[1]]/b0[[1]]  ],1/Abs[b0[[1]]  ]  ]
,
Circle[CToR2[w0[[2]]/b0[[2]]  ],1/b0[[2]]  ]
,
Circle[CToR2[w0[[3]]/b0[[3]]  ],1/b0[[3]]  ]
,
Circle[CToR2[w0[[4]]/b0[[4]]  ],1/b0[[4]]  ]
,
Circle[CToR2[w0[[5]]/b0[[5]]  ],1/b0[[5]]  ]
,
Circle[CToR2[w0[[6]]/b0[[6]]  ],1/b0[[6]]  ]
(*
,
Text["2",CToR2[w0[[2]]/b0[[2]]  ]]
,
Text["3",CToR2[w0[[3]]/b0[[3]]  ]]
,
Text["1'",CToR2[w0[[4]]/b0[[4]]  ]]
,
Text["2'",CToR2[w0[[5]]/b0[[5]]  ]]
,
Text["3'",CToR2[w0[[6]]/b0[[6]]  ]]
*)
};
(* JUST FOR TESTING *
]
,
Graphics[
{Dotted,Blue,Circle[{X0,Y0},K  ],
Point[{X0,Y0}]
}
]
,
stdPic
]
* *)

For[j=2,j≤6,j++,
If[labels,AppendTo[pic0,
Text[
Style[Abs[b0[[j]]],Floor[fscale/Abs[b0[[j]]]]],{ Re[((w0[[j]])/(b0[[j]]))], Im[((w0[[j]])/(b0[[j]]))]}
]
];
]
];
 (* Show[Graphics[pic0]] *)

(* All generating matrices *)
(* 123 *)
m123={
{1,0,0,0,0,0},
{0,1,0,0,0,0},
{0,0,1,0,0,0},
{2,4,4,-1,0,0},
{4,2,4,0,-1,0},
{4,4,2,0,0,-1}
};
(* 123' *)
m123p={
{1,0,0,0,0,0},
{0,1,0,0,0,0},
{4,4,-1,0,0,2},
{2,4,0,-1,0,4},
{4,2,0,0,-1,4},
{0,0,0,0,0,1}
};
(* 12'3 *)
m12p3={
{1,0,0,0,0,0},
{4,-1,4,0,2,0},
{0,0,1,0,0,0},
{2,0,4,-1,4,0},
{0,0,0,0,1,0},
{4,0,2,0,4,-1}
};
(* 12'3' *)
m12p3p={
{1,0,0,0,0,0},
{4,-1,0,0,2,4},
{4,0,-1,0,4,2},
{2,0,0,-1,4,4},
{0,0,0,0,1,0},
{0,0,0,0,0,1}
};
(* 1'23 *)
m1p23={
{-1,4,4,2,0,0},
{0,1,0,0,0,0},
{0,0,1,0,0,0},
{0,0,0,1,0,0},
{0,2,4,4,-1,0},
{0,4,2,4,0,-1}
};
(* 1'23' *)
m1p23p={
{-1,4,0,2,0,4},
{0,1,0,0,0,0},
{0,4,-1,4,0,2},
{0,0,0,1,0,0},
{0,2,0,4,-1,4},
{0,0,0,0,0,1}
};
(* 1'2'3 *)
m1p2p3={
{-1,0,4,2,4,0},
{0,-1,4,4,2,0},
{0,0,1,0,0,0},
{0,0,0,1,0,0},
{0,0,0,0,1,0},
{0,0,2,4,4,-1}
};
(* 1'2'3p *)
m1p2p3p={
{-1,0,0,2,4,4},
{0,-1,0,4,2,4},
{0,0,-1,4,4,2},
{0,0,0,1,0,0},
{0,0,0,0,1,0},
{0,0,0,0,0,1}
};



(* Testing 
b=m.b0;
w=m.w0;
pic=pic0;
For[j=1,j≤6,j++,
pic=Union[pic,
{Circle[CToR2[w[[j]]/b[[j]]  ],1/b[[j]]  ]}
];
If[labels,
pic=Union[pic,{
Text[
Style[Abs[b[[j]]],Floor[fscale/Abs[b[[j]]]]],{ Re[((w[[j]])/(b[[j]]))], Im[((w[[j]])/(b[[j]]))]}
]}
];
];
];
Show[Graphics[pic]]
*)

(* Collect all matrices *)
Smat={
m123,
m123p,
m12p3,
m12p3p,
m1p23,
m1p23p,
m1p2p3,
m1p2p3p
};
MatrixForm/@Smat;

(* Which circles are new corr to matrix? *)
new={
{4,5,6},
{4,5,3},
{4,2,6},
{4,2,3},
{1,5,6},
{1,5,3},
{1,2,6},
{1,2,3}
};


(* The main routine *)
morecircs[from_,level_]:=If[level<levelmax,
If[debug,Print["in level ",level];];
Module[{j,k},For[j=1,j≤8,j++,If[
from≠j,

If[debug,Print["flipping by matrix ",j];];

b=Smat[[j]].b;
w=Smat[[j]].w;
If[
(* smallest new circle isn't too big *)
Sort[b][[4]]<sizemax,

curvs=Union[curvs,b];

(* Loop through new circles *)
For[k=1,k≤3,k++,
AppendTo[pic,
Circle[CToR2[w[[ new[[j,k]] ]]/b[[ new[[j,k] ] ]]  ],1/b[[ new[[j,k]] ]]  ]
];
If[labels,
AppendTo[pic,
Text[
Style[Abs[b[[ new[[j,k]]  ]]],Floor[fscale/Abs[b[[ new[[j,k]] ]]]]],{ Re[((w[[ new[[j,k]] ]])/(b[[new[[j,k]] ]]))], Im[((w[[ new[[j,k]] ]])/(b[[ new[[j,k]] ]]))]}
]
];
];
];


morecircs[j,level+1];
];
b=Smat[[j]].b;
w=Smat[[j]].w;
]
]
];
];


(* Run the program *)
pic=pic0;
b=b0;
curvs=b0;
w=w0;
morecircs[0,1];
Length[pic]
Show[Graphics[pic]]

(* Then draw a dual configuration: *)

(* Parameters *)
sizemax=50;  (* Largest curvature to draw *)
levelmax=40;  (* Largest level to go down *)
labels=False; (* Whether to label the circle by its curvature *)
fscale=700; (* Size of the text inside a circle *)
fsize=20;
debug=False;


CirThru[a_,b_,c_]:=Solve[
((a[[1]]-x0)^2+(a[[2]]-y0)^2==(r)^2)
&&
((b[[1]]-x0)^2+(b[[2]]-y0)^2==(r)^2)
&&
((c[[1]]-x0)^2+(c[[2]]-y0)^2==(r)^2)

,{x0,y0,r}];

b0={-1,4,4,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)
b0={-1,2,2,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)
b0={-2,4,5,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)
b0={-4,7,10,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)

 (* Check that this is an integer *)
Sqrt[(b0[[1]]b0[[2]]+b0[[1]]b0[[3]]+b0[[2]]b0[[3]])/2];


(* Code *)

R2ToC[v_]:=v[[1]]+I v[[2]];
CToR2[z_]:={Re[z],Im[z]};

(* Initializations *)
(* Find circle center and rad so that reflection of "standard" pic is the given one; use last solution found *)
{X0,Y0,K}={x0,y0,k}/.Flatten[Take[Solve[
-1/b0[[1]]==k^2/(2(y0-1))
&&
1/b0[[2]]==k^2/(2(y0+1))
&&
1/b0[[3]]==k^2/(x0^2+y0^2-1)
&&k>0&&x0>0
,
{x0,y0,k}],1]];

(* To flip a point *)
circInvP[cir_,p_]:=Simplify[{cir[[1]],cir[[2]]}+cir[[3]]^2(p-{cir[[1]],cir[[2]]})/Norm[(p-{cir[[1]],cir[[2]]})]^2];
(* To flip a circle, given {x,y,rad} *)
circInvCir[cirFix_,cirToMove_]:=Module[{near,far,nearTil,farTil,unitV,ans},
unitV=(Drop[cirToMove,-1]-Drop[cirFix,-1])/Norm[(Drop[cirToMove,-1]-Drop[cirFix,-1])];

near=Drop[cirToMove,-1]-cirToMove[[3]] unitV;
far=Drop[cirToMove,-1]+cirToMove[[3]] unitV;

nearTil=circInvP[cirFix,near];
farTil=circInvP[cirFix,far];

ans=Flatten[{(nearTil+farTil)/2,Norm[nearTil-farTil]/2}]//FullSimplify;

ans
];

(* Get "standard" sextuplet configuration *)
xx0=xx/.Flatten[Take[Solve[ xx^2+(1/2)^2==(3/2)^2&&xx>0,xx],1]];

stdPic=Show[Graphics[
{Dotted, Red,
Line[{{-2,-1},{5,-1}}]
,
Line[{{-2,1},{5,1}}]
,
Circle[{0,0},1  ]
,
Circle[{xx0,1/2},1/2  ]
,
Circle[{xx0,-1/2},1/2  ]
,
Circle[{2xx0,0},1  ]
,
Point[{0,0}]
}
]
];

(* Fill out b0 = curvatures  *)
b0[[4]]=1/circInvCir[{X0,Y0,K},{xx0,-1/2,1/2}][[3]];
b0[[5]]=1/circInvCir[{X0,Y0,K},{xx0,1/2,1/2}][[3]];
b0[[6]]=1/circInvCir[{X0,Y0,K},{2 xx0,0,1}][[3]];

(* Fill out w0 = curv centers *)
w0={
(* flip line y=1 *)
b0[[1]]
R2ToC[
{X0,Y0}
-
1/2 ({X0,Y0}-{X0,1}) K^2/ Norm[{X0,Y0}-{X0,1}]^2
]

,
(* flip line y=-1 *)
b0[[2]]
R2ToC[
{X0,Y0}
-
1/2 ({X0,Y0}-{X0,-1}) K^2/ Norm[{X0,Y0}-{X0,-1}]^2
]
,
b0[[3]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{0,0,1}],-1]
]
,
b0[[4]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{xx0,-1/2,1/2}],-1]
]
,
b0[[5]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{xx0,1/2,1/2}],-1]
]
,
b0[[6]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{2xx0,0,1}],-1]
]
}//FullSimplify;

pic0={
Circle[CToR2[w0[[1]]/b0[[1]]  ],1/Abs[b0[[1]]  ]  ]
,
Circle[CToR2[w0[[2]]/b0[[2]]  ],1/b0[[2]]  ]
,
Circle[CToR2[w0[[3]]/b0[[3]]  ],1/b0[[3]]  ]
,
Circle[CToR2[w0[[4]]/b0[[4]]  ],1/b0[[4]]  ]
,
Circle[CToR2[w0[[5]]/b0[[5]]  ],1/b0[[5]]  ]
,
Circle[CToR2[w0[[6]]/b0[[6]]  ],1/b0[[6]]  ]
};

pic2={
Text[Style[Subscript[κ,2],fsize],CToR2[w0[[2]]/b0[[2]]  ]]
,
Text[Style[Subscript[κ,3],fsize],CToR2[w0[[3]]/b0[[3]]  ]]
,
Text[Style[Subscript[κ,1]',fsize],CToR2[w0[[4]]/b0[[4]]  ]]
,
Text[Style[Subscript[κ,2]',fsize],CToR2[w0[[5]]/b0[[5]]  ]]
,
Text[Style[Subscript[κ,3]',fsize],CToR2[w0[[6]]/b0[[6]]  ]]

};

pic1=
Graphics[
{Dotted,Blue,Circle[{X0,Y0},K  ],
Point[{X0,Y0}]
}
];
(* *)


For[j=2,j≤6,j++,
If[labels,AppendTo[pic0,
Text[
Style[Abs[b0[[j]]],Floor[fscale/Abs[b0[[j]]]]],{ Re[((w0[[j]])/(b0[[j]]))], Im[((w0[[j]])/(b0[[j]]))]}
]
];
]
];
 (* Show[Graphics[pic0]] *)

(* All generating matrices *)
(* 123 *)
m123={
{1,0,0,0,0,0},
{0,1,0,0,0,0},
{0,0,1,0,0,0},
{2,4,4,-1,0,0},
{4,2,4,0,-1,0},
{4,4,2,0,0,-1}
};
(* 123' *)
m123p={
{1,0,0,0,0,0},
{0,1,0,0,0,0},
{4,4,-1,0,0,2},
{2,4,0,-1,0,4},
{4,2,0,0,-1,4},
{0,0,0,0,0,1}
};
(* 12'3 *)
m12p3={
{1,0,0,0,0,0},
{4,-1,4,0,2,0},
{0,0,1,0,0,0},
{2,0,4,-1,4,0},
{0,0,0,0,1,0},
{4,0,2,0,4,-1}
};
(* 12'3' *)
m12p3p={
{1,0,0,0,0,0},
{4,-1,0,0,2,4},
{4,0,-1,0,4,2},
{2,0,0,-1,4,4},
{0,0,0,0,1,0},
{0,0,0,0,0,1}
};
(* 1'23 *)
m1p23={
{-1,4,4,2,0,0},
{0,1,0,0,0,0},
{0,0,1,0,0,0},
{0,0,0,1,0,0},
{0,2,4,4,-1,0},
{0,4,2,4,0,-1}
};
(* 1'23' *)
m1p23p={
{-1,4,0,2,0,4},
{0,1,0,0,0,0},
{0,4,-1,4,0,2},
{0,0,0,1,0,0},
{0,2,0,4,-1,4},
{0,0,0,0,0,1}
};
(* 1'2'3 *)
m1p2p3={
{-1,0,4,2,4,0},
{0,-1,4,4,2,0},
{0,0,1,0,0,0},
{0,0,0,1,0,0},
{0,0,0,0,1,0},
{0,0,2,4,4,-1}
};
(* 1'2'3p *)
m1p2p3p={
{-1,0,0,2,4,4},
{0,-1,0,4,2,4},
{0,0,-1,4,4,2},
{0,0,0,1,0,0},
{0,0,0,0,1,0},
{0,0,0,0,0,1}
};


(* Collect all matrices *)
Smat={
m123,
m123p,
m12p3,
m12p3p,
m1p23,
m1p23p,
m1p2p3,
m1p2p3p
};
MatrixForm/@Smat;



(* Which circles are new corr to matrix? *)
new={
{4,5,6},
{4,5,3},
{4,2,6},
{4,2,3},
{1,5,6},
{1,5,3},
{1,2,6},
{1,2,3}
};

(*  pjkS = standard pt of intersection of bj and bk *)
p13S={0,1};
p23S={0,-1};
p12pS={xx0,1};
p1p2S={xx0,-1};
p13pS={2xx0,1};
p23pS={2xx0,-1};
p2p3S={2xx0/3,1/3};
p1p3S={2xx0/3,-1/3};
p1p2pS={xx0,0};
p2p3pS={4xx0/3,1/3};
p1p3pS={4xx0/3,-1/3};

(* get dual circs *)
Dual23={x0,y0,r}/.CirThru[
{X0,Y0},
circInvP[{X0,Y0,K},p13S],
circInvP[{X0,Y0,K},p23S]
][[2]];
Dual23p={x0,y0,r}/.CirThru[
{X0,Y0},
circInvP[{X0,Y0,K},p13pS],
circInvP[{X0,Y0,K},p23pS]
][[2]];
Dual2p3={x0,y0,r}/.CirThru[
circInvP[{X0,Y0,K},p12pS],
circInvP[{X0,Y0,K},p2p3S],
circInvP[{X0,Y0,K},p13S]
][[2]];
Dual2p3p={x0,y0,r}/.CirThru[
circInvP[{X0,Y0,K},p12pS],
circInvP[{X0,Y0,K},p2p3pS],
circInvP[{X0,Y0,K},p13pS]
][[2]];
Dualp23={x0,y0,r}/.CirThru[
circInvP[{X0,Y0,K},p1p2S],
circInvP[{X0,Y0,K},p23S],
circInvP[{X0,Y0,K},p1p3S]
][[2]];
Dualp23p={x0,y0,r}/.CirThru[
circInvP[{X0,Y0,K},p1p2S],
circInvP[{X0,Y0,K},p23pS],
circInvP[{X0,Y0,K},p1p3pS]
][[2]];
Dualp2p3={x0,y0,r}/.CirThru[
circInvP[{X0,Y0,K},p1p2pS],
circInvP[{X0,Y0,K},p2p3S],
circInvP[{X0,Y0,K},p1p3S]
][[2]];
Dualp2p3p={x0,y0,r}/.CirThru[
circInvP[{X0,Y0,K},p1p2pS],
circInvP[{X0,Y0,K},p2p3pS],
circInvP[{X0,Y0,K},p1p3pS]
][[2]];

circ[v_]:=Circle[Drop[v,-1],v[[3]]  ];


Show[Graphics[pic0]
,
(*
Graphics[
Point[
{
p13S,
p23S,
p12pS,
p1p2S,
p13pS,
p23pS,
p2p3S,
p1p3S,
p1p2pS,
p2p3pS,
p1p3pS
}
]
]
,*)
Graphics[
{Red,Thick,Dotted,
circ[Dual23],
circ[Dual23p],
circ[Dual2p3],
circ[Dual2p3p],
circ[Dualp23],
circ[Dualp23p],
circ[Dualp2p3],
circ[Dualp2p3p]
}
]
]

(* Now draw a labelling *)

fsize=20;
pic={};
m=Smat[[1]];
b=m.b0;
w=m.w0;
j=1;
(* Testing *)
For[k=1,k≤3,k++,
AppendTo[pic,

Circle[CToR2[w[[  new[[j,k]]  ]]/b[[ new[[j,k]] ]]  ],1/b[[ new[[j,k]] ]]  ]
];

];

(* *)


j=4;
AppendTo[pic0,Text[Style[Subscript[κ,1]'',fsize],{ Re[((w[[ j ]])/(b[[ j  ]]))], Im[((w[[ j  ]])/(b[[ j ]]))]}]];
j=5;
AppendTo[pic0,Text[Style[Subscript[κ,2]'',fsize],{ Re[((w[[ j=5  ]])/(b[[ j  ]]))], Im[((w[[ j  ]])/(b[[ j ]]))]}]];
j=6;
AppendTo[pic0,Text[Style[Subscript[κ,3]'',fsize],{ Re[((w[[ j=6  ]])/(b[[ j  ]]))], Im[((w[[ j  ]])/(b[[ j ]]))]}]];
Clear[j,k,b,m,w];
(*
Show[Graphics[pic]]
*)



(* Dual Circle 123 *)
DualCirc123={x0,y0,r}/.CirThru[
circInvP[{X0,Y0,K},{0,1}],
circInvP[{X0,Y0,K},{0,-1}],
{X0,Y0}
][[2]];

Show[Graphics[pic0],Graphics[ pic2], 
Graphics[Text[ Style[Subscript[κ,1],fsize],CToR2[ w0[[1]]  ]/b0[[1]]  + .8{1,-1}/b0[[1]]  ] ]
,
Graphics[{Blue,Thick,pic}]
,
Graphics[{Red,PointSize[Large]
,
Point[  circInvP[{X0,Y0,K},{0,1}] ]
,
Point[  circInvP[{X0,Y0,K},{0,-1}] ]
,
Point[ {X0,Y0} ]

}]
,
Graphics[{Red,Dotted,Thick,Circle[Drop[DualCirc123,-1],DualCirc123[[3]],{3Pi/2+.8,2Pi+.05}]}]
]

(* Finally draw an orbit in H^3 *)

sphInv[p_,sphJ_]:=Simplify[{sphJ[[1]],sphJ[[2]],0}+sphJ[[3]]^2(p-{sphJ[[1]],sphJ[[2]],0})/Norm[(p-{sphJ[[1]],sphJ[[2]],0})]^2];
p0={1,3,1};

sph=N[{
Dual23,
Dual23p,
Dual2p3,
Dual2p3p,
Dualp23,
Dualp23p,
Dualp2p3,
Dualp2p3p
}];
numSph=Length[sph];
heightMin=1/20000;

p=p0;

points={p0};

morePts[from_]:=
Module[{j},
For[j=1,j≤numSph,j++,
If[from≠j,
p=sphInv[p,sph[[j]]];
AppendTo[points,p];
If[p[[3]]>heightMin,
morePts[j];
];
p=sphInv[p,sph[[j]]];
]
]
];

(* Run the Code *)
morePts[0]

pts=Select[points,#[[3]]<1/500&];

ptsBigg=Select[points,#[[3]]>1/50&];

ptsBig=Select[points,1/50>#[[3]]>1/500&];


Show[
Graphics3D[{PointSize[0.0002],
Point[pts]}]
,
Graphics3D[{PointSize[0.002],
Point[ptsBig]}]
,
Graphics3D[{Opacity[1],PointSize[.01],
Point[ptsBigg]}]
,
Graphics3D[{PointSize[Large],
Point[p0]}]
,
Graphics3D[
{
Opacity[0.4],
Table[
Sphere[{sph[[l,1]] , sph[[l,2]],0  }, sph[[l,3]] ]
,
{l,1,numSph}]
}
]
,
Table[
ParametricPlot3D[{sph[[l,1]]+sph[[l,3]]Cos[t],sph[[l,2]]+sph[[l,3]]Sin[t],0},{t,0,2Pi},PlotStyle->{Red,Dotted}]
,
{l,1,numSph}]
,


Graphics3D[Polygon[{{-1,2,0},{2,2,0},{2,5,0},{-1,5,0}}]],

Axes -> False, Boxed -> False, 
  ImageSize -> {1147.3651724317267, 715.}, PlotRange -> {{-0.5, 1.9}, {2.5, 4.5}, {0, 1}}, 
  ViewPoint -> {0.3266363398197495, -3.1177426077002073, 1.2739661430508549}, 
  ViewVertical -> {0.057708985506855104, -0.28729236448743445, 2.326084629943947}]

(* Now for an Apollonian 5-gasket: *)

fscale=3 600; (* Size of the text inside a circle *)
labels=True; (* Whether to label the circle by its curvature *)

b0={-3,4,15,blank,blank,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)
b0={-1,2,5,blank,blank,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)
b0=3{-3,6,7,blank,blank,blank,blank,blank}; (* Three curvatures determine everything; first must be <0 *)
(* Check that this is an integer *)
Sqrt[(b0[[1]]b0[[2]]+b0[[1]]b0[[3]]+b0[[2]]b0[[3]])/3];

R2ToC[v_]:=v[[1]]+I v[[2]];
CToR2[z_]:={Re[z],Im[z]};

(* To flip a point *)
circInvP[cir_,p_]:=Simplify[{cir[[1]],cir[[2]]}+cir[[3]]^2(p-{cir[[1]],cir[[2]]})/Norm[(p-{cir[[1]],cir[[2]]})]^2];
(* To flip a circle, given {x,y,rad} *)
circInvCir[cirFix_,cirToMove_]:=Module[{near,far,nearTil,farTil,unitV,ans},
unitV=(Drop[cirToMove,-1]-Drop[cirFix,-1])/Norm[(Drop[cirToMove,-1]-Drop[cirFix,-1])];

near=Drop[cirToMove,-1]-cirToMove[[3]] unitV;
far=Drop[cirToMove,-1]+cirToMove[[3]] unitV;

nearTil=circInvP[cirFix,near];
farTil=circInvP[cirFix,far];

ans=Flatten[{(nearTil+farTil)/2,Norm[nearTil-farTil]/2}]//FullSimplify;

ans
];
(* Get "standard" sextuplet configuration *)
xx0=xx/.Flatten[Take[Solve[ xx^2+(2/3)^2==(4/3)^2&&xx>0,xx],1]];

stdPic=Show[Graphics[
{Dotted, Red,
Line[{{-2,-1},{5,-1}}]
,
Line[{{-2,1},{5,1}}]
,
Circle[{0,0},1  ]
,
Circle[{xx0,1/6},1/6  ]
,
Circle[{xx0,-1/6},1/6  ]
,
Circle[{xx0,2/3},1/3  ]
,
Circle[{xx0,-2/3},1/3  ]
,
Circle[{2xx0,0},1  ]
,
Point[{0,0}]
}
]
];

(* Initializations *)
(* Find circle center and rad so that reflection of "standard" pic is the given one; use last solution found *)
{X0,Y0,K}={x0,y0,k}/.Flatten[Take[Solve[
-1/b0[[3]]==k^2/(2(y0-1))
&&
1/b0[[2]]==k^2/(2(y0+1))
&&
1/b0[[1]]==k^2/(x0^2+y0^2-1)
&&k>0&&x0>0
,
{x0,y0,k}],1]];
(* Fill out b0 = curvatures  *)
b0[[4]]=1/circInvCir[{X0,Y0,K},{xx0,2/3,1/3}][[3]];
b0[[5]]=1/circInvCir[{X0,Y0,K},{xx0,1/6,1/6}][[3]];
b0[[6]]=1/circInvCir[{X0,Y0,K},{xx0,-1/6,1/6}][[3]];
b0[[7]]=1/circInvCir[{X0,Y0,K},{xx0,-2/3,1/3}][[3]];
b0[[8]]=1/circInvCir[{X0,Y0,K},{2 xx0,0,1}][[3]];


(* Fill out w0 = curv centers *)
w0={
b0[[1]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{0,0,1}],-1]
]
,
(* flip line y=-1 *)
b0[[2]]
R2ToC[
{X0,Y0}
-
1/2 ({X0,Y0}-{X0,-1}) K^2/ Norm[{X0,Y0}-{X0,-1}]^2
]
,(* flip line y=1 *)
b0[[3]]
R2ToC[
{X0,Y0}
-
1/2 ({X0,Y0}-{X0,1}) K^2/ Norm[{X0,Y0}-{X0,1}]^2
]
,
b0[[4]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{xx0,2/3,1/3}],-1]
]
,
b0[[5]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{xx0,1/6,1/6}],-1]
]
,
b0[[6]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{xx0,-1/6,1/6}],-1]
]
,
b0[[7]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{xx0,-2/3,1/3}],-1]
]
,
b0[[8]]
R2ToC[
Drop[circInvCir[{X0,Y0,K},{2xx0,0,1}],-1]
]
}//FullSimplify;

pic0=(*Show[Graphics[*)
{
Table[
Circle[CToR2[w0[[j]]/b0[[j]]  ],1/Abs[b0[[j]]  ]  ],{j,1,8}]
};
Show[Graphics[pic0],
Graphics[
Table[
Text[Style[Subscript[κ,j],Floor[fscale/b0[[j]]]],CToR2[w0[[j]]/b0[[j]]  ]],{j,2,8}]]
,
Graphics[
Text[Style[Subscript[κ,1],Floor[fscale/Abs[b0[[2]]]]],CToR2[w0[[1]]/b0[[1]]  ]-1/b0[[1]] (shift={-1.2,.7})]
]
]


(* Code *)
For[j=2,j≤8,j++,
If[labels,AppendTo[pic0,
Text[
Style[Abs[b0[[j]]],Floor[fscale/Abs[b0[[j]]]]],{ Re[((w0[[j]])/(b0[[j]]))], Im[((w0[[j]])/(b0[[j]]))]}
]
];
]
];
 Show[Graphics[pic0],
Graphics[
Text[
Style[b0[[1]],Floor[fscale/Abs[b0[[2]]]]],CToR2[w0[[1]]/b0[[1]]  ]-1/b0[[1]] shift
]
]
]