It's easier (I think) to do the whole thing using vectors and dot products. At least it's more compact, which makes it a little less prone to error.
First, the sphere. p is a point at (x,y,z), r is the sphere radius, pc
is the sphere center at (xc, yc, zc).
The equation for a sphere centered at the origin is
x2 + y2 + z2 = p . p = r2 or p . p - r2 = 0The equation for a sphere centered at pc (the vector version of the equation I gave in class) is
(p - pc) . (p - pc) - r2 = 0Meanwhile, we've also got a ray, in direction d = (dx, dy, dz), parameterized by t, starting at p0 = (x0, y0, z0):
p = t d + p0Plugging the ray equation into the sphere equation, we get
(t d + p0 - pc) . (t d + p0 - pc) - r2 = 0To solve, we'll need to get this into a t2 + b t + c = 0 form to apply the quadratic formula to find the two solutions t = (-b +/- sqrt(b2 - 4 a c)) / (2 a). Getting it into that form is a simple application of the distributive property of dot products:
d . d t2 + 2 d . (p0 - pc) t + (p0 - pc) . (p0 - pc) - r2 = 0Or, to make it even more explicit:
a = d . dAs a reminder of one thing I said in class, the sign of the discriminant, b2 - 4 a c, inside the square root tells you how many roots the equation has. If this is negative, there are no intersections, and the ray misses the sphere. This is the one that will happen most often. If it's zero, there's one glancing intersection right on the edge of the sphere. If it's positive, there are two intersections, one where the ray enters the sphere, and one where it exits again.
b = 2 d . (p0 - pc)
c = (p0 - pc) . (p0 - pc) - r2