 {
   myspherer.r = exp(myspherer.roffset + (threadIdx.z + blockIdx.z * blockDim.z) * myspherer.rfactor);
   myspherer.offsetx = (threadIdx.x + blockIdx.x * blockDim.x) + myspherer.baseoffsetx;
-  myspherer.offsety = (blockIdx.y + myspherer.baseoffsety);
+  myspherer.offsety = /*(blockIdx.y +*/ myspherer.baseoffsety/*)*/;
   myspherer.lfactor = 1.0 / sqrt(lsum) * threadIdx.y + 1.0;
 
   int idx =  (threadIdx.x + blockIdx.x * blockDim.x)  + (threadIdx.y + blockIdx.y * blockDim.y) * gridDim.x * blockDim.x + (threadIdx.z + blockIdx.z * blockDim.z)  * gridDim.x * blockDim.x * gridDim.y * blockDim.y;

	}
	if (hIntensity[k] < minval) minval = hIntensity[k];
}
	int maxR = maxidx / 40 / 50;
	int maxX = maxidx % 50;
	int maxY = (maxidx / 50) % 40;

      printf("%d %d %lf %g %g %g %d %d %d %d %g\n", img, psum, lsum, minval, maxval, hIntensity[0], maxR, maxX, maxY, cudaGetLastError(), maxint);
      fflush(stdout);

	    }
	}
      
      dim3 grid(1, 400, 24);
      dim3 block(1, 7, 50);

      spherer.rfactor = 0.005;
      spherer.roffset = -10;
      spherer.baseoffsetx = NX / 2 + 70 - 0.5; // good val -10
      spherer.baseoffsety = NY / 2 + 70 - 0.5; // good val +10
      dPhotons.assign(photonVals.data(), photonVals.data() + NY * NX);
      dLambdas.assign(lambdaVals.data(), lambdaVals.data() + NY * NX);
      
