The aim of this post is to calculate the coordinates of touching points of a line tangent on two circles. As illustrated on the figure bellow, four configurations may exist in the general case.
We first consider the fact that each point lies on a circle: $$ \left \{ \begin{array}{r c l} (x_1 - x_A)^2 + (y_1 - y_A)^2 = {R_A}^2 \\ (x_2 - x_B)^2 + (y_2 - y_B)^2 = {R_B}^2 \end{array} \right . $$
We now consider the fact that each Tangent is perpendicular to the radius: $$ \left \{ \begin{array}{r c l} (x_1 - x_2)^2 + (y_1 - y_2)^2 + {R_A}^2 = (x_2 - x_A)^2 + (y_2 - y_A)^2 \\ (x_1 - x_2)^2 + (y_1 - y_2)^2 + {R_B}^2 = (x_1 - x_B)^2 + (y_1 - y_B)^2 \end{array} \right . $$
It become possible to solve this four equations and thus to find the four unknow. Unfortunatly, it is not linear and solving such system is complex.
Fortunatly, Thales configurations make it simpler and alows us to easily calculate \( \left | P_1 P_2 \right |\).
In this configuration, the point \( U \) is added on the radius of \(C_B\) in such a way that \( \left | UP_B \right | = R_B - R_A \).
The triangle \( \widehat {P_A U P_B} \) is rectangle and \( \left | P_AU \right |^2 + \left | P_BU \right |^2 = \left | P_AP_B \right |^2\). Based on this triangle, we can deduce that:
$$ \left | P_1P_2 \right | = \sqrt { \left | P_AP_B \right |^2 - (R_A - R_B)^2 }$$
From the previous equation, we can deduce that the first configuration tangent may exist only if \( \left | P_AP_B \right |^2 \ge (R_A - R_B)^2 \) ie. none of the circle is fully included in the other.
In this new configuration, points \(A'\) and \(P_1'\) are added in such a way that \( \left | P_AP_B \right | = \left | MA' \right | \) and \( \left | P_1P_2 \right | = \left | MP_1' \right | \). Note that according to Thales \( \left | P_1'A' \right | = R_A + R_B \).
As previously, the triangle \( \widehat {M P_1' A'} \) is rectangle. $$ \left | MP_1' \right |^2 + \left | P_1'A' \right |^2 = \left | MA' \right |^2 $$ $$ \left | P_1P_2 \right |^2 + (R_A + R_B)^2 = \left | P_AP_B \right |^2 $$
It becomes thus easy to deduce that: $$ \left | P_1P_2 \right | = \sqrt { \left | P_AP_B \right |^2 - (R_A + R_B)^2 }$$ From the previous equation, we can deduce that the second configuration tangent may exist only if \( \left | P_AP_B \right |^2 \ge (R_A + R_B)^2 \) ie. none of the circle is fully included in the other.
Once \( \left | P_1P_2 \right |\) is known, equations can be reformulated. Let \(L\) be equal to \( \left | P_1P_2 \right |\), equations can be rewriten:
$$ \left \{ \begin{array}{ll} (x_1 - x_A)^2 + (y_1 - y_A)^2 = {R_A}^2 \\ L^2 + {R_B}^2 = (x_1 - x_B)^2 + (y_1 - y_B)^2 \end{array} \right . $$
$$ \left \{ \begin{array}{ll} (x_2 - x_B)^2 + (y_2 - y_B)^2 = {R_B}^2 \\ L^2 + {R_A}^2 = (x_2 - x_A)^2 + (y_2 - y_A)^2 \end{array} \right . $$
We now have two independant non-linear systems to solve. Futhermore, due to the symetry, solving one will solve the whole problem.
As the systems are equivalent, we will focus on solving the first one: $$ \left \{ \begin{array}{ll} (x_1 - x_A)^2 + (y_1 - y_A)^2 = {R_A}^2 \\ (x_1 - x_B)^2 + (y_1 - y_B)^2 = L^2 + {R_B}^2 \end{array} \right .$$
It is clear that the geometrical solution is the intersection of two circles of centers \(P_A\) and \(P_B\) with respective radius of \(R_A\) and \(\sqrt{L^2+{R_B}^2}\).
The coordinates of \(P_1\) are given by:
$$ \left \{ \begin{array}{ll} x_1= \frac{x_A+x_B}{2} + \frac{(x_B-x_A)({R_A}^2-{R_1}^2)}{2D^2} \pm 2\frac{y_A-y_B}{D^2}\sigma_1 \\ y_1= \frac{y_A+y_B}{2} + \frac{(y_B-y_A)({R_A}^2-{R_1}^2)}{2D^2} \pm 2\frac{x_A-x_B}{D^2}\sigma_1 \end{array} \right .$$
with
$$ \begin{array}{ll} D = \sqrt{ (x_B-x_A)^2 + (y_B-y_A)^2 } \\ L = \sqrt { D^2 - (R_A \pm R_B)^2 } \\ R_1= \sqrt{L^2+{R_B}^2} \\ \sigma_1=\frac{1}{4}\sqrt{ (D+R_A+R_1)(D+R_A-R_1)(D-R_A+R_1)(-D+R_A+R_1) } \end{array} .$$
In the same way, it is possible to deduce the coordinates of \(P_2\):
$$ \left \{ \begin{array}{ll} x_2= \frac{x_B+x_A}{2} + \frac{(x_A-x_B)({R_B}^2-{R_2}^2)}{2D^2} \pm 2\frac{y_B-y_A}{D^2}\sigma_2 \\ y_2= \frac{y_B+y_A}{2} + \frac{(y_A-y_B)({R_B}^2-{R_2}^2)}{2D^2} \pm 2\frac{x_B-x_A}{D^2}\sigma_2 \end{array} \right .$$
with
$$ \begin{array}{ll} D = \sqrt{ (x_B-x_A)^2 + (y_B-y_A)^2 } \\ L = \sqrt { D^2 - (R_B \pm R_A)^2 } \\ R_2= \sqrt{L^2+{R_A}^2} \\ \sigma_2=\frac{1}{4}\sqrt{ (D+R_B+R_2)(D+R_B-R_2)(D-R_B+R_2)(-D+R_B+R_2) } \end{array} .$$
A quick Matlab test:
Matlab source code (example on this page) can be download here: