void aapx(const Eigen::Ref<const Eigen::MatrixXcd>& uh, Eigen::Ref<const Eigen::MatrixXcd>& vh, Eigen::Ref<const Eigen::MatrixXcd>& ph){ 
//some calculations

}

static const int nx = 10;
static const int ny = 10; 
static const int mm = nx* 3/2;

Eigen::Matrix<std::complex<double>, (ny+1), (nx)> uh;
uh.setZero();

Eigen::Matrix<std::complex<double>, (ny+1), (nx)> vh;
vh.setZero();

Eigen::Matrix<std::complex<double>, (ny+1), nx> ph;  
ph.setZero();

aapx(uh,vh,ph); //ERROR 

aapx(const Eigen::Ref<const Eigen::MatrixXcd>& uh, const Eigen::Ref<const Eigen::MatrixXcd>& vh, Eigen::Ref<const Eigen::MatrixXcd>& ph){

	//define a new Eigen matrix
	Eigen::Matrix<std::complex<double>, (ny+1), mm> uhp;
	uhp.setZero();
	//pad the first input with zeros 
	uhp(all,seqN(0,nx/2)) = uh(all,seqN(0,nx/2));
	uhp(all,seq(last-nx/2+1, last)) = uh(all,seq(last+1-nx/2, last)); 

	//define a new Eigen matrix
	Eigen::Matrix<std::complex<double>, (ny+1), mm> vhp;
	vhp.setZero();
	//pad the first input with zeros
	vhp(all,seqN(0,nx/2)) = vh(all,seqN(0,nx/2)); 
	vhp(all,seq(last-nx/2+1, last)) = vh(all,seq(last+1-nx/2, last));

	//define a c++ array in order to use FFTW with Eigen matrices
	fftw_complex *uhpk; //dummy variable
	uhpk = (fftw_complex*) fftw_malloc((((mm)*(ny+1))*nyk)* sizeof(fftw_complex)); 
	memset(uhpk, 42, (((mm)*(ny+1))*nyk)* sizeof(fftw_complex));

	memcpy(uhpk,uhp.data(),uhp.size() * sizeof(fftw_complex));
	//######################################################
	//use 1dfft function (note: write in a separate function)

	cplx_buffer dummyII2 = my_fftw_allocate_cplx((ny+1), mm);
	for (int i = 0; i < mm; i++){ 
		for (int j = 0; j < ny+1; j++){ 
			for (int k = 0; k < ncomp; k++){
				dummyII2(j,i)[k] = uhpk[i*(ny+1)+j][k];
			}
		}
	
	}
	real_buffer dummyRR2 = my_fftw_allocate_real((ny+1), mm);
	//real(ifft(duxk,[],2))
	{ // Transform all the rows
	fftw_execute(fftw_plan_many_dft_c2r(1, &nx, dummyII2.rows, dummyII2.a, &dummyII2.cols, 1, dummyII2.cols, dummyRR2.a, &dummyII2.cols, 1, dummyRR2.cols, FFTW_ESTIMATE));
	//std::cout << "Inverse FT of input, row-wise\n"; 
	//print_matrixR(dummyRR2);
	}
	double *up;
	up = (double*) fftw_malloc((((ny+1)*mm))*sizeof(double)); 
	for (int i = 0; i < ny+1; i++){ 
		for (int j = 0; j < mm; j++){ 
			up[j*(ny+1)+i] = dummyRR2(i,j);
		}
	}

	fftw_complex *vhpk; 
	vhpk = (fftw_complex*) fftw_malloc((((mm)*(ny+1))*nyk)* sizeof(fftw_complex));
	memset(vhpk, 42, (((mm)*(ny+1))*nyk)* sizeof(fftw_complex)); 

	memcpy(vhpk,vhp.data(),vhp.size() * sizeof(fftw_complex));
	cplx_buffer dummyIi2 = my_fftw_allocate_cplx((ny+1), mm);

	for (int i = 0; i < mm; i++){ 
		for (int j = 0; j < ny+1; j++){ 
			for (int k = 0; k < ncomp; k++){
				dummyIi2(j,i)[k] = vhpk[i*(ny+1)+j][k];
			}
		}
	
	}
	
	real_buffer dummyRr2 = my_fftw_allocate_real((ny+1), mm);
	
	{ // Transform all the rows
	fftw_execute(fftw_plan_many_dft_c2r(1, &nx, dummyIi2.rows, dummyIi2.a, &dummyIi2.cols, 1, dummyIi2.cols, dummyRr2.a, &dummyIi2.cols, 1, dummyRr2.cols, FFTW_ESTIMATE));
	}

	double *vp;
	vp = (double*) fftw_malloc((((ny+1)*mm))*sizeof(double)); 
	//use for loop to set dummyR with vp and flip rows and cols
	for (int i = 0; i < ny+1; i++){ 
		for (int j = 0; j < mm; j++){ 
			vp[j*(ny+1)+i] = dummyRr2(i,j);
		}
	}

	double *w;
	w = (double*) fftw_malloc((((ny+1)*mm))*sizeof(double));
	
	//take the product in real space: w=up.*vp;
	for (int i = 0; i < ny+1; i++){ 
		for (int j = 0; j < mm; j++){ 
			w[j*(ny+1)+i] = vp[j*(ny+1)+i] * up[j*(ny+1)+i];
		}
	}

	
	real_buffer inW = my_fftw_allocate_real((ny+1), mm);
	for (int i = 0; i < ny+1; i++){ 
		for (int j = 0; j < mm; j++){ 
			inW(i,j) = w[j*(ny+1)+i];
		}
	}

	// FFT along the x direction:  wh=fft(w,[],2);
	cplx_buffer outW = my_fftw_allocate_cplx((ny+1), mm); //((ny+1), nx);

	 { // Transform all the rows
    fftw_execute(fftw_plan_many_dft_r2c(1, &nx, inW.rows, inW.a, &inW.cols, 1, inW.cols, outW.a, &inW.cols, 1, outW.cols, FFTW_ESTIMATE));
	}

	fftw_complex *whk;
	whk = (fftw_complex*) fftw_malloc((((mm)*(ny+1))*nyk)* sizeof(fftw_complex));  
	memset(whk, 42, (((mm)*(ny+1))*nyk)* sizeof(fftw_complex)); 

	for (int i = 0; i < mm; i++){ 
		for (int j = 0; j < ny+1; j++){ 
			for (int k = 0; k < ncomp; k++){ 
				whk[i*(ny+1)+j][k] = outW(j,i)[k];
			}
		}
	
	}
	Eigen::Map<Eigen::MatrixXcd, Eigen::Unaligned> wh(reinterpret_cast<std::complex<double>*>(whk),(ny+1),mm); 

	ph(all,seqN(0,nx/2)) = wh(all,seqN(0,nx/2));
	ph(all,seq(last+1-nx/2, last)) = wh(all,seq(last+1-nx/2, last));
	ph = 1.5 * ph;
//free stuff
}

static const int nx = 10;
static const int ny = 10; 
static const int mm = nx* 3/2;

struct cplx_buffer
{
  fftw_complex* a;
  int rows;
  int cols;

  fftw_complex& operator()(int i, int j) const { return a[i * cols + j]; }
};

struct real_buffer
{
  double* a;
  int rows;
  int cols;

  double& operator()(int i, int j) const { return a[i * cols + j]; }
};

real_buffer my_fftw_allocate_real(int x, int y) { return real_buffer{ fftw_alloc_real(x * y), x, y }; }
cplx_buffer my_fftw_allocate_cplx(int x, int y) { return cplx_buffer{ fftw_alloc_complex(x * y), x, y }; }

static const int nx = 10;
static const int ny = 10;

using namespace std;
using namespace Eigen;

void aapx(const Eigen::Ref<const Eigen::MatrixXcd>& uh, 
          const Eigen::Ref<const Eigen::MatrixXcd>& vh,
          Eigen::Ref<Eigen::MatrixXcd> ph)

#ifndef spectralFUN_H
#define spectralFUBN_H