在笛卡尔坐标系中计算双重积分，而非(R,Theta)

Question

我之前的问题（Python中强度函数的积分）

你可以在下面的图片中看到衍射模型：

我想计算每个像素（方块）中的强度积分，所以我不能用R和Theta作为变量。我该如何在X-Y坐标中做到这一点呢？

我们的函数是：

我们可以用以下方式代替sin(theta)：

sintheta= (np.sqrt((x)**2 + (y)**2)/(np.sqrt((x)**2 + (y)**2 + d**2)))

其他常数：

lamb=550*10**(-9)
k=2.0*np.pi/lamb
a=5.5*2.54*10**(-2)
d=2.8

当你绘制这个函数时，结果大概是这样的：（上面的图片是从顶部看的视图）

在之前的话题中：计算函数在(0.0, dist)的积分，然后乘以(2*np.pix)，其中x = ka*np.sin(theta)，但现在我想在每个像素中进行积分。之前的方法不适用，因为这是X-Y坐标，而不是极坐标。

Answer 1

其实，在笛卡尔坐标系下进行积分是相对简单的。现在你已经有了强度函数，你需要用坐标 x 和 y 来表示半径 r。这其实在你的问题中已经做过了。

所以，待积分的函数（不包括一些常数）是：

from scipy import special as sp

# Fraunhofer intensity function (spherical aperture)
def f(x,y):
    r = np.sqrt(x**2 + y**2)
    return (sp.j1(r)/r)**2

或者利用这个事实：2 J1(x)/x = J0(x) + J2(x) [谢谢你，Jaime!]：

def f(x,y):
    r = np.sqrt(x**2 + y**2)
    return (sp.j0(r) + sp.jn(2,r))**2

这种形式更好，因为它在任何地方都没有奇点。

现在，我没有使用任何常数因子。如果你想的话可以加上，但我觉得用无限区域的积分结果来归一化更简单。否则很容易忘记某个常数（我通常会忘）。

可以使用 scipy.integrate.nquad 来进行积分。它接受一个多维函数进行积分。所以在这个例子中：

import scipy.integrate

integral = scipy.integrate.nquad(f, ([-d/2, d/2], [-d/2, d/2]))[0]

不过，由于你的函数非常对称，你可以考虑只在一个象限内积分，然后乘以四：

4. * scipy.integrate.nquad(f, ([0, d/2], [0, d/2]))[0]

使用这些，完整的强度是：

>>> 4. * scipy.integrate.nquad(f, [[0,inf],[0,inf]])[0]
12.565472446489999

（顺便说一下，这看起来很像 4π。）当然，你也可以使用极坐标来计算完整的值，因为这个函数具有圆对称性（如在Python中的强度函数积分中所述）。不同的值是由于不同的缩放（在极坐标积分中省略了 2π，这里使用的是贝塞尔函数的和形式，所以有个 2）。

例如，对于一个从 -1 到 1 的正方形区域，归一化后的（除以上完整功率值）功率在正方形区域内是：

>>> 4*scipy.integrate.nquad(f, [[0,1],[0,1]])[0] / 12.565472446489999
0.27011854108867

所以，大约 27% 的入射光照射到这个正方形光电探测器上。

---

关于你的常数，似乎缺少了一些东西（至少是单位）。我猜测：

• 波长：550 纳米
• 圆形光圈直径：0.0055" = 0.14 毫米
• 光圈到传感器的距离：2.8 毫米
• 正方形传感器尺寸 5.4 微米 x 5.4 微米

最后一个我只是根据图像猜测的。由于传感器的尺寸远小于距离，sin(ϴ) 非常接近 y/d，其中 d 是距离，y 是光轴的偏移。使用这些数字 x = ka sin(ϴ) = kay / d ≈ 1.54。对于这个数字，强度积分大约给出 0.52（或 52%）。

如果你在与某个实验值比较，请记住有很多误差来源。成像平面上的图像是光圈的傅里叶变换。如果光圈边缘有小缺陷，它们可能会改变结果的光斑。艾里环通常没有天文学家想象的那么美丽……

只是为了好玩：