Geometry problem involving infinite number of circles What is the sum of the areas of the grey circles? I have not made any progress so far.

 A: Consider the picture below.

On the left hand side we have the original picture which we have put in the complex plane, and on the right hand side is its image under the Möbius transformation $f(z) = \frac{z}{2-z}$. The inverse transformation is $g(z) = \frac{2z}{z+1}$. In the new coordinates the green circles have radius $\frac{1}{2}$ and center at $C_k = \frac{1}{2} + i(k + \frac{1}{2})$, $k=0,1,2,\dots$.
We have $\text{Jac}\, g(z) = |g'(z)|^2 = \frac{4}{|z+1|^4}$. Therefore the area of the disks in the original picture is
$$\sum_{k=0}^\infty \int_{B_k} \frac{4}{|(x+iy)+1|^4} \, dx \, dy,$$
where $B_k = \{z : |z - C_k| \le \frac{1}{2}\}$ are the discs in the second picture. Calculating the integral inside the sum turned out to be tedious, however, so I opted for another route. EDIT: achille hui showed how to calculate this integral in an answer to my question at Integral related to a geometry problem. This yields a shorter way to get the answer.
The points where a green circle touches the red line or the blue line are at $i(k + \frac{1}{2})$ and $1 + i(k + \frac{1}{2})$ respectively. Therefore in the original picture they are at
$$A = g(i(k+\frac{1}{2})) = \frac{2(k+\frac{1}{2})^2}{1 + (k + \frac{1}{2})^2} + i \frac{2(k+\frac{1}{2})}{1 + (k + \frac{1}{2})^2}$$
and
$$B = g(1 + i(k + \frac{1}{2})) = \frac{4 + 2(k+\frac{1}{2})^2}{4 + (k + \frac{1}{2})^2} + i \frac{2(k+\frac{1}{2})}{4 + (k+\frac{1}{2})^2}.$$
Now to find the center of the green circle we calculate the intersection of the lines $1 + t(A-1)$ and $\frac{3}{2} + s(B - \frac{3}{2})$. Real and imaginary part give us two linear equations for $s$ and $t$, and we end up with the solution
$$s = \frac{4k^2 + 4k + 17}{4k^2 + 4k + 9}, \quad t = \frac{4k^2 + 4k + 5}{4k^2 + 4k + 9}.$$
Thus the center is at
$$1 + t(A-1) = \frac{8k^2 + 8k + 6}{4k^2 + 4k + 9} + i \frac{8k + 4}{4k^2 + 4k + 9}.$$
The radius is then
$$|A - 1 - t(A-1)| = |A-1| |t-1| = \frac{4}{4k^2 + 4k + 9}.$$
Hence the answer to the problem is
$$\sum_{k=0}^\infty \pi \frac{16}{(4k^2 + 4k + 9)^2},$$
for which Wolfram Alpha gives the closed form
$$\frac{1}{16} \pi^2 (\sqrt{2} \tanh(\sqrt{2}\pi) - 2 \pi \text{sech}^2(\sqrt{2}\pi)) \approx 0.8699725\dots$$
A: This is an alternate approach to derive the area of the circles
using Descartes four circle theorem.
WOLOG, assume $R = 1$. Let us call


*

*the outer circle (with radius $r_a = 1$) as $C_a$.

*the inner green circle (with radius $r_b = \frac12$) as $C_b$.

*the $1^{st}$ gray circle (the largest one at bottom) as $C_0$.

*the $2^{nd}$ gray circle (the one above $C_b$ and $C_0$) as $C_1$.

*in general, for any $k > 0$, we will call the gray circle above $C_b$ and $C_{k-1}$ as $C_k$.

*We can reflect the figure in question vertically. For $k < 0$, we will let $C_k$ be the mirror image of the circle  $C_{-(1+k)}$.


The key is for any $k\in\mathbb{Z}$, the four circles $C_a, C_b, C_k$ and $C_{k\pm 1}$ are kissing each other. Let $r_k$ be the radius of $C_k$. We can apply Descartes four circle theorem to obtain:
$$\left(\frac{1}{r_{k\pm 1}} + \frac{1}{r_k} + \frac{1}{r_b} - \frac{1}{r_a} \right)^2
= 2 \left(\frac{1}{r_{k\pm 1}^2} + \frac{1}{r_k^2} + \frac{1}{r_b^2} + \frac{1}{r_a^2}\right)\tag{*1}$$
Since $r_a = 1$ and $r_b = \frac12$, this implies $\displaystyle\;\frac{1}{r_{k \pm 1}}\;$ are the two roots of quadratic equation
$$\left(\rho + \frac{1}{r_k} + 2 - 1 \right)^2
= 2 \left(\rho^2 + \frac{1}{r_k^2} + 2^2 + 1^2 \right)
$$
and hence
$$\frac{1}{r_{k+1}} + \frac{1}{r_{k-1}} = 2\left(\frac{1}{r_k} + 1\right)
\quad\iff\quad
\frac{1}{r_{k+1}} + \frac{1}{r_{k-1}} - \frac{2}{r_k} = 2$$
The RHS is an inhomogeneous linear recurrence relation on $\displaystyle\;\frac{1}{r_k}\;$ 
with constant term. Since the characteristic polynomial $(\lambda-1)^2$ has a double root at $1$, we know its solution must have the form $\displaystyle\;\frac{1}{r_k} = k^2 + \lambda k + \mu\;$.  By symmetry, 
$$r_k = r_{-(1+k)}\quad\implies\quad \lambda = 1 \quad\implies\quad
\frac{1}{r_k} = k(k+1) + \mu$$
To fix $\mu$, apply $(*1)$ to the case $C_a, C_b, C_0$ and $C_{-1}$, we get
$$\left(2\mu + 1\right)^2 = 2\left(2\mu^2 + 5\right)
\quad\implies\quad \mu = \frac{9}{4}
$$
From this, we get 
$$r_k = \frac{4}{4k^2 + 4k + 9}
\quad\implies\quad
\text{Area} = \sum_{k=0}^\infty \frac{16\pi}{(4k^2+4k+9)^2}$$
reproducing what's in J.J's answer. 
Update
To evaluate the sum, we will use following infinite product expansion of $\cosh(\pi x)$,
$$\cosh\pi x = \prod_{k=0}^\infty \left( 1 + \frac{x^2}{(k+\frac12)^2}\right)
$$
Taking logarithm, differentiate w.r.t $x$ and divide by $2x$, we get
$$\sum_{k=0}^\infty \frac{1}{(k+\frac12)^2 + x^2} = \frac{\pi}{2x} \tanh\pi x$$
Differentiate w.r.t $x$ and divide by $-2x$ once more, we get
$$\sum_{k=0}^\infty \frac{1}{((k+\frac12)^2 + x^2)^2} 
= \frac{\pi}{4x^3}\tanh(\pi x) - \frac{\pi^2}{4x^2\cosh^2(\pi x)}
$$
With this, we find
$$\text{Area} 
=  \sum_{k=0}^\infty \frac{\pi}{((k+\frac12)^2 + 2)^2}
=  \frac{\pi^2}{16}\left[\sqrt{2}\tanh(\pi\sqrt{2}) - \frac{2\pi}{\cosh^2(\pi\sqrt{2})}\right]
$$
Once again, this agrees with J.J's result.
