p: the point of origin of the ray.

rd: ray direction.

si: OUTPUT, the closest sphere hit.

t is returned, the distance to the closest intersection point.

Assumption: spheres are passed in as vec4s into a uniform, in the normal GLSL manner. Each vec4 is (x,y,z,r2): centre and radius squared.

uniform vec4 sph[4]; // in this case 4 spheres

float isect (in vec3 p, in vec3 rd, out vec4 si){

float t=999.9,tnow,b,disc;

for (int i=0; i<4; i++) { //each sphere

tnow = 9999.9; //temporary closest hit is far away

//next 4 lines are the intersect routine for a sphere

vec3 sd=sph[i].xyz-p;

b = dot ( rd,sd );

disc = b*b + sph[i].w - dot ( sd,sd );

if (disc>0.0) tnow = b - sqrt(disc);

// hit, so compare and store if this is the closest

if ((tnow>0.0001)&&(t>tnow)) {t=tnow; si=sph[i];}

}

return t;

}

A couple of explanations:

disc = b*b + sph[i].w - dot ( sd,sd );

sph[i].w is already a radius squared so saves some code here. Also dot(sd,sd) is short for length(sd) - less bytes.

if ((tnow>0.0001)&&(t>tnow))...

tnow is the current intersection parameter, the distance along rd of the current intersection. If tnow is very small, we may be intersection with the sphere we are already on and must ignore it. 0.0001 is arbitrary.

Of course to get the actual intersection you need to do this in the main routine:

if (t!=999.9) inter=p+t*rd;

or something similar.