/*-----------------------------------------------------------------------------------------
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
 *
 * This file is part of IFOS.
 *
 * IFOS is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
 *
 * IFOS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
 -----------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
 *   2D-FFT preparation
 *  ----------------------------------------------------------------------*/

#include "fd.h"

void	fft2(float **array, float **arrayim, int NYG, int NXG, int dir) {

		extern FILE *FP;
		
		COMPLEX	*C;
		int	j, i, k, nx, ny, nxy;
		float **RE, **IM;
	
		fprintf(FP,"\n Testposition fft2_filt 1\n ");
		usleep(900000);
	
		/****************************************************/
		/* recompute the dimensions to become values of 2^m */
		nx = (int)(pow(2,(ceil(log(NXG)/log(2)))));
		ny = (int)(pow(2,(ceil(log(NYG)/log(2)))));
		nxy = max(nx,ny);
		nxy *= 2;
		
		fprintf(FP,"\n Testposition fft2_filt 1, nxy: %i\n ",nxy);
		usleep(900000);
		
		RE = matrix(1, nxy, 1, nxy);
		IM = matrix(1, nxy, 1, nxy);
		C = cplxvector(1, nxy*nxy);
		

				
		/* copy array data to the temp-array */
		for (i=1;i<=NXG;i++)
			for (j=1;j<=NYG;j++){
			  RE[j][i] = array[j][i];
			  IM[j][i] = arrayim[j][i];
			}
			
		fprintf(FP,"\n Testposition fft2_filt 3\n ");
		usleep(900000);
	
	if (dir==1) {
		k = 1;
		for (i=1;i<=nxy;i++)
			for (j=1;j<=nxy;j++) {
				C[k].re = RE[j][i];
				C[k].im = 0.0;
				k++;
			}
		forward_fft2f(C, nxy, nxy);
		k = 1;
		for (i=1;i<=nxy;i++)
			for (j=1;j<=nxy;j++) {
				RE[j][i] = C[k].re;
				IM[j][i] = C[k].im;
				k++;
			}
	} else {
		k = 1;
		for (i=1;i<=nxy;i++)
			for (j=1;j<=nxy;j++) {
				C[k].re = RE[j][i];
				C[k].im = IM[j][i];
				k++;
			}
		inverse_fft2f(C, nxy, nxy);
		k = 1;
		for (i=1;i<=nxy;i++)
			for (j=1;j<=nxy;j++) {
				RE[j][i] = C[k].re;
				IM[j][i] = 0.0;
				k++;
			}
	}
	/* write result back into array*/
	for (i=1;i<=NXG;i++)
		for (j=1;j<=NYG;j++){
			  array[j][i]=RE[j][i];
			  arrayim[j][i]=IM[j][i];
		}
				
	free_cplxvector(C, 1, nxy*nxy);
	free_matrix(RE, 1,nxy,1,nxy);
	free_matrix(IM, 1,nxy,1,nxy);
}



/********************************************************/
/* 2D-FFTSHIFT						*/
/********************************************************/
// void	fft2shift(float **RE, float **IM, int ny, int nx) {
// 
// 	int	j, i;
// 	float	**A;
// 		
// 	
// 	A = matrix(1, ny, 1, nx);
// 	
// 	/* shift the real part */
// 	for (i=1;i<=nx;i++)
// 		for (j=1;j<=ny;j++)	A[j][i] = RE[j][i];
// 	for (i=1;i<=nx/2;i++)
// 		for (j=1;j<=ny/2;j++) {
// 			RE[j][i]		= A[j+ny/2][i+nx/2];
// 			RE[j+ny/2][i+nx/2]	= A[j][i];
// 			RE[j+ny/2][i]		= A[j][i+nx/2];
// 			RE[j][i+nx/2]		= A[j+ny/2][i];
// 		}
// 
// 	/* shift the imaginary part */
// 	for (i=1;i<=nx;i++)
// 		for (j=1;j<=ny;j++)	A[j][i] = IM[j][i];
// 	for (i=1;i<=nx/2;i++)
// 		for (j=1;j<=ny/2;j++) {
// 			IM[j][i]		= A[j+ny/2][i+nx/2];
// 			IM[j+ny/2][i+nx/2]	= A[j][i];
// 			IM[j+ny/2][i]		= A[j][i+nx/2];
// 			IM[j][i+nx/2]		= A[j+ny/2][i];
// 		}
// 	
// 	free_matrix(A, 1, ny, 1, nx);
// }