#ifndef CUDA_CUFFT_SHIFT_H
#define CUDA_CUFFT_SHIFT_H

#define BLOCK_SIZE 16

__global__ void devCufftShift(char* gpuSource, int sE, int sX, int sY, int sZ)
{
	int ix = blockIdx.x * blockDim.x + threadIdx.x;
	int blocks_per_slice = sY/blockDim.y + 1;
	int iy = (blockIdx.y % blocks_per_slice)*blockDim.y + threadIdx.y;
	int iz = blockIdx.y/blocks_per_slice;

	if(ix >= sX || iy >= sY || iz >= sZ/2)
		return;

	//array index
	int i = iz*sX*sY*sE + iy*sX*sE + ix*sE;

	int iSwap = (iz + sZ/2)*sX*sY*sE;
	
	if(iy < sY/2)
		iSwap += (iy + sY/2)*sX*sE;
	else
		iSwap += (iy - sY/2)*sX*sE;

	if(ix < sX/2)
		iSwap += (ix + sX/2)*sE;
	else
		iSwap += (ix - sX/2)*sE;
	//iSwap = (iz + sZ/2)*sX*sY*sE + (iy + sY/2)*sX*sE + (ix + sX/2)*sE;

	char temp;
	for(int e=0; e<sE; e++)
	{
		temp = gpuSource[i + e];
		gpuSource[i + e] = gpuSource[iSwap + e];
		gpuSource[iSwap + e] = temp;
	}
	
}

void rts_cudaCufftShift(char* gpuSource, int elementSize, int sX, int sY, int sZ)
{
	//This function duplicates the Matlab fftshift function, making it easier to apply certain filters in
	//the frequency domain.

	//currently only usable for even numbers
	if(sX%2 != 0 || sY%2 != 0 || sZ%2 != 0)
	{
		printf("Error: cudaCufftShift dimensions must be even (for now)\n");
		exit(1);
	}

	dim3 block(BLOCK_SIZE, BLOCK_SIZE);
	dim3 grid(sX/block.x + 1, (sY/block.y + 1)* sZ);

	//call the kernel
	devCufftShift<<<grid, block>>>(gpuSource, elementSize, sX, sY, sZ);



}


#endif