       dPhotons.assign(photonVals.data(), photonVals.data() + NY * NX);
       dLambdas.assign(lambdaVals.data(), lambdaVals.data() + NY * NX);
       
-      computeintensity<<<grid, block>>>(dIntensity.data().get(), spherer, psum, lsum);
+      computeintensity<<<grid, block>>>(dIntensity.data().get(), dIntensity2.data().get(), spherer, psum, lsum);
       hIntensity.assign(dIntensity.begin(), dIntensity.end());
+      hIntensity2.assign(dIntensity2.begin(), dIntensity2.end());
 
       float minval = 1e30;
 float maxval = -1e30;
       int maxidx = 0;
-      for (int k = 0; k < 280 * 50 * 1200; k++)
+      float maxint = 0;
+      for (int k = 0; k < 40 * 50 * 1200; k++)
 {
 	if (hIntensity[k] > maxval)
 	{
 		maxval = hIntensity[k];
+		maxint = hIntensity2[k];
 		maxidx = k;
 	}
 	if (hIntensity[k] < minval) minval = hIntensity[k];
 }
-	int maxR = maxidx / 280 / 50;
+	int maxR = maxidx / 40 / 50;
 	int maxX = maxidx % 50;
-	int maxY = (maxidx / 50) % 280;
+	int maxY = (maxidx / 50) % 40;
 
-      printf("%d %d %lf %g %g %g %d %d %d %d\n", img, psum, lsum, minval, maxval, hIntensity[0], maxR, maxX, maxY, cudaGetLastError());
+      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);
     }
 }

  thrust::device_vector<float> dLambdas(NY * NX);
  thrust::device_vector<float> dIntensity(50 * 300 * 1200);
  thrust::host_vector<float> hIntensity(50 * 300 * 1200);

  idealsphere spherer(dPhotons, dLambdas);


	}
      
      dim3 grid(25, 40, 24);
      dim3 block(2, 1, 50);

      spherer.rfactor = 0.005;
      spherer.roffset = -10;
