   __shared__ float dq[MAXDEPTH * THREADS5/WARPSIZE];
 
   if (0 == threadIdx.x) {
-    tmp = radiusd * 2;
+    // tmp = radiusd * 2;
     // precompute values that depend only on tree level
-    dq[0] = tmp * tmp * itolsqd;
+    dq[0] = radiusd * theta; //tmp * tmp * itolsqd;
     for (i = 1; i < maxdepthd; i++) {
-      dq[i] = dq[i - 1] * 0.25f;
-      dq[i - 1] += epssqd;
-    }
-    dq[i - 1] += epssqd;
+      dq[i] = dq[i - 1] * 0.5f; // radius is halved with every level of the tree
+      // dq[i - 1] += epssqd;
+    }
+    // dq[i - 1] += epssqd;
 
     if (maxdepthd > MAXDEPTH) {
       *errd = maxdepthd;

          if (n >= 0) {
            dx = posxd[n] - px;
            dy = posyd[n] - py;
            tmp = dx*dx + (dy*dy + epssqd); // distance squared + softening (how does the softening actually interact with things?
            // tmp = dx*dx + (dy*dy + epssqd) (why softening?)
            if ((n < nbodiesd) || __all(tmp >= dq[depth])) {  // check if all threads agree that cell is far enough away (or is a body)
            //   tmp = rsqrtf(tmp);  // compute distance

    float norm;
    
    // These variables currently govern the tolerance (whether it recurses on a cell)
    float theta = 0.5f;
    InitializationKernel<<<1, 1>>>(thrust::raw_pointer_cast(errl.data()));
    gpuErrchk(cudaDeviceSynchronize());
    
