Bessel function using scipy.quad integration not working

178 Views Asked by At

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()
0

There are 0 best solutions below