Circle and Sphere Equation
I am trying to find an equation for the intersection of a
circle and a sphere in 3 space. Here's the situation.
I'm working on a computer program for drawing free flowing drapery in real
time. I need to compute the 2 points on a plane that passes through the
origin that are 1 unit distance from the origin and a second "center" not
on the plane. I know the possibilities are that either there is 0, 1 or 2
intersects, but because of the rest of my algorithm I'm guaranteed that
there will always be 2 intersections.
I've conceived the problem in 2 different ways. First I thought of the
intersects as being the intersection of a unit radius circle on my plane
and sphere around the center not on my plane. I solves for x,y and z in my
plane and substituted these into a sphere at the origin to give me a
formula for the circle on the plane. Then I'd solve for the x,y and zs of
the sphere and substitute them into the circle. When I did so the algebra
blew up into a horrendous problem with no solution in sight.
I also conceived of the problem as 2 unit radius spheres, one at the
origin one at the other center that intersect in a circle, then this
circle is intersected by my plane. It also became very cumbersome.
I know I could probably solve this if I rotate my centers into alignment
with my 3D space, then rotate it back when I'm finished. I'm avoiding this
solution as being to slow or inefficient to get the job done for a real
time program. I need to solve this for thousands of points in a tiny
fraction of second. That's why I'm hoping to find a straight forward
formula for the x, y and z of the 2 intercepts. Any ideas?
Well, this is a very interesting question, and I am inclined to spend a lot
more time working on it than I ought to. Given that you're a computer
programmer and I'm a scientist, you're probably paid a lot more than I am.
However, I know from my own programming experience that the actual
mathematical guts of a program take up only about 5% of the programming
work, and the remaining 95% of the effort goes into I/O. So I'll give you
an idea of how I'd address the problem, and leave it to you to work it out.
What you have is a fairly simple geometric problem made fairly opaque by its
representation in Cartesian coordinates. You're trying to figure out the
problem in the language of algebra, which is how you will eventually need to
communicate it to your computer. But first, it would be best for you to
think of it in terms of shapes, distances, and spaces, which is, after all,
how you eventually want your computer to communicate it back to its user.
You are trying to avoid rotating the coordinates into alignment, because you
think that transforming coordinates may be too computation intensive. True,
transforming coordinates can be slow, especially if you have many different
coordinate systems to address (one for each point, I imagine). But in this
case you don't need to use trigonometric or exponential functions, which are
the real hogs of CPU time, and in fact you won't actually need to transform
coordinates in the code at all. You may still, however, be served well by
at least thinking in terms of transformed coordinates when you set up your
Your system is fairly simple. (However, you might want to sketch along on a
piece of paper to follow the discussion, as I'm not including pictures in my
explanation.) You have two points in 3-space, A and B, and a plane that
passes through A, but not through B. Around both A and B are spheres of
unit radius. You already know that the distance |AB| < 2 (otherwise their
spheres wouldn't intersect), and that the plane passes through the circle of
intersection between these two spheres. (Here I am denoting the length of
vector x as |x|.) You want to find the coordinates of the two points of
intersection between the plane and this circle.
My recommendation is to first find the circle of intersection between the
plane and the sphere around point B. (Don't panic! You don't need to write
an x-y equation for it!) The first thing to do is express the vector
between points A and B as a sum of two vectors, one of which is in the plane
of your problem, and the other of which is orthogonal to the plane. (Bear
with me, it will eventually become clear why I am doing this.) The vector
in the plane points from point A to point B', which is the point in the
plane nearest to point B. (Point B' is the projection of point B onto the
plane; vector AB' is the projection of vector AB onto the plane.)
You can find the vector B'B by first finding the unit vector perpendicular
to the plane, k. If you have expressed the plane as an equation ax + by +
cz = d, then a (non-unit) vector perpendicular to this plane is (a,b,c). The
unit vector then is k = (a,b,c)/|(a,b,c)|. Once you have k, it is very
simple to find the vector B'B: it is B'B = (k dot AB)k, where "x dot y"
denotes the inner (dot) product between vectors x and y, e.g.,
(x1,x2,x3)dot(y1,y2,y3) = x1y1 + x2y2 + x3y3. (This is a very simple
Once you have vector B'B, all you must do to find AB' is subtract it from
AB. Thus, AB' = AB - B'B.
The shortest distance between point B and the plane is the length of the
vector B'B. Point B' is the center of the circle defined by the
intersection between the plane and the sphere about B. The radius r of this
circle can be found by the Pythagorean Theorem: |B'B|^2 + r^2 = 1. (You
can use this formula because the vectors AB' and B'B are perpendicular.)
Now your problem consists of two points in a plane, A and B', a circle of
radius 1 about A, and a circle of radius r about B'. You want to find the
points of intersection between these two circles. This problem reduces to
finding the position of the third vertex of a triangle in which one side is
vector AB', with length |AB'|, and the other two sides have lengths 1 and r.
To get this in a manner that makes sense, it is wise to define a new
coordinate system. RELAX! This won't be very computation intensive. You
want your coordinate system to be such that point A lies at (0,0) and point
B' lies at (a,0), where a = |AB'|, the distance from A to B'. This means
that the x-axis goes along the vector AB'. The y-axis of your new
coordinates must be perpendicular to the x-axis, and must fall within the
plane that defines your problem. Thus, the x and y axes of this coordinate
system define the plane or your problem.
It will be computationally easiest if you find the orthonormal basis vectors
i and j of this coordinate system. In the new coordinates, they are simply
i = (1,0) and j = (0,1). In the old coordinates, they are i = AB'/|AB'|,
and j must be defined as the unit vector perpendicular to i that falls
within the plane of your problem. This is easy to find. It is just the
cross product of vectors i and k (recall that k is the unit vector
perpendicular to your plane): j = k x i. (The "x" here denotes vector cross
product, not the letter x. Everywhere else in this discussion it is the
Why did I have you find i and j? Because you will use them to convert
answers from the simple coordinate system into the original one. In the
new, simple coordinate system, your geometry problem is this: Find the two
points of intersection between the circle of radius 1 about (0,0) and the
circle of radius r about (a,0). This is done by solving the two
x^2 + y^2 = 1 and
(a-x)^2 + y^2 = r^2
which gives the easy answers of x = 0.5*(a + (1 - r^2)/a) and y^2 = 1 - x^2.
(We are operating in these transformed coordinates so that this algebra is
easy.) The points of intersection are A + ix + jy and A + ix - jy. This is
true in both coordinate systems. So, without thinking any more about
converting between coordinates, you can just apply these formulas to get the
points you need. The unit vectors i and j will have fewer zeroes in their
components in the original coordinate system, but that will be no trouble at
all for your computer. You should also be happy to know that the toughest
mathematical operation you need to use is a square root.
If this doesn't answer your question, or if it's not clear to you, ask again
if you want.
Richard E. Barrans Jr., Ph.D.
PG Research Foundation, Darien, Illinois
Click here to return to the Mathematics Archives
Update: June 2012