Circle Back

This is somewhat in the same spirit as Pack It In and Eyes Wide, in that it shows the stereographic projection of a collection of non-overlapping circles on the sphere. In this case, though, it's an infinite family of circles on the sphere (though of course I'm only showing a finite subcollection).

This family is defined recursively, starting from 16 equally-spaced points on the equator: each is the center of a circle of (spherical) radius $\pi/16$ (i.e., half the distance between adjacent points). Next, I want to put 16 equally-spaced points on a circle of longitude in the northern hemisphere so that the circles centered at those points which are just tangent will also be tangent to the corresponding circles from the previous iteration.

Now, the spherical distance between $n$ equally-spaced points at height $z$ is given by

FullSimplify[
 VectorAngle[{Sqrt[1 - z^2] 1, 0, z}, {Sqrt[1 - z^2] Cos[2 ?/n], Sqrt[1 - z^2] Sin[2 ?/n], z}], 
 z ? Reals && -1 < z < 1 && n ? Integers && n >= 2]

which comes out to be $\arccos\left(z^2-\left(z^2-1\right) \cos \left(\frac{2 \pi }{n}\right)\right)$. So we want to work out when half of that quantity (the radius of the circles on the second row) plus $\pi/16$ (the radius of the circles on the first row) is equal to the distance between points on the first row and the corresponding point on the second row.

For each successive row we do the same computation, so we can recursively compute the correct $z$-values of the rows up to level $k$ using NestList[]:

NestList[
 Re[x] /.
   FindRoot[
    ArcCos[#^2 - (-1 + #^2) Cos[(2 ?)/n]]/2 + 
      1/2 ArcCos[x^2 - (-1 + x^2) Cos[(2 ?)/n]] ==
     VectorAngle[{Sqrt[1 - x^2], 0, x}, {Sqrt[1 - #^2], 0, #}],
    {x, 1}] &,
 0, k]

So now we can put it all together (to produce a not-completely-unresponsive Manipulate, I use depth = 6 in the below code, which leaves big gaps near the poles of the sphere; for the final animation I used depth = 16):

DynamicModule[{p, b, n = 16, depth = 6, zList, verts,
  cols = RGBColor /@ {"#14FFEC", "#323232"}},
 zList = NestList[
   Re[x] /.
     FindRoot[
      ArcCos[#^2 - (-1 + #^2) Cos[(2 ?)/n]]/2 + 
        1/2 ArcCos[x^2 - (-1 + x^2) Cos[(2 ?)/n]] ==
       VectorAngle[{Sqrt[1 - x^2], 0, x}, {Sqrt[1 - #^2], 0, #}],
      {x, 1}] &,
   0, depth];
 verts = Flatten[
   Table[
    {Sqrt[1 - zList[[i]]^2] Cos[s + ?/n], Sqrt[1 - zList[[i]]^2] Sin[s + ?/n], sgn*zList[[i]]},
    {i, 1, Length[zList]}, {sgn, {-1, 1}}, {s, 2. ?/n, 2 ?, 2 ?/n}],
   2];
 Manipulate[
  p = verts.RotationMatrix[-?, {1, -1, 0}];
  b = NullSpace[{#}] & /@ p;
  Graphics[
   {EdgeForm[None],
    Table[
     {cols[[1]],
      Polygon[
       Table[
        With[{t = ArcCos[verts[[i, 3]]^2 - (-1 + verts[[i, 3]]^2) Cos[(2 ?)/n]]},
         Stereo[
          Cos[t/2] p[[i]] + Sin[t/2] (Cos[s] b[[i, 1]] + Sin[s] b[[i, 2]])]],
        {s, 0., 2 ?, 2 ?/200}]]},
     {i, 1, Length[p]}]},
   PlotRange -> 3.5, Background -> cols[[-1]], ImageSize -> 540],
  {?, 0, ?}]
 ]