Coordinates of intersection of circle and sphere

How do I calculate the x, y, & z coordinates of the two points where a sphere intersects a circle? The coordinates of the center of both are known and the radii are known.

The center of the circle and the center of the sphere are not the same point.

The circle and the sphere are not tangent to each other.

The center of the sphere is not on the perpendicular through the center of the circle. I.e., the solution is not the circle, only two points.

You may orient the circle as a rotation around an axis (but the sphere is not on the axis of rotation of course).

DSG

algorithmic geometry solution

For this, you need to have already solved the general intersection of two 2D circles as a 2D algorithm (able to handle all cases):

Given: SPH = [ **c** r ] CIR = [ **c orient** r ]

Solve: numIntersectionPts (0,1,2,infinite), ** i1 i2 **

1. Translate coordinates adopting CIR.c as newOrigin (compute SPH --> SPH')

SPH' = Translate(CIR.c, SPH)

2. Rotate coordinates adopting CIR.orient as the newZaxis (compute SPH' --> SPH")

Rotator R = RotatorForNewZAt(CIR.orient)

SPH" = Rotate(R, SPH')

3. Now, CIR" lies flat in x"-y" plane centered at origin.

4. Visualize 2nd circle where x"-y" plane cuts thru SPH"

if (abs(SPH".c.z) > SPH.r) numIntersectionPts = 0; return; // SPH too distant from x'=y" plane

5. Solve i1" i2" in 2D as the intersection of these two circles:

cir1 = [ 0 0 ] CIR.r

cir2 = [ SPH".c.x SPH".c.y ] sqrt(sqr(SPH.r) - sqr(SPH".c.z))

numIntersectionPts = IntersectionOf(cir1, cir2, /*returns*/ i1", i2")

6. If numIntersectionPts == 1 or 2,

back-transform solution points into original coordinates:

i1 = Untranslate(CIR.c, Unrotate(R, i1"))

i2 = Untranslate(CIR.c, Unrotate(R, i2"))

This is a general solution (works in all cases).

Illustrated solution in: *Flexing the Power of Algorithmic Geometry* by Pierre Bierre