python - ベッセル積分を使用して回折パターンを作成するPythonを使用した密度プロットですが、実行が停止しません

Question

一連のリングに囲まれた中心スポットを持つ円形の回折パターンを作成しようとしています。それを行うには、コードで定義されているベッセル積分が必要です。

私の問題は、コードが実行されるのを 10 分待ったのに何も表示されなかったように、時間がかかりすぎることです。私のベッセル積分にはポイントごとに1000回の反復があるため、これを手伝ってくれる人はいますか?

私は正しい軌道に乗っていますか？

Mark Newmans の本 Computational Physics で Python と計算物理学を独学しようとしています。演習は計算物理学の 5.4 です。この章へのリンクは次のとおりです。9ページにあります。 http://www-personal.umich.edu/~mejn/cp/chapters/int.pdf

これが私が作ろうとしているイメージです。

.

私のコード:

import numpy as np
import pylab as plt
import math as mathy

#N = number of slicices 
#h = b-a/N 

def J(m,x): #Bessel Integral
    def f(theta):
        return (1/mathy.pi)*mathy.cos(m*theta - x*mathy.sin(theta)) #I replaced np. with mathy. for this line

    N = 1000
    A = 0
    a=0
    b=mathy.pi
    h = ((b-a)/N)

    for k in range(1,N,2):

        A +=  4*f(a + h*k)

    for k in range(2,N,2):

        A +=  2*f(a + h*k)

    Idx =  (h/3)*(A + f(a)+f(b))

    return Idx

def I(lmda,r): #Intensity
    k = (mathy.pi/lmda)    
    return ((J(1,k*r))/(k*r))**2

wavelength = .5        # microm meters
I0 = 1
points = 500           
sepration = 0.2  

Intensity = np.empty([points,points],np.float)

for i in range(points):
    y = sepration*i
    for j in range(points):
        x = sepration*j
        r = np.sqrt((x)**2+(y)**2)

        if r < 0.000000000001:
            Intensity[i,j]= 0.5 #this is the lim as r  -> 0, I -> 0.5
        else: 
            Intensity[i,j] = I0*I(wavelength,r)

plt.imshow(Intensity,vmax=0.01,cmap='hot')
plt.gray()
plt.show()