Nhiễu xạ sóng qua lỗ tròn – Vật lý mô phỏ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 tròn.

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

Hình 2: Nhiễu xạ 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

function Diffraction_Circ_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:
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));

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 = sum(r_lo.*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','Circular 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:
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)

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

figure('name','Circular 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:
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)

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