Chùm hạt với hố 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
function de_Broglie_wave_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;
U0 = 10*q_element;

xmin = -20*Angstrom;
xmax = 20*Angstrom;
x_step = 0*Angstrom;

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

[psi_incident,psi_reflective,psi] = step(E,U0,x_step,x,t);

%% FIGURE
figure('name','Wavepaket - Particle Flow','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_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 -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_incident,psi_reflective,psi] = step(E,U0,x_step,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',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 [psi_incident,psi_reflective,psi] = step(E,U0,x_step,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;

B1 = (k1-k2)/(k1+k2)*exp(2*j*k1*x_step);
S2 = 2*k1/(k1+k2)*exp(j*(k1-k2)*x_step);

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

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