python - 带有“quad”和“args”的积分 - 杨氏双缝干涉仪

Question

我必须使用 Python 中的“quad”和“args”积分来解决“Fringes of Young”问题

对于源大小 R，屏幕上 M(X,Y) 的强度公式如下：

源点 S 的坐标(xs=0,ys)为-R/2<=ys<=R/2

我需要创建一个函数来I(X,Y,R)使用“quad”的“args”计算强度。然后，绘制I(0,Y,10e-6)介于 -0.01 和 0.01 之间的 Y，还有I(0,Y,0.002),I(0,Y,0.003),I(0,Y,0.004). 知道我的错在哪里吗？

我的代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

y_min   = -0.01
y_max   = +0.01
R       = y_max-y_min
y       = np.linspace(y_min, y_max, 100)
X       = 0
Y       = 0
d       = 1
D       = 10
s       = 10
Lambda  = 0.5e-3

delta_s   = lambda ys,X,Y : np.sqrt(X**2+(Y-d/2)**2+D**2)+np.sqrt((ys-d/2)**2+s**2)- \
                            np.sqrt(X**2+(Y+d/2)**2+D**2)-np.sqrt((ys+d/2)**2+s**2)
def integrand(y_s,x,y):
    value =  2*(1+np.cos(2*np.pi*delta_s(x,y,y_s)/Lambda))
    return value

def calcul_XYR(X,Y,R):
    compteur  = 0
    I_XYR    = []               # array for I(X,Y,R)
    while compteur < len(y-1):
        Y = y[compteur]
        print(Y)
        I_XYR.append(1/R*quad(integrand, -R/2, R/2, args=(X,Y))[0])
        compteur+=1
    return I_XYR

plt.figure(figsize=(7, 5))
plt.title("Franges de Young - Figure 3")
plt.axis([y_min, 0.015, 0, 4])
plt.xlabel("Y (mm)")
plt.ylabel("Intensity (a.u.)")
plt.plot(y, calcul_XYR(0,Y,1e-6),  '-', color="red",   label=r'$R=10^{-6}$')
plt.plot(y, calcul_XYR(0,Y,0.002), '-', color="blue",  label=r'$R=0.002$')
plt.plot(y, calcul_XYR(0,Y,0.003), '-', color="black", label=r'$R=0.003$')
plt.plot(y, calcul_XYR(0,Y,0.004), '-', color="green", label=r'$R=0.004$')
plt.legend(loc='right', bbox_to_anchor=(1.00, 0.3))
plt.savefig('question 3 figure.pdf', format='pdf')
plt.show()

结果 ：

预期的 ：

我还想绘制（使用带有参数的imshow：cmp（gray），vmin，vmax）对应于I（X，Y，1e-06）的2D图像。（X 在 -10 到 10 之间）。

Answer 1

主要错误是参数的顺序delta_s。它被定义为delta_s = lambda ys, X, Y，但称为delta_s(X, Y, ys)。

另外，calcul_XYR不使用其参数Y，因此最好将其删除。循环可以写成for Y in y。

这是生成所需图的修改后的代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

y_min = -0.01
y_max = +0.01
#R = y_max - y_min
y = np.linspace(y_min, y_max, 100)
X = 0
Y = 0
d = 1
D = 10
s = 10
Lambda = 0.5e-3

def delta_s(X, Y, ys):
    return np.sqrt(X ** 2 + (Y - d / 2) ** 2 + D ** 2) + np.sqrt((ys - d / 2) ** 2 + s ** 2) - \
           np.sqrt(X ** 2 + (Y + d / 2) ** 2 + D ** 2) - np.sqrt((ys + d / 2) ** 2 + s ** 2)

def integrand(y_s, x, y):
    return 2 * (1 + np.cos(2 * np.pi * delta_s(x, y, y_s) / Lambda))

def calcul_XR(X, R):
    I_XYR = []  # array for I(X,Y,R)
    for Y in y:
        I_XYR.append(1 / R * quad(integrand, -R / 2, R / 2, args=(X, Y))[0])
    return I_XYR


plt.figure(figsize=(7, 5))
plt.title("Franges de Young - Figure 3")
plt.axis([y_min, 0.015, 0, 4])
plt.xlabel("Y (mm)")
plt.ylabel("Intensity (a.u.)")
plt.plot(y, calcul_XR(0, 1e-6), '-', color="red", label=r'$R=10^{-6}$')
plt.plot(y, calcul_XR(0, 0.002), '-', color="blue", label=r'$R=0.002$')
plt.plot(y, calcul_XR(0, 0.003), '-', color="black", label=r'$R=0.003$')
plt.plot(y, calcul_XR(0, 0.004), '-', color="green", label=r'$R=0.004$')
plt.legend(loc='right', bbox_to_anchor=(1.00, 0.3))
plt.savefig('question 3 figure.pdf', format='pdf')
plt.show()

以下代码显示了函数的图像：

x_min = -10
x_max = 10
x = np.linspace(x_min, x_max, 100)
R = 1e-6
plt.figure(figsize=(7, 5))
graphe1 = []
for xi in x:
    graphe1.append(calcul_XYR(xi, R))
graphe1 = np.array(graphe1).T  #convert to numpy array and transpose

# imshow normally starts displaying at the top `origin='lower'` reverses this;
# the extent is used to tell imshow what the x and y limits of the image are, to correctly put the ticks
# without `aspect='auto'` imshow seems to want to display x and y with the same scale
# interpolation='bilinear' tells to smooth out the pixels
plt.imshow(graphe1, cmap=plt.cm.gray, vmin=None, vmax=None,
           extent=[x_min, x_max, y_min, y_max],
           aspect='auto', origin='lower', interpolation='bilinear')
plt.xlabel('X')
plt.ylabel('Y')
plt.title(f'R={R}')
plt.show()