Chùm hạt với rào thế chữ nhật

  • 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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
function de_Broglie_wave_rectangular_potential_barrier
% Created by Tran Hai Cat
% 2018.05.14
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 = 15*q_element;
U0 = 10*q_element;

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

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

[psi_in1,psi_ref1,psi_in2,psi_ref2,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_incident1 = plot(psi_in1(1,:)/Angstrom,real(psi_in1(2,:)),'linewidth',1,'color','b');
line_psi_reflective1 = plot(psi_ref1(1,:)/Angstrom,real(psi_ref1(2,:)),'linewidth',1,'color','magenta');
line_psi_incident2 = plot(psi_in2(1,:)/Angstrom,real(psi_in2(2,:)),'linewidth',1,'color','b');
line_psi_reflective2 = plot(psi_ref2(1,:)/Angstrom,real(psi_ref2(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 -0 2*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_in1,psi_ref1,psi_in2,psi_ref2,psi] = barrier(E,U0,x1,x2,x,t);
set(line_psi_incident1,'ydata',U0/q_element+real(psi_in1(2,:)));
set(line_psi_reflective1,'ydata',U0/q_element+real(psi_ref1(2,:)));
set(line_psi_incident2,'ydata',U0/q_element+real(psi_in2(2,:)));
set(line_psi_reflective2,'ydata',U0/q_element+real(psi_ref2(2,:)));
set(line_psi_psi,'ydata',psi.*conj(psi));
set(line_psi,'ydata',U0/q_element+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_in1,psi_ref1,psi_in2,psi_ref2,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;
B1 = Koef(1); % S1 = 1
S2 = Koef(2);
B2 = Koef(3);
S3 = Koef(4); % B3 = 0;

psi_in1 = [x(1);exp(j*k1.*x(1)-j*omega*t)];
psi_ref1 = [x(1);B1*exp(-j*k1.*x(1)-j*omega*t)];
psi_in2 = [x1;S2*exp(j*k2.*x1-j*omega*t);];
psi_ref2 = [x1;B2*exp(-j*k2.*x1-j*omega*t)];

for i = 1:lenx
if(x(i)<=x1)
psi_incident_i = exp(j*k1.*x(i)-j*omega*t);
psi_reflective_i = B1*exp(-j*k1.*x(i)-j*omega*t);
psi_in1 = [psi_in1 [x(i);psi_incident_i]];
psi_ref1 = [psi_ref1 [x(i);psi_reflective_i]];

psi(i) = psi_incident_i+psi_reflective_i;
elseif(x(i)<=x2)
psi_incident_i = S2*exp(j*k2.*x(i)-j*omega*t);
psi_reflective_i = B2*exp(-j*k2.*x(i)-j*omega*t);
psi_in2 = [psi_in2 [x(i);psi_incident_i]];
psi_ref2 = [psi_ref2 [x(i);psi_reflective_i]];

psi(i) = psi_incident_i+psi_reflective_i;
else
psi(i) = S3*exp(j*k1.*x(i)-j*omega*t);
end
end