Problem Description

Write a program to investigate the math behind to circles which intersect, specifically the points of intersection .

Background & Techniques

This problem came to mind while was writing demo code to draw the Olympic rings and needed to redraw small arcs around the points where the rings intersect to give illusion of one ring passing over or under another.  It turns out that redrawing by quadrant was good enough for that application, but the problem of determining exact points of intersection   sounded,  and has been,  interesting.   A good review of your high school trigonometry if nothing else.  

The code for this program is based on a description by Paul Bourke and C code by Tim Voght  found at  http://astronomy.swin.edu.au/~pbourke/geometry/2circle/.     All I added was the conversion of the code to Delphi, a test program and  an expanded explanation of the math involved.    Mr. Bourke made several leaps that were not obvious to me, but after several hours of enjoyable study, I think I understand it.

I modified Paul’s diagram with a few extra labels like this:

P₁ and P₂ are circle centers located at (x₁,y₁) and (x₂,y₂).  dx and dy are the x and y offsets of P₂ from P₁.

dx = x₂-x₁
dy = y₂-y₁

Pythagoras tells us that the distance between P₁ and P₂ is d = sqrt(dx*dx+dy*dy).   If  d is less than then sum of the radii of the circles, r₁+r₂, then the circles may intersect.  (There is still the chance that P₂ is so close to P₁ that it is entirely contained within and still does not intersect.  This is true if and only  d<|r₁-r₂| .  Imagine that the circles a concentric, share the same center.  Then the distance from the rim of the inner circle to the outer is the absolute value of the difference in radii.  If d becomes greater than this distance the circles again intersect until it exceeds the sum of the radii.)

Assuming the circles intersect, we want to find the coordinates of the points or points of intersection, IP₁ and IP₂.  ( For clarity, I left the label off of IP₂ on the diagram above).

So we have enough information given to solve for a and h.    Considering the right triangle IP₁P₁P_3,  we can write cos(a)= a/r₁ where a is the angle at P_{1 .}

Also by considering triangle P₁P₂IP₁ , we can apply the "Cosine law" to get another expression for cos(a).

  Law of Cosines:  Cos(a)=(r₁²+d²-r₂²)/(2r₁d)

Therefore   cos(a)= a/r₁= (r₁²+d²-r₂²)/(2r₁d) and multiplying both sides by r₁ leaves

a= (r₁² - r₂² + d² ) / (2 d)

If the circles touch at one point, you can substitute d=r₁+r₂ in the above equation and prove that a=r₁ at the single point of intersection.

Once we know a, we can find the value for h as h=sqrt(r₁² - a²).

Whew, just one more critical step to find the actual coordinates of the intersection points    The trick here is to conceptually rotate point P₃ down so that it lies on the x-axis.  We don't know what that angle is yet, but we'll call it,  -beta, -b.     Then the coordinates of P₃ relative to P₁ would be (a, h).    Now if we rotate through the point (a, h) back up by angle b, it will be at the upper  intersection point.   The other point if it exists, will be the result of rotating the point (a, -h).   

So if we look up 2D rotation equations (at this site , for example, we'll find the equations for rotating at point through a given angle. 

x' = x cos b - y sin b
y' = y cos b + x sin b

 To apply the rotation up by b, we need its value, or at least the values for sin( b)  and cos(b).    These are defined by the triangle formed by sides d, dx and dy in the above diagram  so Sin(b) = dy/d and Cos(b) = dx/d.    One final consideration, since the above equations assume that we are rotating about the origin( (0,0) and we are actually rotating about P₁  at (x₁,y₁), we need to translate the final coordinates back by adding (x₁,y₁) to the rotated values. .

Our intersection points are therefore:   IP₁.x = x₁ + a*dx/d - h*dy/d  and  IP₁.y = y₁ + h*dx/d + a*dy/d  

or

IP₁.x = x₁ + (a*dx - h*dy)/d  and , IP₁.y = y₁ + (h*dx+a*dy)/d   

and for the lower intersection (a,-h) rotated by b,

IP₂.x  = x₁ +  (a*dx + h*dy)/d  and  IP₂.y =  y₁ + (-h*dx+a*dy)/d

As they say in Latin.  "Q.E.D.", which loosely translates as "We have arrived at our destination" .

I'll skip coding notes on this one - not much special and the math description wore me out!

Running/Exploring the Program 

Suggestions for Further Explorations

???