 //    if (rawx < 233 - extrax * 3) rawx += /*29*/ -18 + extrax * 3;
 //    if (rawy < 226 - extray * 2) rawy += /*3*/ -12 + extray * 2;*/
 
-    float x = rawx - offsetx;
-    float y = rawy - offsety;
-    
-    float q = sqrt(x * x + y * y);
-float val;
-    float den = (2 * ((float) CUDART_PI_F) * q * r);
-    if (r > 0.47e-4)
-{
-    if ( r < 1e-4)
-{
-    val = 3 * den;
-    if (r < 0.5e-4)
-    {
-	val *= den;
-    }
-    if (r < 0.48e-4)
-    {
-       val /= sqrt(den);
-    }
-}
-else
-    val = 3 * (sinpif(2 * q * r)  - 2 * ((float) CUDART_PI_F) * q * r * cospif(2 * q * r));
-
-    den = den * den * den;
-    
-    val /= den;
-}
-else
-val = 1.0f;
-
-// Global spherical convolution 
-if (false)
-{   
-    float r = tid * 1e-4;
-    den = (2 * ((float) CUDART_PI_F) * q * r);
-    val *= 3 * (sinpif(2 * q * r)  - 2 * ((float) CUDART_PI_F) * q * r * cospif(2 * q * r));
-
-    den = den * den * den;
-    val /= den;
-}
-// Glocal circular convolution
-if (false)
-{   
-    float r = tid * 1e-4;
-    den = (2 * ((float) CUDART_PI_F) * q * r);
-    val *= 2 * j1(den);
-    
-    val /= den;
-}
-// Global normal convolution
-/*{
-    float sigma = NX * 0.1 * tid;
-    val *= exp(-q*q/(2*sigma*sigma));
-}*/
-
-
-    val = val * val;
-
-    return val;
+    int x = rawx - offsetx;
+    int y = rawy - offsety;
+    if (x < 0) x += FX;
+    if (y < 0) y += FY;
+
+    return dPattern[y * FX + x];
   }
 };
 

#include <boost/multi_array.hpp>
#include <H5Cpp.h>

#include <thrust/execution_policy.h>
#include <thrust/device_vector.h>
#include <thrust/copy.h>

const int NX = 414;
const int NY = 414;

const int FX = 2048;
const int FY = 2048;

struct realsphere : public thrust::unary_function<int, cufftComplex>
{
  float r;
  float sigma;
  float x;
  float y;

  __host__ __device__ cufftComplex operator () (const int& data) {
    int rawx = data & (FX - 1);
    int rawy = data / FX;
    float qx = rawx - x - FX * 0.5;
    float qy = rawy - y - FY * 0.5;

    float r2 = r * r - qx * qx + qy * qy;
    if (r2 < 0) r2 = 0;

    r2 = sqrt(r2);
    
    qx += x;
    qy += y;
    float gr = qx * qx + qy * qy;
    cufftComplex val;
    val.y = 0;
    val.x = r2 * exp(- gr / (2 * sigma * sigma)); 
    return val;
  }
};

struct absfunc : public thrust::unary_function<cufftComplex, float>
{
  __host__ __device__ float operator () (const cufftComplex& val) {
      return /* sqrt*/ (val.x * val.x + val.y * val.y);
  }
} sqabser;

struct idealsphere : public thrust::unary_function<int, float>
{
  int offsetx, offsety;
  float baseoffsetx, baseoffsety;
//  float arraybase;
  float r;
  float lfactor;
  float ebase;
  int tid;
  int extrax;
  int extray;
  thrust::device_ptr<short> dPhotons;
  thrust::device_ptr<float> dPattern;
  thrust::device_ptr<float> dLambdas;

  idealsphere(thrust::device_vector<short>& dPhotons, thrust::device_vector<float>& dLambdas) :
