 }
 
 // Solution parameters.
-// For better quality, use eg:
-// maxstep = 0.5
-// iterations = 200
-// maxincrease = 1.1
-const float maxstep = 2.0;     // The largest step that can be taken.
-const float minstep = 0.0001;  // The smallest step
+//#define FAST
+//#define QUALITY
+
+#if defined FAST
 const int iterations = 100;    // Maximum number of iterations
 const float maxincrease = 1.1; // Largest allowed step increase.
-
+#elif defined QUALITY
+const int iterations = 300;    // Maximum number of iterations
+const float maxincrease = 1.03; // Largest allowed step increase.
+#else
+const int iterations = 150;    // Maximum number of iterations
+const float maxincrease = 1.06; // Largest allowed step increase.
+#endif
+
+const float maxstep = 1.0;     // The largest step that can be taken.
+const float minstep = 0.001;  // The smallest step
+const float initstep = 1.0;
+const float camera = 10.0;
 const float radius = 6.0; // Restrict view to this (unused)
-
+const float horizon = 20.0; // limit on k
 void solve(float x0, float y0, float z0,
            float a, float b, float c) {
   float k0 = 0.0, k1;
   float a0 = Fun(x0,y0,z0,1.0), a1;
   bool bracketed = false;
-  float step = 1.0;
+  bool found = false;
+  float step = initstep;
   float x, y, z, w;
+  float expected = 0.0;
   for (int i = 0; i < iterations; i++) {
     if (bracketed) {
       // Once we are bracketed, just use bisection
-      if (k1-k0 < minstep) break;
+      if (k1-k0 < minstep) {
+        found = true;
+        break;
+      }
       float k2 = (k0 + k1)/2.0;
       x = x0+k2*a, y = y0+k2*b, z = z0+k2*c, w = 1.0;
       float a2 = Fun(x,y,z,w);

      }
    } else {
      k1 = k0 + step;
      //if (k1 > 20.0) break;
      x = x0+k1*a, y = y0+k1*b, z = z0+k1*c, w = 1.0;
      //if (x*x+y*y+z*z > radius*radius) break;
      a1 = Fun(x,y,z,w);
      // We can hit exactly 0 - this counts as bracketed.
      if (a0*a1 <= 0.0) {
        bracketed = true;
      } else {
        float step0 = step;
        step = a1*step/(a0-a1);
        //step = perp(k0,a0,k1,a1); // Nice idea...
        step = abs(step);
        step = min(step,maxstep);
        // Don't grow step by more than 10%

}

float Chmutov14(float x, float y, float z, float w) {
  return T14(x)+T14(y)+T14(z)+1.0;
}

float Chmutov18(float x, float y, float z, float w) {
