         if(y!=0)            { c1 = p[idx-w]; c2 = p[idx+(size_t)h*w-w]; c3 = p[idx+(size_t)2*h*w-w]; }
         if((y+1)!=h)        { d1 = p[idx+w]; d2 = p[idx+(size_t)h*w+w]; }
         if(y!=0 && x!=0)    {                                           e3  = p[idx+(size_t)2*h*w-w-1]; }
+        //if((y+1)!=h && (x+1)!=h)    {                                   e3  = p[idx+(size_t)2*h*w+w+1]; }
 
         d[idx] = sqrtf(1.0f/3.0f)*( a1 + b1 + c1 + d1 - 4*pp )
                + sqrtf(2.0f/3.0f)*( c2 + d2 - a2 - b2 )
                + sqrtf(8.0f/3.0f)*( pp + e3 - a3 - c3 );
-
-    }
-}
-
-
-__global__ void ComputeU(float *u, float *v, float *d, int w, int h, float theta)
+               //+ sqrtf(8.0f/3.0f)*( pp + e3 - b3 - d3 );
+    }
+}
+
+
+__global__ void ComputeU(float *v, float *d, float *u, int w, int h, float theta)
 {
     int x = threadIdx.x + blockIdx.x*blockDim.x;
     int y = threadIdx.y + blockIdx.y*blockDim.y;

    getParam("tau", tau, argc, argv);
    cout << "tau: " << tau << endl;

    float decay = 0.95f;
    getParam("decay", decay, argc, argv);
    cout << "decay: " << decay << endl;

    int N1 = 20;
    getParam("N1", N1, argc, argv);
    cout << "N1: " << N1 << endl;

    int N2 = 10;
    getParam("N2", N2, argc, argv);
    cout << "N2: " << N2 << endl;


        if(y!=0)                    { c = u[idx-w]; }
        if((y+1)!=h)                { d = u[idx+w]; }
        if((y+1)!=h && (x+1)!=w)    { e = u[idx+w+1]; }
        //if(y!=0 && x!=0)            { e = u[idx-w-1]; }

        dd[idx]               = sqrtf(1.0f/3.0f)*( a + b + c + d - 4*u[idx] );
        dd[idx+(size_t)h*w]   = sqrtf(2.0f/3.0f)*( c + d - a - b );
        dd[idx+(size_t)2*h*w] = sqrtf(8.0f/3.0f)*( u[idx] + e - b - d );
        //dd[idx+(size_t)2*h*w] = sqrtf(8.0f/3.0f)*( u[idx] + e - a - c );
    }
}

