# 2D Fourier Synthesis (IDFT) not yielding expected result

I am trying to recover a 2D signal using inverse DFT, to my understanding the IDFT outputs the coefficients of the fourier series of the original function up to the Nyquist frequency.

So for example I'm trying with the function $$f(x) = \exp(-x^Tx),\quad x \in\mathbb{R}^2$$ $$\hat{f}(z) = \mathcal{F}[f(x)](z) = \frac{1}{2} \cdot \exp \left(-\frac{1}{4} z^Tz\right), \quad z \in\mathbb{R}^2$$

Then the 2D FS should be (and I'm not sure here) $$s_f(x, y) = \sum_{j=-N}^N \sum_{k=-N}^N c_{j,k} \cdot \exp\left(\frac{2 \pi ijx}{l_x}\right) \cdot \exp\left(\frac{2 \pi iky}{l_y}\right)$$

In the following code I start from $$\hat{f}$$ and try to get to $$f$$

import numpy as np
from numpy.fft import ifft2, fftshift

size = 51
shape = (size, size)
bound = np.pi

base = np.linspace(-bound, bound, size)
domx, domy = np.meshgrid(base, base)
domain = np.array([domx, domy])

# the einsum is to perform the dot product over the whole domain
transform = 0.5 * np.exp(-0.25 * np.einsum('jki,jki->jk', domain.T, domain.T))

n = np.arange(size)-size//2
nx, ny = np.meshgrid(n, n)

coeffs = fftshift(ifft2(transform))
fs = np.vectorize(
lambda x, y: np.sum(coeffs *
np.exp(2 * np.pi * 1j * nx * x / (2 * bound)) *
np.exp(2 * np.pi * 1j * ny * y / (2 * bound))
))

reconstruction = fftshift(np.real(fs(domx, domy)))
expected = np.exp(-np.einsum('jki,jki->jk', domain.T, domain.T))


The thing is that the reconstruction is not as close to the expected result as I thought it should be, I'm not sure if it is the expected behavior of the approximation or if there is some conceptual error somewhere.

Here is the result I got, expected on the left and reconstruction on the right

My idea is that if it was an approximation error and I increase the resolution, the difference should decrease, but here is the same with size 151 (previous was with size 51)

So I think it is a conceptual issue that I am not getting right, I'd like to understand what is happening here.