Nhiễu xạ hai chùm sáng qua lỗ tròn

Giới thiệu

Chương trình được viết trên Matlab, mô phỏng hiệu ứng nhiễu xạ khi cho hai chùm sáng đơn sắc, đặt lệch góc nhau. Mỗi chùm sáng sẽ tạo nên một vết tròn nhiễu xạ trên màn.

Hình 1: Nhiễu xạ hai chùm sáng qua lỗ tròn – đồ thị phân bố cường độ

Hình 2: Nhiễu xạ hai chùm sáng qua lỗ tròn – hình ảnh nhiễu xạ trên màn

Code chương trình Matlab

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121

function Diffraction_2circ_Aperture
% Author: Tran Hai Cat
% Lecturer in Physics, HCM University of Technology and Education
% - Dai hoc Su pham Ky thuat Tp. Ho Chi Minh
% Created: 2019.05.26
clc
clear variables
close all

%% Hien thi:
show = 1; % 1:3D graph, 0:intensive

%% Nguon:
lambda = 0.6e-3;
k = 2*pi/lambda;
% Goc toi cua song:
alpha = [-pi/50,pi/50];

%% Tham so lo tron:
R = 1e-2; % ban kinh lo tron
L = 10e-2; % khoang cach den man
N_R = 50;
N_phi = 100;
r = linspace(0,R,N_R);
phi = linspace(0,2*pi,N_phi);

%% Calculation:
r_screeen = 2e-2;
xmin = -r_screeen;
xmax = r_screeen;
ymin = -r_screeen;
ymax = r_screeen;

Nx = 100;
Ny = 100;

x = linspace(xmin,xmax,Nx);
y = linspace(ymin,ymax,Ny);
[X,Y] = meshgrid(x,y);
Z_map = zeros(size(X));
E = zeros(size(X));

n_in_1 = [sin(alpha(1)),0,-cos(alpha(1))];
n_in_2 = [sin(alpha(2)),0,-cos(alpha(2))];
for i = 1:Ny
for j = 1:Nx
r0 = [x(j),y(i),-L];
len_r0 = sqrt(sum(r0.*r0));
n0 = r0/len_r0;
for i_lo = 1:N_R
for j_lo = 1:N_phi
x_lo = r(i_lo)*cos(phi(j_lo));
y_lo = r(i_lo)*sin(phi(j_lo));
r_lo = [x_lo,y_lo,0];
delta_r_1 = -sum(r_lo.*n_in_1)+sum(r_lo.*n0);
delta_r_2 = -sum(r_lo.*n_in_2)+sum(r_lo.*n0);
E(i,j) = E(i,j)+exp(-1i*k*delta_r_1)+exp(-1i*k*delta_r_2);
end
end
end
end
I = abs(E).^2;
Imax = max(I(:));
I = 0.3*I/Imax*L;

figure('name','Circular Hole Diffration','color','w','numbertitle','off');
if show == 0
img = getimage(imagesc(x,y,I));
img = 255*img/max(img(:));
wp = warp(Z_map,X,Y,img,colormap(parula));
wp.AlphaData = 1;
else
H = mesh(I,X,Y);
set(H, 'CData', get(H, 'XData'));
H.AlphaData = 0.8;
end
hold on
% Display hole:
x_lo = R*cos(phi);
y_lo = R*sin(phi);
z_lo = L*ones(size(x_lo));
plot3(z_lo,x_lo,y_lo,'linewidth',2)

% Display source:
N_R = 5;
N_phi = 15;
r = linspace(0,R,N_R);
phi = linspace(0,2*pi,N_phi);
n_in_1_draw = n_in_1/50;
n_in_2_draw = n_in_2/50;
for i_lo = 1:N_R
for j_lo = 1:N_phi
x2 = r(i_lo)*cos(phi(j_lo));
y2 = r(i_lo)*sin(phi(j_lo));
z2 = L;

x1 = x2-n_in_1_draw(1);
y1 = y2-n_in_1_draw(2);
z1 = L-n_in_1_draw(3);

plot3([z1,z2],[x1,x2],[y1,y2],...
'linewidth',0.5,'color','b');

z1 = L-n_in_2_draw(3);
x1 = x2-n_in_2_draw(1);
y1 = y2-n_in_2_draw(2);

plot3([z1,z2],[x1,x2],[y1,y2],...
'linewidth',0.5,'color','r');
end
end

axis equal
axis([-inf inf xmin xmax ymin ymax]);
rotate3d on
view(150,30);
grid off
xlabel('X');
ylabel('Y');
zlabel('Z');
% set(gca,'color','k','xcolor','w','ycolor','w','zcolor','w')