### 夫琅和费衍射图样的三维强度分布 夫琅和费衍射描述了当光通过一个小孔或障碍物时，在远处屏幕上形成的干涉图案。对于不同形状的孔径，其对应的衍射花样会有所不同。 为了实现这一目标，可以采用数值计算的方法来模拟这些现象，并利用 Python 中强大的科学计算库如 NumPy 和 Matplotlib 来完成可视化工作。下面是一个具体的例子，展示了如何创建一个函数 `fraunhofer_diffraction` ，该函数接受参数定义孔径类型以及尺寸大小并返回相应的三维强度分布图像： #### 圆形孔径 圆形孔径下的远场（夫琅禾费）衍射可以通过 Airy 函数表示[^1]。这里提供了一个简单的 Python 实现方法用于生成圆形单缝衍射的效果: ```python import numpy as np from scipy.special import jv import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D def airy_intensity(x, y, D=0.001): # 默认直径为1mm k = 2 * np.pi / (589e-9) # 波长设为钠黄光波长589nm r = np.sqrt((x ** 2 + y ** 2)) I_airy = (jv(1, k*D*r/2)/(k*D*r/2))**2 return I_airy # 创建网格数据集 X = Y = np.linspace(-0.005, 0.005, 500) X, Y = np.meshgrid(X, Y) Z_circular = airy_intensity(X, Y) fig = plt.figure(figsize=(10,7)) ax = fig.add_subplot(projection='3d') surf = ax.plot_surface(X*1e6, Y*1e6, Z_circular, cmap='viridis') plt.title('Circular Aperture Fraunhofer Diffraction Pattern Intensity Distribution', pad=20) ax.set_xlabel(r'$\mu m$') ax.set_ylabel(r'$\mu m$') ax.set_zlabel('Intensity') plt.show() ``` #### 方形孔径 方形孔径的情况相对复杂一些，因为这涉及到两个方向上的独立傅里叶变换过程。但是基本思路仍然是相似的——我们只需要改变输入到 FFT 的二维数组即可得到正确的结果。 ```python def square_aperture(u, v, a=0.001): sinc_u = np.sinc(a*u/(np.pi*(a))) sinc_v = np.sinc(a*v/(np.pi*(a))) intensity_square = abs(sinc_u*sinc_v)**2 return intensity_square U = V = np.fft.fftfreq(len(X), d=X[0][1]-X[0][0]) UU,VV=np.meshgrid(U,V) Z_square = square_aperture(UU, VV) fig = plt.figure(figsize=(10,7)) ax = fig.add_subplot(projection='3d') surf = ax.plot_surface(UU*1e6, VV*1e6, Z_square, cmap='plasma') plt.title('Square Aperture Fraunhofer Diffraction Pattern Intensity Distribution', pad=20) ax.set_xlabel(r'Frequency ($\mu$m)') ax.set_ylabel(r'Frequency ($\mu$m)') ax.set_zlabel('Normalized Intensity') plt.show() ``` #### 正六边形孔径 正多边形孔径则更加具有挑战性；然而，如果考虑的是规则图形，则仍然能够找到解析解。对于正六边形而言，可以将其分解成六个相同的三角区域来进行处理。 ```python def hexagon_aperture(x, y, side_length=0.001): def inside_hexagon(p_x, p_y): dist_to_center = lambda px, py: ((px - center[0])**2+(py-center[1])**2)**0.5 angles = [] for i in range(6): angle_between_points = np.arctan2(vertices[i%6][1], vertices[i%6][0]) \ - np.arctan2(vertices[(i+1)%6][1],vertices[(i+1)%6][0]) while(angle_between_points<0):angle_between_points+=2*np.pi angles.append(angle_between_points) total_angle_sum = sum([abs(ang) for ang in angles]) point_vector = [p_x-center[0], p_y-center[1]] vector_magnitude = dist_to_center(*point_vector) if all(abs(np.dot(point_vector,[np.cos(theta), np.sin(theta)]) /vector_magnitude)<=side_length/2 for theta in angles)\ and total_angle_sum<=2*np.pi: return True else: return False centers = [[0]*2]*6 sides = [-side_length/2]+list(range(side_length//2,(int)(side_length)+1)) global_vertices = [(sides[i % len(sides)], sides[(i+1)//len(sides)]) for i in range(len(sides)*2)] mask = np.zeros_like(x).astype(bool) height,width=x.shape for row in range(height): for col in range(width): mask[row,col]=inside_hexagon(col,row) return mask.astype(float) center=[0,0] vertices=[[side_length*np.cos(i*np.pi/3), side_length*np.sin(i*np.pi/3)]for i in range(6)] mask=hexagon_aperture(X,Y) I_hexagonal_fftshifted = np.abs(np.fft.fftshift(np.fft.fft2(mask)))**2 extent = (-width/2*width/wavelength*focal_length, width/2*width/wavelength*focal_length, -height/2*height/wavelength*focal_length, height/2*height/wavelength*focal_length) fig=plt.figure(figsize=(10,7)) ax=fig.gca(projection="3d") surf=ax.imshow(I_hexagonal_fftshifted,cmap='hot', extent=extent,vmin=I_hexagonal_fftshifted.min(), vmax=I_hexagonal_fftshifted.max()) plt.colorbar(surf,label="Intensity",fraction=.04,pad=-0.02) plt.xlabel("Distance from Center (m)") plt.ylabel("Distance from Center (m)") plt.title("Hexagonal Aperture Fraunhofer Diffraction Pattern\n(Intensity Distribution)",pad=20) plt.tight_layout() plt.show() ```