 intens[idx] = 0;*/
 intensityfactor = 1e-6;
 }
-  intensityfactor *= pow(1.01, blockIdx.y - gridDim.y * 0.5);
+  //intensityfactor *= pow(1.005, blockIdx.y - gridDim.y * 0.5);
   likelihood likelihooder(myspherer, intensityfactor);
+
+  return likelihooder;
+}
+
+__global__ void getExpLambdas(float* target, float* target2, int maxidx, idealsphere myspherer, float psum, float lsum, dim3 block, dim3 grid)
+{
+  int tx = maxidx % block.x;
+  int bx = (maxidx / block.x) % grid.x;
+  int ty = (maxidx / block.x / grid.x) % block.y;
+  int by = (maxidx / block.x / grid.x / block.y) % grid.y;
+  int zval = (maxidx / block.x / grid.x / block.y / grid.y);
+  likelihood likelihooder = getObjects(myspherer, tx, ty, zval, bx, by, lsum, psum);
+  intensitygetter intensities(likelihooder);
+
+  thrust::tabulate(thrust::device, target, target + NX * NY, intensities);
+  thrust::tabulate(thrust::device, target2, target2 + NX * NY, likelihooder);
+}
+
+__global__ void computeintensity(float* target, float* intens, idealsphere myspherer, float psum, float lsum)
+{
+  likelihood likelihooder = getObjects(myspherer, threadIdx.x, threadIdx.y, threadIdx.z + blockIdx.z * blockDim.z, blockIdx.x, blockIdx.y, lsum, psum);
   float likelihood1 = thrust::reduce(thrust::seq,
 				   thrust::make_transform_iterator(thrust::make_counting_iterator(0), likelihooder),
 				   thrust::make_transform_iterator(thrust::make_counting_iterator(NY * NX), likelihooder));
 
+  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;	
   target[idx] = likelihood1;
   intens[idx] = likelihooder.factor;
 }

  
  hsize_t count[3] = {1, NY, NX};
  hsize_t offset[3] = {0, 0, 0};
  lambdaSpace.selectHyperslab(H5S_SELECT_SET, count, offset);
  photonSpace.selectHyperslab(H5S_SELECT_SET, count, offset);
  H5::DataSpace memSpace(3, count);

  boost::multi_array<float, 2> lambdaVals(extents[NY][NX]);
  boost::multi_array<short, 2> photonVals(extents[NY][NX]);
  thrust::device_vector<short> dPhotons(NY * NX);
  thrust::device_vector<float> dLambdas(NY * NX);
  thrust::device_vector<float> dIntensity(175000000);
  thrust::host_vector<float> hIntensity(175000000);
  thrust::device_vector<float> dIntensity2(175000000);

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

  __host__ __device__ float getIntensity(const int& data) {
  	float intensity = spherer(data) * factor + spherer.dLambdas[data] * spherer.lfactor;

	return intensity;
  }

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

    float intensity = getIntensity(data);
    float val = 0;
    if (spherer.dPhotons[data])
      {
      val += log(1.0f - exp(-intensity));
      }
      else
    val -= intensity;

    return val;
  }
};

struct intensitygetter : public thrust::unary_function<int, float>
{
  likelihood& likelihooder;

  __host__ __device__ intensitygetter(likelihood& likelihooder) : likelihooder(likelihooder) {}

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

    return likelihooder.getIntensity(data);
  }
};

__device__ likelihood getObjects(idealsphere& myspherer, unsigned int tx, unsigned int ty, unsigned int zval, unsigned int bx, unsigned int by, float lsum, float psum)
{
  myspherer.r = exp(myspherer.roffset + (zval) * myspherer.rfactor);
  myspherer.offsetx = (tx) * 3 + myspherer.baseoffsetx;
  myspherer.offsety = (ty * 3 + myspherer.baseoffsety);
  myspherer.lfactor = 1.0 / sqrt(lsum) * (((int) bx)) + 1.0;
//    myspherer.lfactor = 1.0;

  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));
