夫琅禾费衍射MATLAB代码：小孔、矩孔、双缝

绘图结果

代码

%夫琅禾费衍射例子
%%1.圆孔
clc;             % 清除命令行窗口
clear all;       % 清除工作区中的所有变量
close all;       % 关闭所有打开的图形窗口

% 参数设置
L = 2e-1;         % 计算区域的物理尺寸 (m)
N = 2^10;         % 网格点的数量
dx = L / N;      % 网格的空间分辨率
x = linspace(-L/2, L/2, N); % x方向的坐标
y = x;                      % y方向的坐标，与x相同
[X, Y] = meshgrid(x, y);    % 创建二维网格坐标，用于计算

% 圆孔参数和波长设置
w = 1e-3;            % 圆孔的半径 (m)
lambda = 632.8e-9;   % 光波的波长 (m)，红光
z = 50;              % 圆孔到屏幕的距离 (m)
k = 2 * pi / lambda; % 波数，k = 2π/λ
lz = lambda * z;     % λz，用于后续计算
XY = X.^2 + Y.^2; % 网格上每一点到圆心的距离的平方
R = sqrt(X.^2 + Y.^2); % 到圆心的距离

% 圆孔场分布和夫琅禾费衍射公式
t1 = double(R <= w);                % 圆孔
T1 = w^2 .* (jinc(w / lz * R));     % 圆孔初始场分布

% 夫琅禾费衍射公式的实现
H = exp(1i * k * z) / (1i * lambda * z) .* exp(1i * k * (XY) / (2 * z)); % 衍射传递函数
U1 = T1 .* H;       % 圆孔后场的远场表示

% 计算光强分布
% I2 = (U1 ./ lz).^2; % 计算光场的强度分布
I1 = abs(U1).^2;

% 绘制圆孔结果
figure('Color', 'w', 'Name', '圆孔');
mm = 1;
nn = 3;

subplot(mm, nn, 1);
imagesc(x, y, t1);
colormap('gray');
xlabel('x (m)'); ylabel('y (m)');
axis square; axis xy; title('圆孔');

subplot(mm, nn, 2);
imagesc(x, y, nthroot(I1, 3));
colormap('gray');
xlabel('x (m)'); ylabel('y (m)');
axis square; axis xy; title('通过小孔的夫琅禾费衍射');

subplot(mm, nn, 3);
plot(x, I1(N/2, :));
xlabel('x (m)'); ylabel('光强度 (W/m^2)');
title('x轴轮廓'); axis square;

%% 2.矩孔
a = 5e-3; %a、b有一个值非常小或去除，可以为单缝衍射
b = 2e-3;
% t2 = rect(X / b);         %单缝
% T2 = a * sinc(a * X / lz);
t2 = rect(X / a) .* rect(Y / b);         % 矩孔，a为x方向长度，b为y方向长度
T2 = a * b * sinc(a * X / lz) .* sinc(b * Y / lz); % 矩孔初始场分布，对t2傅里叶变换
U2 = T2 .* H;
I2 = abs(U2).^2;

% 绘制矩孔结果
figure('Color', 'w', 'Name', '矩孔');

subplot(mm, nn, 1);
imagesc(x, y, t2);
colormap('gray');
xlabel('x (m)'); ylabel('y (m)');
axis square; axis xy; title('矩孔');

subplot(mm, nn, 2);
imagesc(x, y, nthroot(I2, 3)); % 用三次方根增强对比度
colormap('gray');
xlabel('x (m)'); ylabel('y (m)');
axis square; axis xy; title('通过矩孔的夫琅禾费衍射');

subplot(mm, nn, 3);
plot(x, I2(N/2, :));
xlabel('x (m)'); ylabel('光强度 (W/m^2)');
title('x轴轮廓'); axis square;

%% 3.双缝衍射
a = 2e-3; 
d = 6e-2;
t3 = rect((X - d/2) / a) + rect((X + d/2) / a);         % 矩孔，a为x方向长度，b为y方向长度
T3 = 2 * a * sinc(a * X / lz) .* cos(pi * X / lz * d); % 矩孔初始场分布，对t2傅里叶变换
U3 = T3 .* H;
I3 = abs(U3).^2;

% 绘制矩孔结果
figure('Color', 'w', 'Name', '双缝');

subplot(mm, nn, 1);
imagesc(x, y, t3);
colormap('gray');
xlabel('x (m)'); ylabel('y (m)');
axis square; axis xy; title('双缝');

subplot(mm, nn, 2);
imagesc(x, y, nthroot(I3, 3)); % 用三次方根增强对比度
colormap('gray');
xlabel('x (m)'); ylabel('y (m)');
axis square; axis xy; title('通过双缝的夫琅禾费衍射');

subplot(mm, nn, 3);
plot(x, I3(N/2, :));
xlabel('x (m)'); ylabel('光强度 (W/m^2)');
title('x轴轮廓'); axis square;

%% 定义 jinc 和 rect 函数
% 定义 jinc 函数：J1(2πx) / x
function [out] = jinc(x)
    mask = (x ~= 0);          % 找到 x 中所有不等于 0 的元素
    out = pi * ones(size(x)); % 初始化输出为 π
    out(mask) = besselj(1, 2 * pi * x(mask)) ./ x(mask); % 计算 jinc 函数值
end

% 定义矩形函数 rect(x)
function y = rect(x)
    y = zeros(size(x));
    y(abs(x) <= 0.5) = 1; % |x| <= 0.5 时，函数值为 1
end

原理（傅里叶光学第三版-吕乃光 p94  3.5.2一些简单孔径的夫琅禾费衍射）