菲涅尔传播2（Fresnel Propagation）

前言

在上一篇文章中，我们介绍了如何模拟菲涅尔传播，并用两个形状（圆形和方形）来观察相应的菲涅尔衍射光斑。但是，有一个问题我们还需要解决，那就是如何知道我们模拟的正确不正确呢？在这篇文章中，我们进一步研究如何验证我们写的程序正确不正确。

方形孔的菲涅尔衍射分析解

数学上很难求得不同形状的菲涅尔衍射分析解，目前具有较好分析解的是对于方形孔的菲涅尔衍射，这里我们给出它的具体数学表达式

U 2 ( x 2 , y 2 ) = e j k z 2 j { [ C ( α 2 ) − C ( α 1 ) ] + j [ S ( α 2 ) − S ( α 1 ) ] } × { [ C ( β 2 ) − C ( β 1 ) ] + j [ S ( β 2 ) − S ( β 1 ) ] } , (1) \begin{aligned} U_{2}(x_2, y_2)= & \frac{e^{j k z}}{2 j}\left\{\left[C\left(\alpha_{2}\right)-C\left(\alpha_{1}\right)\right]+j\left[S\left(\alpha_{2}\right)-S\left(\alpha_{1}\right)\right]\right\} \\ & \times\left\{\left[C\left(\beta_{2}\right)-C\left(\beta_{1}\right)\right]+j\left[S\left(\beta_{2}\right)-S\left(\beta_{1}\right)\right]\right\} \end{aligned},\tag{1} U2​(x2​,y2​)=​2jejkz​{[C(α2​)−C(α1​)]+j[S(α2​)−S(α1​)]}×{[C(β2​)−C(β1​)]+j[S(β2​)−S(β1​)]}​,(1)
其中 $C$ 和 $S$ 被称为菲涅尔积分，其具体表达式为
$$C\left(\chi \right)={\int }_0^{\chi }\cos \left(\frac{\pi }{2}{t}^2\right)dt,S\left(\chi \right)={\int }_0^{\chi }\sin \left(\frac{\pi }{2}{t}^2\right)dt,\text{\hspace{1em}(2)}$$
公式（1）中的 ${\alpha }_1$, ${\alpha }_2$, ${\beta }_1$, 和 ${\beta }_2$ 的表达式为
$${\alpha }_1=-\sqrt{\frac{2}{\lambda z}}\left(\frac{a}{2}+x_2\right),{\alpha }_2=\sqrt{\frac{2}{\lambda z}}\left(\frac{a}{2}-x_2\right),\text{\hspace{1em}(3)}$$
$${\beta }_1=-\sqrt{\frac{2}{\lambda z}}\left(\frac{a}{2}+y_2\right),{\beta }_2=\sqrt{\frac{2}{\lambda z}}\left(\frac{a}{2}-y_2\right),\text{\hspace{1em}(4)}$$

方形孔的菲涅尔衍射分析解代码实现

function [U2] = analytical_square_fresnel_propagation(x2,y2,a,lambda,z)
% x2, y2: coordinates in the observation plane
% a: width of the square
% lambda: wavelength
% z:propagation distance
alpha1 = -sqrt(2/(lambda*z)) * (a/2 + x2);
alpha2 = sqrt(2/(lambda*z)) * (a/2 - x2);
beta1 =  -sqrt(2/(lambda*z)) * (a/2 + y2);
beta2 =  sqrt(2/(lambda*z)) * (a/2 - y2);
% Fresnel sine and cosine integrals
C_alpha1 = fresnelc(alpha1);
C_alpha2 = fresnelc(alpha2);
C_beta1 = fresnelc(beta1);
C_beta2 = fresnelc(beta2);
S_alpha1 = fresnels(alpha1);
S_alpha2 = fresnels(alpha2);
S_beta1 = fresnels(beta1);
S_beta2 = fresnels(beta2);
% observation plane field
k = 2 * pi / lambda; % wave number
U2 = exp(1j * k * z) / (2 * 1j) * ((C_alpha2 - C_alpha1) + 1j * (S_alpha2 - S_alpha1))...
    .* ((C_beta2 - C_beta1) + 1j * (S_beta2 - S_beta1));
end

在上述代码中，fresnelc 和 fresnels 是matlab中自带的计算菲涅尔积分（2）的函数命令，可以直接调用。上述代码是完全对应公式（2）到（4）的，不难理解。

方形孔的菲涅尔衍射分析解与模拟解的比较

具体代码如下所示

clear;close all;clc
%% parameters
a = 1e-3;    % width of rectangular shape unit:m
N = 1024;    % number of points along each side (row and column)
L = 1e-2;    % physical size of the source plane, unit: m
delta1 = L / N; % sampling interval in the source plane
lambda = 532e-9; % wavelength,unti: m
z = 0.1;        % propagation distance, unit:m
%% rectangular shape
% cooridates in the source plane
[x1, y1] = meshgrid((-N/2 : 1: N/2 -1) * delta1, (-N/2 : 1: N/2 -1) * delta1);
rect_shape = rect(x1, a) .* rect(y1, a);
figure(1) 
imagesc(x1(1,:),y1(:,1),rect_shape); %image display 
colormap('gray'); %linear gray display map 
axis square;      % square figure 
axis xy           % y positive up 
xlim([-0.005,0.005])  % x axis range
ylim([-0.005,0.005])  % y axis range
set(gca,'FontSize',16) % Fontsize of x, y ticks
% Fontsize of x,y labels
xlabel('x [m]', 'FontSize', 16); ylabel('y [m]', 'FontSize', 16);  
%% Fresnel propagation of circle shape
[Fres_rect,x2, y2] = direct_fresnel_propagation(rect_shape,lambda,delta1,z);
%%
figure(2),
imagesc(x2(1,:),y2(:,1),abs(Fres_rect)); %image display 
colormap('gray'); %linear gray display map 
axis square;      % square figure 
axis xy           % y positive up 
set(gca,'FontSize',16) % Fontsize of x, y ticks
% Fontsize of x,y labels
xlabel('x [m]', 'FontSize', 16); ylabel('y [m]', 'FontSize', 16);  
%%
U2 = analytical_square_fresnel_propagation(x2,y2,a,lambda,z);
%%
figure(3),
imagesc(x2(1,:),y2(:,1),abs(U2)); %image display 
colormap('gray'); %linear gray display map 
axis square;      % square figure 
axis xy           % y positive up 
set(gca,'FontSize',16) % Fontsize of x, y ticks
% Fontsize of x,y labels
xlabel('x [m]', 'FontSize', 16); ylabel('y [m]', 'FontSize', 16); 
%% compare the 1D result
% compare the simulation result with the analytical result based the 513 row 
sim = abs(Fres_rect(513,:));
ana = abs(U2(513,:));
figure(4),
plot(x2(1,:),sim,'--r'),hold on
plot(x2(1,:),ana,'b')
xlim([min(x2(1,:)),max(x2(1,:))]);ylim([0,1.5]);
set(gca,'FontSize',16) % Fontsize of x, y ticks
% Fontsize of x,y labels
xlabel('x [m]', 'FontSize', 16); ylabel('Magnitude', 'FontSize', 16); 
legend('Simulation Result','Analytical Result')

上述代码中，分析解的计算部分耗时较长，需要耐心等待。运行结果如下图所示

(a) 光源平面的方形孔；
(b) 模拟菲涅尔传播的结果，基于上一篇文章中的公式（2）；

© 菲涅尔传播的分析解；

(d) 菲涅尔传播的模拟结果和分析结果的比较，二者的比较是基于（b）和（c）的第513行进行比较的。
图1.(a) 光源平面的方形孔；(b) 模拟菲涅尔传播的结果，基于上一篇文章中的公式（2）；© 菲涅尔传播的分析解；(d) 菲涅尔传播的模拟结果和分析结果的比较，二者的比较是基于（b）和（c）的第513行进行比较的。

结论

从图1中的结果可以看到，我们模拟的计算结果和分析解的结果，吻合度较高，进而证明我们的模拟结果是有效的。