Gridline cutouts after FFT and iFFT on Python

Question

EDIT: I think I messed up on the coordinates of $(p,q)$. Num was missing a multiple of $2\pi/N$. Assuming my interpretation of DFT isn't wrong.

I am currently using FFT to run Fresnel Diffraction as given by:

However, due to the fact that using Fourier transform would put the equation into "Frequency space," I would need to use iFFT to put it back to the original $(x,y)$ coordinates. However, the plot of $E^2$ I get after the FFT and iFFT has grid lines that scale with my z value (the thickness and spacing of the lines scale with z). For example:

1. Original plot without FFTs

2. z=1 plot

3. z=20 plot

I can understand that DFT would mean information loss, but why does it scale with my variable z? Also is there a way to mitigate this? Or am I doing something wrong? Code below:

def Fresnel1(E, z, x_space, y_space):
    ii = len(x_space)
    kk = len(y_space)
    num = np.array(range(ii))
    U1 = np.zeros((ii, kk), dtype=complex)
    for i in range(ii):
        for k in range(kk):
            U1[i,k] = E[i,k]*np.exp(1j*np.pi*(x_space[i]**2+y_space[k]**2)/(lamb*z))
    F1 = fft.fft2(U1)
    D1 = np.zeros((ii, kk), dtype=complex)
    for i in range(ii):
        for k in range(kk):
            D1[i,k] = np.exp(1j*konst*z)/(1j*lamb*z)*np.exp(1j*np.pi*lamb*z*(num[i]**2+num[k]**2))*F1[i,k]
    F2 = fft.ifft2(D1)
    return F2

def find_int(E):
    I1 = np.zeros((len(X1),len(Y1)))
    for i in range(len(X1)):
        for k in range(len(Y1)):
            I1[i,k] = np.real(E[i,k] * np.conj(E[i,k]))
    return np.round(I1, decimals=16)

E2 = Fresnel1(E0, 20, X1, Y1)
intensity_a = find_int(E2)
intensity_a1 = intensity_a/np.max(intensity_a)
fig, ax = plt.subplots(constrained_layout=True)
cs = ax.contourf(X1, Y1, intensity_a1, 100)
cbar = fig.colorbar(cs)