Nhiễu xạ sóng qua lỗ vuông

Giới thiệu

Chương trình được viết trên Matlab, mô phỏng sự hình thành bức tranh nhiễu xạ khi cho một chùm sáng song song đơn sắc, chiếu vuông góc qua một lỗ hình vuông.

Hình 1: Nhiễu xạ qua lỗ vuông – đồ thị phân bố cường độ

 

Hình 2: Nhiễu xạ qua lỗ vuông – 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
function Diffraction_Square_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

%% Wave Source:
lambda = 0.6e-3;
k = 2*pi/lambda;

%% Aperture's parameter:
b = 1e-2; % do rong khe
L = 10e-2; % khoang cach den man
Nap_x = 50;
Nap_y = 50;
x_ap = linspace(-b/2,b/2,Nap_x);
y_ap = linspace(-b/2,b/2,Nap_y);

%% 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));

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 iap = 1:Nap_x
for jap = 1:Nap_y
r_aperture = [x_ap(iap),y_ap(jap),0];
delta_r = sum(r_aperture.*n0);
E(i,j) = E(i,j)+exp(-1i*k*delta_r);
end
end
end
end
I = abs(E).^2;
Imax = max(I(:));
I = 0.3*I/Imax*L;

figure('name','Square Aperture Diffration','color','w','numbertitle','off');
H = mesh(I,X,Y);
set(H, 'CData', get(H, 'XData'));
H.AlphaData = 0.8;
hold on
% Display source:
plot3([L L L L L],[-b/2 b/2 b/2 -b/2 -b/2],[b/2 b/2 -b/2 -b/2 b/2],...
'linewidth',2)

axis equal
axis([-inf inf xmin xmax ymin ymax]);
rotate3d on
view(150,30);
grid off
xlabel('X');
ylabel('Y');
zlabel('Z');

figure('name','Square Aperture Diffration','color','w','numbertitle','off');
img = getimage(imagesc(x,y,I));
img = 255*img/max(img(:));
wp = warp(Z_map,X,Y,img,colormap(parula));
wp.AlphaData = 1;
hold on
% Display source:
plot3([L L L L L],[-b/2 b/2 b/2 -b/2 -b/2],[b/2 b/2 -b/2 -b/2 b/2],...
'linewidth',2)

axis equal
axis([-inf inf xmin xmax ymin ymax]);
rotate3d on
view(150,30);
grid off
xlabel('X');
ylabel('Y');
zlabel('Z');