-
Notifications
You must be signed in to change notification settings - Fork 0
/
phenomenon_fit.m
172 lines (153 loc) · 5.56 KB
/
phenomenon_fit.m
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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
% reflection of a bloch state by a barrier; numerics vs analytic
% 2016.01.11
clear all; close all; clc; myfont = 22;
step = 0.001;
nPeriod = 11.5;
Llist = [200, 200, 100];
kilist = [80, 100, 50];
Ulist = [1.5, 2, 12];
Tlist = zeros(1, 3);
plist = zeros(6, nPeriod/step + 1);
plist2 = plist;
envolop = zeros(3, nPeriod/step + 1);
for sw = 1: 3
L = Llist(sw); N = 2*L+1;
ki = kilist(sw);
U = Ulist(sw);
delta = 2*sin(2*pi/N)*sin(2*pi*ki/N);
g = U/N;
T = 2*pi/delta;
Tlist(sw) = T;
dt = 0.001*T;
Tmax = nPeriod*T;
tlist = 0:dt:Tmax;
theta = 2*atan(g*T);
omega = theta /T;
envolop(sw,:) = abs(1 - exp(-i*omega*tlist)).^2/4;
xlist = -L:L;
xlist = xlist';
psi0 = (1/sqrt(N))*exp(i*(2*pi*ki/N)*xlist);
psif = (1/sqrt(N))*exp(i*(-2*pi*ki/N)*xlist);
% hamiltonian
H = zeros(N, N);
for s= 1:(N-1)
H(s,s+1) = -1; H(s+1,s) = -1;
end
H(1,N) = -1; H(N,1) = -1;
H(L+1, L+1) = U;
[VV,DD] = eig(H);
dd = diag(DD);
psi1 = VV'*psi0;
for s = 1: length(tlist)
time = (s-1)*dt;
psi = VV*(exp(-i*time*dd).*psi1);
plist(1+ (sw-1)*2,s) = abs(psi'*psi0)^2;
plist(2+ (sw-1)*2,s) = abs(psi'*psif)^2;
n1 = floor(time/T);
time2 = time - n1*T;
amp1 = 1/2;
amp2 = 0.5*(((1-i*2*g*pi/delta)/(1+i*2*g*pi/delta))^n1)*(1-i*2*g*(time2- pi/delta))/(1+ i*2*g*pi/delta);
plist2(1+ (sw-1)*2,s) = abs( amp1 + amp2 )^2;
plist2(2+ (sw-1)*2,s) = abs(-amp1 + amp2 )^2;
end
end
% h1 = figure;
% plot(tlist/T, plist, tlist/T, plist2, ':')
% set(gca,'fontsize',myfont)
% xlim([0 nPeriod])
% ylim([0 1])
% xlabel('t/T','fontsize',myfont);
% ylabel('p','fontsize',myfont);
% str = strcat ('U=', num2str(U),', N=',num2str(N),', ki=',num2str(ki));
% title(str,'fontsize',myfont)
h1= figure;
ha = tight_subplot(1,3,[.05 .05],[.7 .01],[.1 .01])
aa = 0.95;
bb = 0.9;
myfont = 12;
axes(ha(1))
plot((0:step:nPeriod)*Tlist(1), plist(1, :), (0:step:nPeriod)*Tlist(1), plist(2, :),':','linewidth',1.5)
xlabel('$t$','fontsize',myfont, 'Interpreter','Latex')
ylabel('$ P_i \;\& \; P_r $','fontsize',myfont,'Interpreter','Latex')
set(gca,'fontsize',myfont)
xlim([0 nPeriod*Tlist(1)])
ylim([0 1])
XL=xlim; YL=ylim;
str=strcat('(a)');
% text(0.5*XL(1)+0.5*XL(2), 0.5*YL(2)+0.5*YL(1),str,'fontsize',myfont)
text(aa*XL(1)+(1-aa)*XL(2), bb*YL(2)+(1-bb)*YL(1),str,'fontsize',myfont)
% h2=figure;
axes(ha(2))
% set(gca,'position',gcaposition);
plot((0:step:nPeriod)*Tlist(2), plist(3, :), (0:step:nPeriod)*Tlist(2), plist(4, :),':','linewidth',1.5)
xlabel('$t$','fontsize',myfont, 'Interpreter','Latex')
% ylabel('$ |\eta|^2 $','fontsize',myfont,'Interpreter','Latex')
set(gca,'fontsize',myfont)
xlim([0 nPeriod*Tlist(2)])
ylim([0 1])
XL=xlim; YL=ylim;
str=strcat('(b)');
text(aa*XL(1)+(1-aa)*XL(2), bb*YL(2)+(1-bb)*YL(1),str,'fontsize',myfont)
% h3=figure;
% set(gca,'position',gcaposition);
axes(ha(3))
plot((0:step:nPeriod)*Tlist(3), plist(5, :), (0:step:nPeriod)*Tlist(3), plist(6, :),':','linewidth',1.5)
xlabel('$t$','fontsize',myfont, 'Interpreter','Latex')
% ylabel('$ |\eta|^2 $','fontsize',myfont,'Interpreter','Latex')
set(gca,'fontsize',myfont)
xlim([0 nPeriod*Tlist(3)])
ylim([0 1])
XL=xlim; YL=ylim;
str=strcat('(c)');
text(aa*XL(1)+(1-aa)*XL(2), bb*YL(2)+(1-bb)*YL(1),str,'fontsize',myfont)
str = strcat('pipf1.eps');
print(h1,'-depsc',str)
h1= figure;
ha = tight_subplot(1,3,[.05 .05],[.7 .01],[.1 .01])
myfont = 12;
axes(ha(1))
plot((0:step:nPeriod)*Tlist(1), plist(1, :), (0:step:nPeriod)*Tlist(1), plist(2, :),':','linewidth',1.5)
hold on
plot((0:step:nPeriod)*Tlist(1), plist2(1, :),':', (0:step:nPeriod)*Tlist(1), plist2(2, :),'linewidth',1.5)
% plot((0:step:nPeriod)*Tlist(1),envolop(1,:),'r--','linewidth',1.5)
xlabel('$t$','fontsize',myfont, 'Interpreter','Latex')
ylabel('$ P_i \;\& \; P_r $','fontsize',myfont,'Interpreter','Latex')
set(gca,'fontsize',myfont)
xlim([0 nPeriod*Tlist(1)])
ylim([0 1])
XL=xlim; YL=ylim;
str=strcat('(a)');
% text(0.5*XL(1)+0.5*XL(2), 0.5*YL(2)+0.5*YL(1),str,'fontsize',myfont)
text(aa*XL(1)+(1-aa)*XL(2), bb*YL(2)+(1-bb)*YL(1),str,'fontsize',myfont)
% h2=figure;
axes(ha(2))
% set(gca,'position',gcaposition);
plot((0:step:nPeriod)*Tlist(2), plist(3, :), (0:step:nPeriod)*Tlist(2), plist(4, :),':','linewidth',1.5)
hold on
plot((0:step:nPeriod)*Tlist(2), plist2(3, :),':', (0:step:nPeriod)*Tlist(2), plist2(4, :),'linewidth',1.5)
% plot((0:step:nPeriod)*Tlist(2),envolop(2,:),'r--','linewidth',1.5)
xlabel('$t$','fontsize',myfont, 'Interpreter','Latex')
% ylabel('$ |\eta|^2 $','fontsize',myfont,'Interpreter','Latex')
set(gca,'fontsize',myfont)
xlim([0 nPeriod*Tlist(2)])
ylim([0 1])
XL=xlim; YL=ylim;
str=strcat('(b)');
text(aa*XL(1)+(1-aa)*XL(2), bb*YL(2)+(1-bb)*YL(1),str,'fontsize',myfont)
% h3=figure;
% set(gca,'position',gcaposition);
axes(ha(3))
plot((0:step:nPeriod)*Tlist(3), plist(5, :), (0:step:nPeriod)*Tlist(3), plist(6, :),':','linewidth',1.5)
hold on
plot((0:step:nPeriod)*Tlist(3), plist2(5, :),':', (0:step:nPeriod)*Tlist(3), plist2(6, :),'linewidth',1.5)
% plot((0:step:nPeriod)*Tlist(3),envolop(3,:),'r--','linewidth',1.5)
xlabel('$t$','fontsize',myfont, 'Interpreter','Latex')
% ylabel('$ |\eta|^2 $','fontsize',myfont,'Interpreter','Latex')
set(gca,'fontsize',myfont)
xlim([0 nPeriod*Tlist(3)])
ylim([0 1])
XL=xlim; YL=ylim;
str=strcat('(c)');
text(aa*XL(1)+(1-aa)*XL(2), bb*YL(2)+(1-bb)*YL(1),str,'fontsize',myfont)
str = strcat('pipf2.eps');
print(h1,'-depsc',str)