Hạt tự do

Hạt tự do trong cơ học lượng tử là một bó sóng được cấu thành từ nhiều sóng de Broglie.

  • 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
function wavepacket_free_particle
% Created by Tran Hai Cat
% 2018.04.30
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 = 20*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;

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

C = packet(x0,delta_x,n_k);
psi = free(C,delta_x,n_k,E,x,t);

%% FIGURE
figure('name','Wavepaket Barrier','color','black','numbertitle','off');
hold on
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 20]);
set(gca,'color','k','xcolor','w','ycolor','w')
xlabel('x [Angstrom]');
ylabel('E, U [eV]');

%% CALCULATION
while t<tmax
t = t+dt;
psi = free(C,delta_x,n_k,E,x,t);
set(line_psi_psi,'ydata',psi.*conj(psi));
set(line_psi,'ydata',10+real(psi));
pause(0.001);
end

function C = packet(x0,delta_x,n_k)
dk = 1/delta_x/n_k;
i_unit = 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(-i_unit*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(-i_unit*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 = free(C,delta_x,n_k,E,x,t)
global h m
i_unit = sqrt(-1);
k0 = sqrt(2*m*E)/h;
dk = 1/delta_x/n_k;
omega = 0.5*k0*k0*h/m;
psi = C(n_k+1).*exp(i_unit*(k0.*x-omega*t));
for ik = 1:n_k
k = k0-ik*dk;
omega = 0.5*k*k*h/m;
psi = psi+C(n_k+1-ik).*exp(i_unit*(k.*x-omega*t));
k = k0+ik*dk;
omega = 0.5*k*k*h/m;
psi = psi+C(n_k+1+ik).*exp(i_unit*(k.*x-omega*t));
end