#### Question

Consider two circles of radii $R$ and $r$ whose centres are separated by a distance $d$.

The area of overlap between them (assuming $|R-r| \le d \le R+r$) is given by the following formula, derived here.

$$
A(d; R, r) = \beta R^2 + \alpha r^2 - \frac{1}{2}r^2\sin 2\alpha - \frac{1}{2}R^2\sin 2\beta,
$$
where
$$
\cos \alpha = \frac{r^2 + d^2 - R^2}{2rd} \quad \mathrm{and} \quad \cos \beta = \frac{R^2 + d^2 - r^2}{2Rd}.
$$

Inverting this formula to find $d$ for a given overlap area $A=A_0$ cannot be achieved analytically, but a numerical solution can be found as the root of the function $A(d; R, r) - A_0$.

Write a Python function which takes arguments `A`

, the target overlap area, and `R`

and `r`

the two circle radii and returns `d`

, the distance between the circle centres giving overlap area `A`

. Use, for example, `brentq`

.

#### Solution

Click here for a solution

Here's one solution. Note that if you are using Python 2 you should ensure that `d`

, `R`

and `r`

are passed as `float`

s, perhaps by adding the line

d, R, r = float(d), float(R), float(r)

at the beginning of the `intersection_area`

function definition.

import numpy as np
from scipy.optimize import brentq
def intersection_area(d, R, r):
"""Return the area of intersection of two circles.
The circles have radii R and r, and their centres are separated by d.
"""
if d <= abs(R-r):
# One circle is entirely enclosed in the other.
return np.pi * min(R, r)**2
if d >= r + R:
# The circles don't overlap at all.
return 0
r2, R2, d2 = r**2, R**2, d**2
alpha = np.arccos((d2 + r2 - R2) / (2*d*r))
beta = np.arccos((d2 + R2 - r2) / (2*d*R))
return ( r2 * alpha + R2 * beta -
0.5 * (r2 * np.sin(2*alpha) + R2 * np.sin(2*beta))
)
def find_d(A, R, r):
"""
Find the distance between the centres of two circles giving overlap area A.
"""
# A cannot be larger than the area of the smallest circle!
if A > np.pi * min(r, R)**2:
raise ValueError("Intersection area can't be larger than the area"
" of the smallest circle")
if A == 0:
# If the circles don't overlap, place them next to each other
return R+r
if A < 0:
raise ValueError('Negative intersection area')
def f(d, A, R, r):
return intersection_area(d, R, r) - A
a, b = abs(R-r), R+r
d = brentq(f, a, b, args=(A, R, r))
return d

For example,

In [x]: find_d(0.35, 1, 0.4))
0.8489711717559927