                 root11 = (-b - sdisc);
             }
         }
-        */
 
         bool check_second = true;
         if( rtPotentialIntersection( root1 + root11 ) ) {

rtDeclareVariable( legion::LocalGeometry, lgeom, attribute local_geom, ); 
rtDeclareVariable( optix::Ray,            ray,   rtCurrentRay, );




RT_PROGRAM void sphereIntersect( int )
{
    // TODO: rename these origin, direction
    float3 O = ray.origin - center;
    float3 D = ray.direction;

    float b = optix::dot(O, D);
    float c = optix::dot(O, O)-radius*radius;

        float sdisc = sqrtf(disc);
        float root1 = (-b - sdisc);


        float root11 = 0.0f;

        // refine root1
        if( fabsf(root1) > 10.f * radius )
        {
            float3 O1 = O + root1 * direction;
            b = optix::dot(O1, D);
            c = optix::dot(O1, O1) - radius*radius;
            disc = b*b - c;

            if(disc > 0.0f)
            {
                sdisc = sqrtf(disc);
                root11 = (-b - sdisc);
            }
        }

        float  t = root1 + root11;
        if( t >= tmin && t <= tmax  ) {

            const float3 normal = (O + t*D)/radius;

            // Fill in a localgeometry
            lg.position         = origin + t*direction;
            lg.position_object  = lg.position; 
            lg.geometric_normal = normal;
            lg.shading_normal   = normal;
            lg.texcoord         = make_float2( 0.0f );

            return true; 
        } 

        float root2 = (-b + sdisc) +  root1;
        t = root2;
        if( t >= tmin && t <= tmax  ) {

            const float3 normal = (O + t*D)/radius;

            // Fill in a localgeometry
            lg.position         = origin + t*direction;
            lg.position_object  = lg.position; 
            lg.geometric_normal = normal;
            lg.shading_normal   = normal;
            lg.texcoord         = make_float2( 0.0f );

            return true; 
        }
    }
    return false;
}


RT_PROGRAM void sphereIntersect( int )
{
    float3 O = ray.origin - center;
    float3 D = ray.direction;

    float b = optix::dot(O, D);
    float c = optix::dot(O, O)-radius*radius;
    float disc = b*b-c;
    if(disc > 0.0f)
    {
        float sdisc = sqrtf(disc);
        float root1 = (-b - sdisc);


        float root11 = 0.0f;

        // refine root1
        if( fabsf(root1) > 10.f * radius )
        {
            float3 O1 = O + root1 * ray.direction;
