Ray-Torus Intersection

int Torus::line_intersect(float ax, float ay, float az, 
			  float bx, float by, float bz,
			  int * num_intersections,
			  float * intersections,
			  float * shade) const
//              
// Intersect ray
//          [ax]     [bx]              
//	    [ay] + t [by]
//	    [az]     [bz]
// with torus at origin, major radius R, minor radius r
//
{

  // This struct is syntactic sugar so that a.b,
  // (the dot product) looks just right (:-)
  struct { float a,b; } a;
  a.b = ax*bx + ay*by + az*bz;
  a.a = ax*ax + ay*ay + az*az;

  // Set up quartic in t:
  //
  //  4     3     2
  // t + A t + B t + C t + D = 0
  //
  float R2 = R*R;
  float K = a.a - r*r - R2;
  float A = 4*a.b;
  float B = 2*(2*a.b*a.b + K + 2*R2*bz*bz);
  float C = 4*(K*a.b + 2*R2*az*bz);
  float D = K*K + 4*R2*(az*az - r*r);

  // Solve quartic...
  double roots[4];
  int nroots = solve_quartic(A,B,C,D,roots,SOLVE_ALL);

  *num_intersections = 0;
  while(nroots--)
    {
      float t = roots[nroots];
      float x = ax + t*bx;
      float y = ay + t*by;
      float l = R*(M_PI/2 - atan2(y,x));
      if (l <= vlength && l >= 0)
        intersections[(*num_intersections)++] = t;
    }
}