////////////////////////////////////////////////////////////
//
//  Sphere::intersect()
//
//  (c) Copyright 2001
//  David C. Banks


////////////////////////////////////////////////////////////
//
//  Compute the sphere's intersection point when given
//  a ray (defined by point p0 and direction dp).
// 
//  Return the number of intersections.
//


////////////////////////////////////////////////////////////
//  The ray is parameterized according to the vector equation
//
//    p(t)  =  p0  +  t dp
//
//  so r(0) is the point p0 defining the ray, and r(1) is 
//  the other point defining the ray.
//
//  A point p on a sphere satisfies the equation
//
//    || p(t) - center || = radius
//
//  Solving these simulataneously leads to a quadratic equation
//
//    t = a0 + a1 t + a2 t^2

int
Sphere::
intersect(SbVec3f &intersection,
	  const SbVec3f &p0, const SbVec3f &dp) const
  {
  float a0, a1, a2;     // coefficients:  a0 + a1 s + a2 s^2 = 0
  float d;              // discriminant of quadratic polynomial
  float sMin, s1, s2;

  //   t = a0 + a1 t + a2 t^2
  //
  // You can derive these coefficients yourself, or search the Web 
  // for ray-sphere intersection equation.

  a0 = -radius*radius + (p0-center).dot(p0-center);
  a1 = 2.0*(p0-center).dot(dp);
  a2 = dp.dot(dp);


  // Discriminant d of the quadratic equation.

  d = a1*a1 - 4.0*a0*a2;

  // Negative discriminant: imaginary roots.
  // The ray does not intersect the sphere.
  if (d < 0.0)
    {
    return (0);
    }


  // Zero discriminant: double roots
  // The ray intersects the sphere one time,
  // tangent to the sphere.

  else if (d == 0.0)
    {
    if (a2 != 0.0)
      {
      sMin = -a1/(2.0*a2);
      if (sMin > 0.0)      // Intersection lies in front.
	{
	intersection = p0 + sMin*dp;
	return (1);
	}
      else                 // Intersection does not lie in front.
	{
	return (0);
	}
      }
    else                   // Divide-by-zero error in quadratic.
      {
      return (0);
      }
    }

  // Positive discriminant: 2 distinct roots
  // The ray intersects the sphere two times.
  // 
  // Find the nearest intersection lying in front
  // of p0 (in the dp direction, with t positive).
  else // (d > 0.0)
    {
    if (a2 != 0.0)
      {
      float muchBigger = 100.0;
      if ( a1*a1 > muchBigger*4.0*a0*a2 )
	{
	// See http://mathworld.wolfram.com/QuadraticEquation.html
	// Eric Weisstein's world of Mathematics
	// "Quadratic Equation"
	float q = -0.5*(a1 + copysign(1.0, a1)*sqrt(d) );
	s1 = q / a2;
	s2 = a0 / q;
	}

      else
	{
	float sqrtd = sqrt(d);
	s1 =  1.0/(2.0*a2) * (-a1 + sqrtd);
	s2 =  1.0/(2.0*a2) * (-a1 - sqrtd);
	}

      sMin = -1.0;            // Flag value: we expect sMin to be positive.
      if ( (s1 > 0.0) && (s1 < s2) )
	{
	sMin = s1;
	}
      else if ( (s2 > 0.0) && (s2 < s1) )
	{
	sMin = s2;
	}
      if (sMin <= 0.0)
	{
	return(0);            // Neither intersection lies in front.
	}
      else
	{
	intersection = p0 + sMin*dp;
	return(2);
	}
      }
    else
      {
      return (0);             // Divide-by-zero error in quadratic.
      }
    }
  }