         float sdisc = sqrtf(disc);
         float root1 = (-b - sdisc);
 
-        bool do_refine = false;
 
         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;

        } 

        if(check_second) {
            float root2 = (-b + sdisc) + (do_refine ? root1 : 0);
            if( rtPotentialIntersection( root2 ) ) {
                const float  t      = root2; 
                const float3 normal = (O + t*D)/radius;


                root11 = (-b - sdisc);
            }
        }

        bool check_second = true;
        if( rtPotentialIntersection( root1 + root11 ) ) {
