菲涅尔传播1（Fresnel Propagation）

菲涅尔积分公式

省去复杂的数学推导，这里我们直接给出菲涅尔积分公式
$$U_2\left(x_2,y_2\right)=\frac{e^{jkz}}{j\lambda z}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{U}_1\left(x_1,y_1\right)e^{j\frac{k}{2z}\left[{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2\right]}d{x}_{1}d{y}_1,\text{\hspace{1em}(1)}$$
这里 $(x_1,y_1)$ 代表光源平面的坐标，$(x_2,y_2)$ 代表观察平面的坐标，$k=2\pi /\lambda$ 表示波数， $\lambda$ 代表光学波长，$z$ 代表光源平面和观察平面中心点的距离。图1展示了光源平面和观察平面的几何构造。

图1. 光源平面和观察平面的几何构造图形示意.
为了方便我们进行模拟菲涅尔传播，我们可以进一步对公式（1）进一步处理。把积分符号里指数部分的平方展开，最终我们可以得到如下的公式，
$$\begin{array}{rl}{U}_2(x_2,y_2)=& \frac{\exp (jkz)}{j\lambda z}\exp \left[j\frac{k}{2z}\left(x_2^2+y_2^2\right)\right]\\ & \times \iint \left\{U_1(x_1,y_1)\exp \left[j\frac{k}{2z}\left(x_1^2+y_1^2\right)\right]\right\}\exp \left[-j\frac{2\pi }{\lambda z}(x_1{x}_2+y_1{y}_2)\right]d{x}_{1}d{y}_1,\end{array}\text{\hspace{1em}(2)}$$
二维傅里叶变换的定义如下
$$G\left(f_{x_1},f_{y_1}\right)=\iint g\left(x_1,y_1\right)\exp \left[-j2\pi \left(f_{x_1}{x}_1+f_{y_1}{y}_1\right)\right]d{x}_{1}d{y}_1,\text{\hspace{1em}(3)}$$
对比公式（2）和公式（3），我们可以发现公式（2）中的积分符号里就是一个傅里叶变换，其中
$$g(x_1,y_1)=U_1(x_1,y_1)\exp \left[j\frac{k}{2z}\left(x_1^2+y_1^2\right)\right],\text{\hspace{1em}(4)}$$
$$f_{x_1}=\frac{x_2}{\lambda z},f_{y_1}=\frac{y_2}{\lambda z}.\text{\hspace{1em}(5)}$$
这样，我们就可以用一个傅里叶变换快速计算积分部分，而快速傅里叶变换（FFT）在大部分软件中都直接可以调用相应的命令。

如何模拟菲涅尔传播

首先，我们需要先构造出一个取样空间，设定好取样间隔。 如图2所示，$x$ 方向上的取样间隔是 $\Delta x$, 取样点数是 $N$, 那么 $x$ 方向上对应的取样物理尺寸就是 $L_X=N\Delta x$。类似的，$y$ 方向上的取样间隔是 $\Delta y$, 取样点数是 $M$, 那么 $y$ 方向上对应的取样物理尺寸就是 $L_Y=M\Delta y$。 通常，我们取 $M=N$, $\Delta x=\Delta y$, 并且设置 $M$ 和 $N$ 为偶数。

图2. 一个二维平面的离散取样。

光源平面的坐标系统

根据图2，我们可以进一步构建光源平面坐标系统，$x_1$ 方向共有 $N$ 个点，从左到右，坐标点数为 $-N/2$ 到 $N/2-1$； $y_1$ 方向共有 $M$ 个点，从下到上，坐标点数为 $-M/2$ 到 $M/2-1$。这里，我用图3进一步阐明这个坐标系统。

图3. 坐标系统。
根据图三所示的坐标系统，我们可以借助软件，例如matlab来实现这一坐标系统。matlab的代码如下

% source plane coordinates
[x1, y1] = meshgrid((-N/2 : 1: N/2 -1) * delta1, (-M/2 : 1: M/2 -1) * delta1);

在上述代码中，我们借助matlab中的命令meshgrid来构造图3所示的一个网格坐标系统。代码中，delta1即是取样间隔，这里我们认为 $\Delta x_1=\Delta y_1={\Delta }_1$。这里 $(x_1,y_1)$ 代表着光源平面的坐标。

观察平面的坐标系统

为了得到观察平面的坐标系统，我们需要知道观察平面的取样间隔。这个取样间隔可以由下面的公式计算得出
$$\Delta x_2=\frac{\lambda z}{N\Delta x_1},\text{\hspace{1em}(6)}$$
在公式（2）中，观察平面的取样间隔可以由传播距离和光源平面的取样间隔得到。因此，我们可以进而得到观察平面的坐标系统，其matlab代码如下

% observation plane coordinates
delta2 = lambda * z / (N * delta1);
[x2, y2] = meshgrid((-N/2 : 1: N/2 -1) * delta2, (-M/2 : 1: M/2 -1) * delta2);

根据公式（2）模拟菲涅尔传播

根据公式（2），我们可以定义如下菲涅尔传播函数。

function [U2, x2, y2] = direct_fresnel_propagation(U1, lambda, delta1,z)
% U1: source plane
% lambda: wavelength
% delta1: sampling interval in the source plane, delta_x1 and delta_y1,
% delta_x1 = delta_y1
% z:distance of the central point between source plane and observation
% plane, i.e., the propagation distance
% U2: output in the observation plane
% x2, y2:the coordinates in the observation plane
[M, N] = size(U1);  % obtain the number of points along each side
k = 2 * pi / lambda; % wave number
% source plane coordinates
[x1, y1] = meshgrid((-N/2 : 1: N/2 -1) * delta1, (-M/2 : 1: M/2 -1) * delta1);
% observation plane coordinates
delta2 = lambda * z / (N * delta1);
[x2, y2] = meshgrid((-N/2 : 1: N/2 -1) * delta2, (-M/2 : 1: M/2 -1) * delta2);
% Fresnel propagation integral
U2 = (exp^(1j * k * z)/1j * lambda * z) * exp(1j * (k/(2*z)) * (x2.^2 + y2.^2))...
    .* fftshift(fft2(fftshift(U1 .* exp(1j * (k/(2*z)) * (x1.^2 + y1.^2))))) * delta1^2;
end

在上述定义的函数中，有几点值得注意：

1. 虚数单位 $j$ 加个数字1是为了增加程序稳定性, 不加这个数字1也是可以的，前提是不能在其它地方再给 $j$ 赋值了，否则引发的错误你很难发现；
2. 程序中点乘的地方需要特别注意，例如
   这里的点乘代表对向量或者矩阵中每一个元素进行操作；
3. 程序中，…代表换行续写；
4. 最后要乘以delta1^2这个比例因子来对应真实的尺寸。
5. 程序中要用fftshift(fft2(fftshift( )))来实现一个标准的数学上的傅里叶变换，用fftshift的根本原因是软件的索引只能从正整数开始。但是我们不需要关心这一点，可以简单的认为用fftshift(fft2(fftshift( )))就是来实现二维傅里叶变换。

模拟一个圆孔和方形孔的菲涅尔衍射

圆孔和方形孔的定义

圆孔定义如下

$$\mathrm{circ}\left(\frac{\sqrt{x^2+y^2}}{a}\right)=\{\begin{array}{ll}1& \sqrt{x^2+y^2}<a\\ \frac{1}{2}& \sqrt{x^2+y^2}=a\\ 0& \sqrt{x^2+y^2}>a\end{array}\text{\hspace{1em}(7)}$$
这里 $a$ 代表圆孔的半径。
根据定义，我们可以用如下代码实现这个函数：

function [circle_shape] = circ(x,y,a)
% x,y: coordinates
% a:radius used to control the size of circle
r = sqrt(x.^2 + y.^2);
circle_shape = double(r<a);
circle_shape(r==a) = 0.5;
end

方形孔定义如下

一维的方形函数定义如下
$$\mathrm{rect}\left(\frac{x}{a}\right)=\{\begin{array}{cc}1& |x|<\frac{a}{2}\\ \frac{1}{2}& |x|=\frac{a}{2}\\ 0& |x|>\frac{a}{2}\end{array}$$
其对应的实现代码如下

function [out] = rect(x,a)
% x:coordinates
% a:used to control the width of rectangular
out = double(abs(x)<a/2);
out(abs(x)==a/2) = 0.5;
end

二维的方形孔，可以用两个一维的方形函数构造

% two dimensional rectangular shape
rect_shape = rect(x, a) .* rect(y, b);

圆孔的菲涅尔衍射

圆孔的菲涅尔衍射代码如下

clear;close all;clc
%% parameters
a = 1e-3;   % radius of circle 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
%% circle shape
% cooridates in the source plane
[x1, y1] = meshgrid((-N/2 : 1: N/2 -1) * delta1, (-N/2 : 1: N/2 -1) * delta1);
circle_shape = circ(x1, y1, a);
figure(1) 
imagesc(x1(1,:),y1(:,1),circle_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_circ,x2, y2] = direct_fresnel_propagation(circle_shape,lambda,delta1,z);
figure(2),
imagesc(x2(1,:),y2(:,1),abs(Fres_circ)); %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);  

产生的圆形孔和相应的衍射光斑如下图所示

（a）

(b)
图4. （a）圆形孔，（b）圆形孔对应的衍射光斑。

方形孔的菲涅尔衍射

方形孔的菲涅尔衍射代码如下

clear;close all;clc
%% parameters
a = 1e-3;   % radius of circle 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);  

产生的方形孔和相应的衍射光斑如下图所示

(a)

(b)
图5. （a）方形孔，（b）方形孔对应的衍射光斑。