Nhiễu xạ sóng trên mép vật cản

Hình ảnh mô phỏng nhiễu xạ sóng trên rìa mép vật cả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

function diffraction_wall
% 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.08.13
clc
clear variables
close all

%% CONST
c = 3e8; % speed of light [m/s]

%% PARAMETERS
wave_length = 550; % nm

% Source:
Ns = 50;

% Enter dimensions:
Nx = 100; % Number of X-grids
Ny = 200; % Number of Y-grids

pattern_type = 0; % 0 - pcolor, 1 - contourf

xmin = 0; xmax = 10e-6; ymin = -10e-6; ymax = 10e-6;
%% CALCULATION:
lambda = wave_length*1e-9;
T = lambda/c;
w = 2*pi/T;
k = 2*pi/lambda;

x = linspace(xmin,xmax,Nx);
y = linspace(ymin,ymax,Ny);
x_plot = x*1e6;
y_plot = y*1e6;

t = 0;
dt = T/20;

xs = zeros(1,Ns);
ys = linspace(0,ymax,Ns);
% Create figure frame:
fig = figure('name','CALCULATION PROCESS',...
'color','w','numbertitle','off');
set(gcf,'Units','normalized');
set(gcf,'Position',[0 0.05 1 0.85]);
axis equal
colormap jet
hold on

N_movie = 20;
for i_movie = 1:N_movie
t = t+dt;
E = zeros(Nx,Ny);
for i = 1:Nx
for j = 1:Ny
for i_s = 1:Ns
r = sqrt((x(i)-xs(i_s))^2+(y(j)-ys(i_s))^2);
if r>lambda
E(i,j) = E(i,j)+exp(1i*(w*t-k*r))/sqrt(r);
end
end
end
end
E1 = E';
I = real(E1);

figure(fig)
if pattern_type==1
contourf(x_plot,y_plot,I)
else
pcolor(x_plot,y_plot,I)
end
axis image
shading interp;
line([0 0],[0 ymin*1e6],'linewidth',3,'color','b');
title(['Wave length \lambda = ',num2str(wave_length),'nm']);
xlabel('x [{\mu}m]','fontsize',14);
ylabel('y [{\mu}m]','fontsize',14);
colorbar('vertical');
set(gca,'fontsize',14);
drawnow
M(i_movie)=getframe(gcf);
end
fig = figure('name','Diffraction pattern',...
'color','w','numbertitle','off');
set(gcf,'Units','normalized');
set(gcf,'Position',[0 0.05 1 0.85]);

movie(fig,M,40)