当前位置：首页 > news >正文

MATLAB双缝干涉实验模拟程序

news 2025/10/13 10:45:02

MATLAB双缝干涉实验模拟程序

% 双缝干涉实验模拟
clear; clc; close all;%% 参数设置
lambda = 632.8e-9;         % 波长，氦氖激光 (m)
d = 0.5e-3;                % 双缝间距 (m)
a = 0.1e-3;                % 缝宽 (m)
D = 1;                     % 屏幕距离 (m)
L = 0.01;                  % 观察屏幕宽度 (m)% 计算参数
k = 2 * pi / lambda;       % 波数%% 方法1: 波动光学方法（连续波）
fprintf('=== 双缝干涉模拟 ===\n');
fprintf('波长: %.1f nm\n', lambda*1e9);
fprintf('双缝间距: %.2f mm\n', d*1e3);
fprintf('缝宽: %.2f mm\n', a*1e3);
fprintf('屏幕距离: %.1f m\n', D);% 观察屏幕坐标
N_points = 1000;
x = linspace(-L/2, L/2, N_points);
y = 0;  % 在y=0平面上观察% 计算干涉图样
I_total = zeros(size(x));for i = 1:length(x)% 从两个缝到观察点的距离r1 = sqrt(D^2 + (x(i) - d/2)^2);r2 = sqrt(D^2 + (x(i) + d/2)^2);% 考虑缝宽的效应（多个点源叠加）N_slit_points = 100;slit_points = linspace(-a/2, a/2, N_slit_points);E1 = 0; E2 = 0;for j = 1:N_slit_points% 缝1的贡献r1j = sqrt(D^2 + (x(i) - d/2 - slit_points(j))^2);E1 = E1 + exp(1i * k * r1j);% 缝2的贡献r2j = sqrt(D^2 + (x(i) + d/2 - slit_points(j))^2);E2 = E2 + exp(1i * k * r2j);end% 总电场和光强E_total = E1 + E2;I_total(i) = abs(E_total)^2;
end% 归一化光强
I_total = I_total / max(I_total);%% 方法2: 光子计数方法（随机过程）
fprintf('正在进行光子计数模拟...\n');N_photons = 10000;         % 光子数量
photon_positions = zeros(1, N_photons);% 概率分布函数（基于波动光学结果）
x_fine = linspace(-L/2, L/2, 10000);
I_fine = interp1(x, I_total, x_fine, 'spline');
I_fine = I_fine / sum(I_fine);  % 归一化为概率分布% 生成随机光子位置
cum_prob = cumsum(I_fine);
for i = 1:N_photonsr = rand();idx = find(cum_prob >= r, 1);photon_positions(i) = x_fine(idx);
end%% 方法3: 单光子双缝干涉（量子力学描述）
fprintf('正在进行单光子干涉模拟...\n');N_single_photons = 5000;
single_photon_positions = zeros(1, N_single_photons);% 波函数描述
psi1 = @(x_pos) exp(1i * k * sqrt(D^2 + (x_pos - d/2).^2));
psi2 = @(x_pos) exp(1i * k * sqrt(D^2 + (x_pos + d/2).^2));% 考虑缝宽
x_slit = linspace(-a/2, a/2, 10);
for i = 1:N_single_photons% 随机选择通过哪个缝（但保持波函数特性）if rand() > 0.5% 通过缝1slit_offset = x_slit(randi(length(x_slit)));effective_x = -d/2 + slit_offset;else% 通过缝2slit_offset = x_slit(randi(length(x_slit)));effective_x = d/2 + slit_offset;end% 计算概率幅x_test = linspace(-L/2, L/2, 1000);prob_amplitude = exp(1i * k * sqrt(D^2 + (x_test - effective_x).^2));probability = abs(prob_amplitude).^2;probability = probability / sum(probability);% 根据概率分布选择探测位置cum_prob = cumsum(probability);r = rand();idx = find(cum_prob >= r, 1);single_photon_positions(i) = x_test(idx);
end%% 可视化结果
figure;% 1. 理论干涉图样
subplot(2,3,1);
plot(x*1000, I_total, 'b-', 'LineWidth', 2);
xlabel('屏幕位置 (mm)');
ylabel('相对光强');
title('双缝干涉理论图样');
grid on;
xlim([-5, 5]);% 添加双缝示意图
hold on;
plot([-d/2*1000, -d/2*1000], [0, 0.1], 'k--', 'LineWidth', 1);
plot([d/2*1000, d/2*1000], [0, 0.1], 'k--', 'LineWidth', 1);
text(-d/2*1000, 0.12, '缝1', 'HorizontalAlignment', 'center');
text(d/2*1000, 0.12, '缝2', 'HorizontalAlignment', 'center');% 2. 光子分布直方图
subplot(2,3,2);
histogram(photon_positions*1000, 100, 'FaceColor', 'blue', 'FaceAlpha', 0.7);
xlabel('屏幕位置 (mm)');
ylabel('光子计数');
title(sprintf('光子分布 (%d 个光子)', N_photons));
grid on;
xlim([-5, 5]);% 3. 单光子累积分布
subplot(2,3,3);
for i = 1:min(1000, N_single_photons)plot(single_photon_positions(1:i)*1000, i, 'b.', 'MarkerSize', 1);if i == 1hold on;end
end
xlabel('屏幕位置 (mm)');
ylabel('光子数');
title('单光子累积分布');
grid on;
xlim([-5, 5]);% 4. 单光子分布直方图
subplot(2,3,4);
histogram(single_photon_positions*1000, 100, 'FaceColor', 'red', 'FaceAlpha', 0.7);
xlabel('屏幕位置 (mm)');
ylabel('光子计数');
title(sprintf('单光子干涉 (%d 个光子)', N_single_photons));
grid on;
xlim([-5, 5]);% 5. 理论 vs 实验对比
subplot(2,3,5);
% 理论曲线
plot(x*1000, I_total, 'b-', 'LineWidth', 2, 'DisplayName', '理论预测');
hold on;
% 光子计数直方图
[counts, edges] = histcounts(photon_positions*1000, 50);
centers = (edges(1:end-1) + edges(2:end)) / 2;
counts_normalized = counts / max(counts);
plot(centers, counts_normalized, 'ro-', 'LineWidth', 1, 'MarkerSize', 3, 'DisplayName', '光子模拟');
xlabel('屏幕位置 (mm)');
ylabel('归一化强度');
title('理论与模拟对比');
legend('Location', 'best');
grid on;
xlim([-5, 5]);% 6. 实验装置示意图
subplot(2,3,6);
% 绘制光源
rectangle('Position', [-0.1, 0.4, 0.2, 0.1], 'FaceColor', 'yellow', 'EdgeColor', 'black');
text(0, 0.52, '光源', 'HorizontalAlignment', 'center');% 绘制双缝
hold on;
rectangle('Position', [-0.02, 0.2, 0.04, 0.15], 'FaceColor', [0.5 0.5 0.5]);
line([-d/2, -d/2], [0.2, 0.35], 'Color', 'white', 'LineWidth', 3);
line([d/2, d/2], [0.2, 0.35], 'Color', 'white', 'LineWidth', 3);
text(0, 0.15, '双缝板', 'HorizontalAlignment', 'center');% 绘制屏幕
line([-L/2, L/2], [0, 0], 'Color', 'black', 'LineWidth', 2);
text(0, -0.05, '观察屏幕', 'HorizontalAlignment', 'center');% 绘制光路
plot([0, -d/2], [0.4, 0.35], 'r--', 'LineWidth', 1);
plot([0, d/2], [0.4, 0.35], 'r--', 'LineWidth', 1);
plot([-d/2, -2e-3], [0.35, 0], 'r--', 'LineWidth', 1);
plot([d/2, 2e-3], [0.35, 0], 'r--', 'LineWidth', 1);axis equal;
xlim([-0.1, 0.1]);
ylim([-0.1, 0.6]);
title('实验装置示意图');
set(gca, 'XTick', [], 'YTick', []);sgtitle('双缝干涉实验模拟', 'FontSize', 14, 'FontWeight', 'bold');%% 计算干涉条纹参数
% 找到极大值位置
[peaks, peak_locs] = findpeaks(I_total, 'MinPeakHeight', 0.5);
peak_positions = x(peak_locs);if length(peak_positions) >= 2fringe_spacing = mean(diff(peak_positions));theoretical_spacing = lambda * D / d;fprintf('\n=== 干涉条纹分析 ===\n');fprintf('观测到的条纹间距: %.3f mm\n', fringe_spacing*1000);fprintf('理论条纹间距: %.3f mm\n', theoretical_spacing*1000);fprintf('相对误差: %.2f%%\n', abs(fringe_spacing - theoretical_spacing)/theoretical_spacing*100);
end%% 交互式参数探索
fprintf('\n=== 交互式探索 ===\n');
fprintf('可以修改以下参数观察干涉图样变化:\n');
fprintf('1. 波长 lambda (当前: %.1f nm)\n', lambda*1e9);
fprintf('2. 双缝间距 d (当前: %.2f mm)\n', d*1e3);
fprintf('3. 缝宽 a (当前: %.2f mm)\n', a*1e3);
fprintf('4. 屏幕距离 D (当前: %.1f m)\n', D);