我正在尝试制作一个圆形衍射图案,它的中心点被一系列环包围。它涉及到代码中定义的贝塞尔积分。
我的问题是它需要很长时间,就像我等待 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()