Chùm hạt đi qua hố thế

  • Code chương trình

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
function de_Broglie_wave_through_square_hole
% Created by Tran Hai Cat
% 2018.05.15
clc;
clear variables
close all

global h m

%% CONSTANTS
h = 1.054e-34;
m = 9.1095e-31;
q_element = 1.6e-19;
Angstrom = 1e-10;

%% INPUT PARAMETRS
E = 5*q_element;
U0 = 10*q_element;

xmin = -20*Angstrom;
xmax = 20*Angstrom;
x1 = 0*Angstrom;
x2 = 5*Angstrom;

%% data processing
x = linspace(xmin,xmax,500);
omega = E/h;
T = 2*pi/omega;
dt = T/500;
tmax = 100*T;
t = 0;

psi = barrier(E,U0,x1,x2,x,t);

%% FIGURE
figure('name','Wavepaket - Particle Flow','color','black','numbertitle','off');
hold on
plot(x/Angstrom,U_barrier(U0,x1,x2,x)/q_element,'w','linewidth',2);
plot(x/Angstrom,Energy(E,x)/q_element,'b','linewidth',2);
% line_psi_incident = plot(psi_incident(1,:)/Angstrom,real(psi_incident(2,:)),'linewidth',1,'color','b');
% line_psi_reflective = plot(psi_reflective(1,:)/Angstrom,real(psi_reflective(2,:)),'linewidth',1,'color','magenta');
line_psi = plot(x/Angstrom,real(psi),'linewidth',1,'color','y');
line_psi_psi = plot(x/Angstrom,psi.*conj(psi),'linewidth',2,'color',[1,0.5,0]);
axis([xmin/Angstrom xmax/Angstrom -U0/q_element U0/q_element]);
set(gca,'color','k','xcolor','w','ycolor','w')
xlabel('x [Angstrom]');
ylabel('E, U [eV]');

%% CALCULATION
while t<tmax
t = t+dt;
psi = barrier(E,U0,x1,x2,x,t);
% set(line_psi_incident,'ydata',U0/q_element+real(psi_incident(2,:)));
% set(line_psi_reflective,'ydata',U0/q_element+real(psi_reflective(2,:)));
set(line_psi_psi,'ydata',-U0/q_element+psi.*conj(psi));
set(line_psi,'ydata',real(psi));
pause(0.001);
end

function y = U_barrier(U0,x1,x2,x)
lenx = length(x);
y = zeros(1,lenx);
for i = 1:lenx
if(x(i)<=x1)||(x(i)>=x2)
y(i) = 0;
else
y(i) = -U0;
end
end

function y = Energy(E,x)
y = E*ones(size(x));

function psi = barrier(E,U0,x1,x2,x,t)
global h m
j = sqrt(-1);
lenx = length(x);
psi = zeros(1,lenx);

omega = E/h;
k1 = sqrt(2*m*E)/h;
k2 = sqrt(2*m*(E+U0))/h;

A_matrix = [...
exp(-j*k1*x1) -exp(j*k2*x1) -exp(-j*k2*x1) 0
-k1*exp(-j*k1*x1) -k2*exp(j*k2*x1) k2*exp(-j*k2*x1) 0
0 exp(j*k2*x2) exp(-j*k2*x2) -exp(j*k1*x2)
0 k2*exp(j*k2*x2) -k2*exp(-j*k2*x2) -k1*exp(j*k1*x2)
];
b_vector = [...
-exp(j*k1*x1)
-k1*exp(j*k1*x1)
0
0
];
Koef = A_matrix\b_vector;
B = Koef(1);
C = Koef(2);
D = Koef(3);
F = Koef(4);

for i = 1:lenx
if(x(i)<=x1)
psi(i) = exp(j*k1.*x(i)-j*omega*t)...
+B*exp(-j*k1.*x(i)-j*omega*t);
elseif(x(i)<=x2)
psi(i) = C*exp(j*k2.*x(i)-j*omega*t)...
+D*exp(-j*k2.*x(i)-j*omega*t);
else
psi(i) = F*exp(j*k1.*x(i)-j*omega*t);
end
end