“四边形”和“args”杨氏双缝干涉仪积分

Question

我必须用Python中带有“quad”和“args”的积分来解决一个问题“fronges of Young”

对于源尺寸R，屏幕上M（X，Y）的强度公式如下：

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

我需要创建一个函数来使用“quad”的“args”计算强度I(X,Y,R)。 然后，在-0.01和0.01之间绘制Y的I(0,Y,10e-6)，也就是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()

结果:

预期：

我还想绘制（使用参数为cmp（gray）、vmin、vmax的imshow）与I（X、Y、1e-06）对应的二维图像。（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()