 
   idealsphere spherer(dPhotons, dLambdas);
 
-  for (int img = 0; img < 2779; img++)
+  for (int img = 0; img < fullsize[0]; img++)
     {
       offset[0] = img;
       lambdaSpace.selectHyperslab(H5S_SELECT_SET, count, offset);
+      expectedLambdaSpace.selectHyperslab(H5S_SELECT_SET, count, offset);
       photonSpace.selectHyperslab(H5S_SELECT_SET, count, offset);
       lambdas.read(lambdaVals.data(), H5::PredType::NATIVE_FLOAT, memSpace, lambdaSpace);
       photons.read(photonVals.data(), H5::PredType::NATIVE_SHORT, memSpace, photonSpace);
-
       int psum = 0;
       double lsum = 0;
       for (int y = 0; y < NY; y++)
 	{
 	  for (int x = 0; x < NX; x++)
 	    {
-	      if (((y > 195 && y < 231) || y < 36 || ((x < 255  || (x < 300 && (y > 124 && y < 325))) && (y > 92 || x > 170 || x < 90)) && !(x<55 && y > 374))
+/*	      if (((y > 195 && y < 231) || y < 36 || ((x < 255  || (x < 300 && (y > 124 && y < 325))) && (y > 92 || x > 170 || x < 90)) && !(x<55 && y > 374))
 	      || y < 105 || y > 343)
+	      {
+		photonVals[y][x] = 0;
+		lambdaVals[y][x] = 0;
+	      }*/
+	      if (y < 233 || (x < 300 && y < 314) || x < 255)
+//	      if (y < 233 || (x < 350 && y < 370) || x < 255)
 	      {
 		photonVals[y][x] = 0;
 		lambdaVals[y][x] = 0;

{
  idealsphere& spherer;
  float factor;
__host__ __device__ likelihood(idealsphere& spherer, float factor) : spherer(spherer), factor(factor) {}

  __host__ __device__ float operator () (const int& data) {
    if (!spherer.dLambdas[data]) return 0;

    float intensity = spherer(data) * factor + spherer.dLambdas[data] * spherer.lfactor;
    float val = 0;
    if (spherer.dPhotons[data])
      {
	val += spherer.dPhotons[data] * log(intensity);
      }
    val -= intensity;

    return val;
  }
};


__global__ void computeintensity(float* target, float* intens, idealsphere myspherer, float psum, float lsum)
{
  myspherer.r = exp(myspherer.roffset + (threadIdx.z + blockIdx.z * blockDim.z) * myspherer.rfactor);
  myspherer.offsetx = (threadIdx.x) * 3 + myspherer.baseoffsetx;
  myspherer.offsety = (threadIdx.y * 3 + myspherer.baseoffsety);
  myspherer.lfactor = 1.0 / sqrt(lsum) * (((int) blockIdx.x)) + 1.0;
//    myspherer.lfactor = 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;
  float intensity = thrust::reduce(thrust::seq,
				   thrust::make_transform_iterator(thrust::make_counting_iterator(0), myspherer),
				   thrust::make_transform_iterator(thrust::make_counting_iterator(NY * NX), myspherer));

	}
	if (hIntensity[k] < minval) minval = hIntensity[k];
}

      dim3 unigrid(1,1,1);
      dim3 uniblock(1,1,1);
      getExpLambdas<<<unigrid, uniblock>>>(dExpLambdas.data().get(), dLogLs.data().get(), maxidx, spherer, psum, lsum, block, grid);
      expLambdaVals.assign(dExpLambdas.begin(), dExpLambdas.end());
      logLsVals.assign(dLogLs.begin(), dLogLs.end());
      expectedLambdas.write(expLambdaVals.data(), H5::PredType::NATIVE_FLOAT, memSpace, expectedLambdaSpace);
      logLs.write(logLsVals.data(), H5::PredType::NATIVE_FLOAT, memSpace, expectedLambdaSpace);
      //expectedLambdaFile.flush(H5F_SCOPE_LOCAL);
      hIntensity2.assign(dIntensity2.begin(), dIntensity2.end());
	int maxR = maxidx / grid.y / block.y / grid.x / block.x;
	int maxX = maxidx % block.x;
	int maxY = (maxidx / (grid.x * block.x)) % block.y;
