I'm trying to plot a Bessel function in 2D using the scipy.quad integration method.
My problem: the plot only has a tiny dot in the middle because my lowest value is far too low.
Here is my code:
import numpy as np
import scipy.integrate as si
import matplotlib.pyplot as plt
wavelength = 5.0
k = 2 * np.pi / wavelength
side = 100.0
points = 100
spacing = side / points
def J(x, m):
def f(t):
return np.cos(m * t - x * np.sin(t)) / np.pi
y = si.quad(f, 0, np.pi)
return y[0]
x1 = side / 2
y1 = side / 2
data = np.empty([points, points], float)
data[0, 0] = 0.5
for i in range(1, points):
y = spacing * i
for j in range(1, points):
x = spacing * j
r = np.sqrt((x - x1) ** 2 + (y - y1) ** 2)
if r == 0:
data[i, j] = 1
else:
data[i, j] = (2 * J(1, r) / r) ** 2
print(data.min(), data.max())
plt.imshow(data, origin="lower", extent=[0, side, 0, side])
plt.show()