 		int particleIdx = threadIdx.x;
 
 		float3 pos = make_float3(seed.x, seed.y, 0);
-		int numSteps = 100000;
+		int numSteps = 1000;
 
 		float3 v = make_float3(0, 0, 0);
 		float valid;

		}
	}

	// Gradient Kernels
	// ~~~~~~~~ Cut with right Particle Cloud ~~~~~~~~ //
	__global__ void CutX(float* cuts, float2* positions)
	{
		// Enough memory for this?
		float2 reference[1024];
		int idx = blockIdx.x + blockIdx.y * Width;
		idx *= blockDim.x;
		int idxR = idx + blockDim.x;
		idx += threadIdx.x;

		reference[threadIdx.x] = positions[threadIdx.x];
		__syncthreads();
		
		
		// Here comes the hammer methode!
		float2 pos = positions[idx];
		unsigned int sum = 0;
		for (int ref = 0; ref < blockIdx.x; ++ref)
		{
			// Is the reference particle within the same "cell"?
			if (abs((pos.x - reference[idxR].x) * (pos.y - reference[idxR].y)) < 0.25f)
				sum++;
			//positions[idx].x = referenceMap[(int)(positions[idx].x * CellToSeedRatio + 0.5) + (int)(positions[idx].y * CellToSeedRatio + 0.5) * WidthCells]; //referenceMap[blockIdx.x*CellToSeedRatio + blockIdx.y * CellToSeedRatio*WidthCells];//
		}
		__syncthreads();
		// Use the same buffer we had before. Reference particles are not needed anymore.

		reference[threadIdx.x].x = (float)sum;
		__syncthreads();
		// Reduce.
		for (int nextSize = blockDim.x / 2; nextSize > 0; nextSize /= 2)
		{
			if (threadIdx.x < nextSize)
				reference[idx].x += reference[idx + nextSize].x;
			__syncthreads();
		}

		// Write data to texture.
		if (threadIdx.x == 0)
		{
			//surf2Dwrite(0.1f, cuts, positions[idx].x * sizeof(float), positions[idx].y, cudaBoundaryModeTrap);
			//surf2Dwrite(reference[0].x / blockDim.x, cuts, blockIdx.x*sizeof(float), blockIdx.y, cudaBoundaryModeTrap);
			cuts[blockIdx.x + blockIdx.y * (Width - 1)] = reference[0].x;
		}
	}
	// ~~~~~~~~ Cut with upper Particle Cloud ~~~~~~~~ //
	__global__ void CutY(float* cuts, float2* positions)
	{

	}

	__global__ void StoreXY(cudaSurfaceObject_t gradsX, cudaSurfaceObject_t gradsY, float* cutsX, float* cutsY)
	{
		int px = blockIdx.x * blockDim.x + threadIdx.x;
		int py = blockIdx.y * blockDim.y + threadIdx.y;
		int linIdx = py * Width + px;
		if (px < Width-1 && py < Height-1)
		{
			surf2Dwrite(cutsX[linIdx], gradsX, px*sizeof(float), py, cudaBoundaryModeTrap);
			surf2Dwrite(cutsY[linIdx], gradsY, px*sizeof(float), py, cudaBoundaryModeTrap);
		}
	}
}
				v.x += gauss.x;
				v.y += gauss.y;

				// Critical point?
				float vLen = v.x*v.x + v.y*v.y + 1;
				vLen = sqrt(vLen);

