   thrust::device_vector<float> dLambdas(NY * NX);
   thrust::device_vector<float> dIntensity(50 * 300 * 1200);
   thrust::host_vector<float> hIntensity(50 * 300 * 1200);
+  thrust::device_vector<float> dIntensity2(50 * 300 * 1200);
+  thrust::host_vector<float> hIntensity2(50 * 300 * 1200);
 
   idealsphere spherer(dPhotons, dLambdas);
 

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

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

      dPhotons.assign(photonVals.data(), photonVals.data() + NY * NX);
      dLambdas.assign(lambdaVals.data(), lambdaVals.data() + NY * NX);
      
      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;
      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 / 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);
    }
}
