Return to Answer

It should be added that the function describing the impulse response ($11$) is called the Jinc function (in analogy to the sinc function):

$$\mathrm{jinc}(x)=\frac{J_1(x)}{x}\tag{12}$$

Using ($12$), the impulse response ($11$) can be written as

$$h[n_1,n_2]=\frac{\omega_c^2}{2\pi}\,\textrm{jinc}(\omega_c\rho),\qquad \rho=\sqrt{n_1^2+n_2^2}\tag{13}$$

from scipy import special
import numpy as np

def circularLowpassKernel(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size.
  kernel = np.fromfunction(lambda x, y: omega_c*special.j1(omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  if N % 2:
    kernel[(N - 1)//2, (N - 1)//2] = omega_c**2/(4*np.pi)
  return kernel

from scipy import special
import numpy as np

def circularLowpassKernel(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size.
  with np.errstate(divide='ignore',invalid='ignore'):
    kernel = np.fromfunction(lambda x, y: omega_c*special.j1(omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  if N % 2:
    kernel[(N - 1)//2, (N - 1)//2] = omega_c**2/(4*np.pi)
  return kernel

from scipy import special
import numpy as np

def circularLowpassKernel(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda x, y: omega_c*special.j1(omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = omega_c**2/(4*np.pi)
  return kernel

import matplotlib.pyplot as plt

kernelN = 11  # Horizontal size of the kernel, also its vertical size. Must be odd.
omega_c = np.pi  # Cutoff frequency in radians <= pi
kernel = circularLowpassKernel(omega_c, kernelN)
plt.imshow(kernel, vmin=-1, vmax=1, cmap='bwr')
plt.colorbar()
plt.show()

kernelN = 41  # Horizontal size of the kernel, also its vertical size. Must be odd.
omega_c = np.pi/4  # Cutoff frequency in radians <= pi
kernel = circularLowpassKernel(omega_c, kernelN)
plt.imshow(kernel, vmin=-np.max(kernel), vmax=np.max(kernel), cmap='bwr')
plt.colorbar()
plt.show()

from scipy import special
import numpy as np

def circularLowpassKernel(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size.
  kernel = np.fromfunction(lambda x, y: omega_c*special.j1(omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  if N % 2:
    kernel[(N - 1)//2, (N - 1)//2] = omega_c**2/(4*np.pi)
  return kernel

import matplotlib.pyplot as plt

kernelN = 11  # Horizontal size of the kernel, also its vertical size.
omega_c = np.pi  # Cutoff frequency in radians <= pi
kernel = circularLowpassKernel(omega_c, kernelN)
plt.imshow(kernel, vmin=-1, vmax=1, cmap='bwr')
plt.colorbar()
plt.show()

kernelN = 41  # Horizontal size of the kernel, also its vertical size.
omega_c = np.pi/4  # Cutoff frequency in radians <= pi
kernel = circularLowpassKernel(omega_c, kernelN)
plt.imshow(kernel, vmin=-np.max(kernel), vmax=np.max(kernel), cmap='bwr')
plt.colorbar()
plt.show()