Bó sóng với rào thế bậc thang

  • 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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
function wavepacket_step_potential_barrier
% Created by Tran Hai Cat
% 2018.05.01
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;
x_step = 0*Angstrom;
U0 = 10*q_element;

xmin = -200*Angstrom;
xmax = 200*Angstrom;

x0 = -150*Angstrom; % Vi tri ban dau cua bo song
delta_x = 10*Angstrom; % Do rong bo song

n_k = 50;

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

C = packet(x0,delta_x,n_k);
psi = step(C,delta_x,n_k,E,U0,x_step,x,t);

%% FIGURE
figure('name','Wavepaket Barrier','color','black','numbertitle','off');
hold on
plot(x/Angstrom,U_step(U0,x_step,x)/q_element,'w','linewidth',2);
plot(x/Angstrom,Energy(E,x)/q_element,'b','linewidth',2);
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 = step(C,delta_x,n_k,E,U0,x_step,x,t);
set(line_psi_psi,'ydata',psi.*conj(psi));
set(line_psi,'ydata',U0/q_element+real(psi));
pause(0.001);
end

function y = U_step(U0,x_step,x)
lenx = length(x);
y = zeros(1,lenx);
for i = 1:lenx
if(x(i)<=x_step)
y(i) = 0;
else
y(i) = U0;
end
end

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

function C = packet(x0,delta_x,n_k)
dk = 1/delta_x/n_k;
j = sqrt(-1);
sum = 0;
C = zeros(1,2*n_k+1);
C(n_k+1) = 1;
for ik = 1:n_k
delta_k = -ik*dk;
C(n_k+1-ik) = exp(-delta_k.^2*delta_x.^2/4)*exp(-j*delta_k*x0);
sum = sum+C(n_k+1-ik)*conj(C(n_k+1-ik));

delta_k = ik*dk;
C(n_k+1+ik) = exp(-delta_k.^2*delta_x.^2/4)*exp(-j*delta_k*x0);
sum = sum+C(n_k+1+ik)*conj(C(n_k+1+ik));
end
norm = dk*sqrt(1./sum);
for i = 1:2*n_k+1
C(i) = C(i)*norm*1.5e-8;
end

function psi = step(C,delta_x,n_k,E,U0,x_step,x,t)
global h m
j = sqrt(-1);
lenx = length(x);
dk = 1/delta_x/n_k;
k0 = sqrt(2*m*E)/h;
psi = zeros(1,lenx);

k1 = k0;
omega = 0.5*k1*k1*h/m;
k2 = sqrt(2*m*(E-U0))/h;
B = (k1-k2)/(k1+k2)*exp(2*j*k1*x_step);
S = 2*k1/(k1+k2)*exp(j*(k1-k2)*x_step);

for i = 1:lenx
if(x(i)<=x_step)
psi(i) = C(n_k+1)*(exp(j*k1.*x(i)-j*omega*t)+B*exp(-j*k1.*x(i)-j*omega*t));
else
psi(i) = C(n_k+1)*S*exp(j*k2.*x(i)-j*omega*t);
end
end

for ik = 1:n_k
k1 = k0-ik*dk;
omega = 0.5*k1*k1*h/m;
E = h*omega;
k2 = sqrt(2*m*(E-U0))/h;
B = (k1-k2)/(k1+k2)*exp(2*j*k1*x_step);
S = 2*k1/(k1+k2)*exp(j*(k1-k2)*x_step);

for i = 1:lenx
if(x(i)<=x_step)
psi(i) = psi(i)+C(n_k+1-ik)*(exp(j*k1.*x(i)-j*omega*t)+B*exp(-j*k1.*x(i)-j*omega*t));
else
psi(i) = psi(i)+C(n_k+1-ik)*S*exp(j*k2.*x(i)-j*omega*t);
end
end

k1 = k0+ik*dk;
omega = 0.5*k1*k1*h/m;
E = h*omega;
k2 = sqrt(2*m*(E-U0))/h;
B = (k1-k2)/(k1+k2)*exp(2*j*k1*x_step);
S = 2*k1/(k1+k2)*exp(j*(k1-k2)*x_step);

for i = 1:lenx
if(x(i)<=x_step)
psi(i) = psi(i)+C(n_k+1+ik)*(exp(j*k1.*x(i)-j*omega*t)+B*exp(-j*k1.*x(i)-j*omega*t));
else
psi(i) = psi(i)+C(n_k+1+ik)*S*exp(j*k2.*x(i)-j*omega*t);
end
end
end