public static float? RayToEllipsoid(ref Ray ray, ref BoundingEllipsoid ellipsoid)
{
    // Source: http://www.ogre3d.org/forums/viewtopic.php?f=2&t=26442&start=0
 
    Ray transformedRay = ray;
    transformedRay.Position -= ellipsoid.Center;
    transformedRay.Direction.Normalize();
 
    float a = ((transformedRay.Direction.X * transformedRay.Direction.X) / (ellipsoid.Radius.X * ellipsoid.Radius.X))
            + ((transformedRay.Direction.Y * transformedRay.Direction.Y) / (ellipsoid.Radius.Y * ellipsoid.Radius.Y))
            + ((transformedRay.Direction.Z * transformedRay.Direction.Z) / (ellipsoid.Radius.Z * ellipsoid.Radius.Z));
 
    float b = ((2 * transformedRay.Position.X * transformedRay.Direction.X) / (ellipsoid.Radius.X * ellipsoid.Radius.X))
            + ((2 * transformedRay.Position.Y * transformedRay.Direction.Y) / (ellipsoid.Radius.Y * ellipsoid.Radius.Y))
            + ((2 * transformedRay.Position.Z * transformedRay.Direction.Z) / (ellipsoid.Radius.Z * ellipsoid.Radius.Z));
 
    float c = ((transformedRay.Position.X * transformedRay.Position.X) / (ellipsoid.Radius.X * ellipsoid.Radius.X))
            + ((transformedRay.Position.Y * transformedRay.Position.Y) / (ellipsoid.Radius.Y * ellipsoid.Radius.Y))
            + ((transformedRay.Position.Z * transformedRay.Position.Z) / (ellipsoid.Radius.Z * ellipsoid.Radius.Z))
            - 1;
 
    float d = ((b * b) - (4.0f * a * c));
 
    if (d < 0)
    {
        return null;
    }
    else
    {
        d = (float)Math.Sqrt(d);
    }
 
    float hit = (-b + d) / (2.0f * a);
    float hitsecond = (-b - d) / (2.0f * a);
 
    if (hit < hitsecond)
    {
        return hit;
    }
    else
    {
        return hitsecond;
    }
}