-
Notifications
You must be signed in to change notification settings - Fork 0
/
test1.py
189 lines (165 loc) · 7.58 KB
/
test1.py
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib import rc
from matplotlib.ticker import ScalarFormatter, FormatStrFormatter
from qutip import *
plt.figure(figsize=(10, 10))
plt.rc('font', family='serif')
rc('text', usetex=False)
ax = plt.gca()
ax.get_yaxis().get_major_formatter().set_useOffset(False)
directory = 'G:\Projects\Fluxonium\Data\Fluxonium #10_New software'
measurement = 'Two_tone_spec_YOKO_41to43mA_Qubit_3to4GHz_5dBm_Cav_10.3039GHz_8dBm_IF_0.05GHz_measTime_500ns_avg_50000'
path = directory + '\\' + measurement
# Read data
current = np.genfromtxt(path + '_CURR.dat')
current = current[1::]-0.04
freq = np.genfromtxt(path + '_FREQ.dat')
freq = freq[1::]
data = np.genfromtxt(path + '_PHASEMAG.dat')
phase = data[1::, 0] # phase is recorded in rad
phase = phase #
mag = data[1::, 0]
# plt.figure(1)
Z = np.zeros((len(current), len(freq)))
for idx in range(len(current)):
temp = np.unwrap(phase[idx * len(freq):(idx + 1) * len(freq)])
Z[idx, :] = temp - np.average(temp)
Z = Z * 180 / (np.pi)
X, Y = np.meshgrid(current, freq[0:len(freq) / 2 + 2])
Z1 = Z.transpose()[0:len(freq) / 2 + 2]
plt.pcolormesh(X, Y, Z1, cmap='GnBu_r', vmin=-3, vmax=-1.5)
X, Y = np.meshgrid(current, freq[len(freq) / 2 + 2:len(freq) - 1])
Z2 = Z.transpose()[len(freq) / 2 + 2:len(freq) - 1]
plt.pcolormesh(X, Y, Z2, cmap='GnBu_r', vmin=-4, vmax=1)
# Define constants
e = 1.602e-19 # Fundamental charge
h = 6.62e-34 # Placnk's constant
phi_o = h / (2 * e) # Flux quantum
"""
First section of the script attempts to plot the energies vs external flux
"""
# Hamiltonian definition
def Ho(N, E_l, E_c, E_j_sum, d, phi_squid, phi_ext):
E_j1 = 0.5 * E_j_sum * (1 + d)
E_j2 = 0.5 * E_j_sum * (1 - d)
a = tensor(destroy(N))
mass = 1.0 / (8.0 * E_c)
w = sqrt(8.0 * E_c * E_l)
phi = (a + a.dag()) * (8 * E_c / E_l) ** (0.25) / np.sqrt(2)
na = 1j * (a.dag() - a) * (E_l / (8 * E_c)) ** (0.25) / np.sqrt(2)
ope1 = 1j * (-phi + phi_ext)
ope2 = 1j * (
phi - phi_ext + phi_squid) # phi_squid and phi_ext here are the external phases, or normalized flux, = flux*2pi/phi_o
H = 4.0 * E_c * na ** 2 + 0.5 * E_l * (phi) ** 2 - 0.5 * E_j1 * (ope1.expm() + (-ope1).expm()) - 0.5 * E_j2 * (
ope2.expm() + (-ope2).expm())
return H.eigenenergies()
def coupled_H(Na, E_l, E_c, E_j_sum, d, phi_squid, phi_ext, Nr, wr, g):
E_j1 = 0.5 * E_j_sum * (1 + d)
E_j2 = 0.5 * E_j_sum * (1 - d)
a = tensor(destroy(Na), qeye(Nr))
b = tensor(qeye(Na), destroy(Nr))
phi = (a + a.dag()) * (8.0 * E_c / E_l) ** (0.25) / np.sqrt(2.0)
na = 1.0j * (a.dag() - a) * (E_l / (8 * E_c)) ** (0.25) / np.sqrt(2.0)
ope1 = 1.0j * (phi_ext - phi)
ope2 = 1.0j * (phi + phi_squid - phi_ext)
H_f = 4.0 * E_c * na ** 2 + 0.5 * E_l * (phi) ** 2 - 0.5 * E_j1 * (ope1.expm() + (-ope1).expm()) - 0.5 * E_j2 * (
ope2.expm() + (-ope2).expm())
H_r = wr * (b.dag() * b + 1.0 / 2)
H_c = -2 * g * na * (b.dag + b)
H = H_f + H_r + H_c
return H.eigenenergies()
def trans_energies(N, E_l, E_c, E_j_sum, d, A_j, A_c, B_coeff, beta_squid, beta_ext, level_num, current, iState):
B_field = current * B_coeff * 1e-4 # in T, this depends on a seperate measurement
phi_squid = B_field * A_j # these are flux, not normalized
phi_ext = B_field * A_c
trans_energy = np.zeros((level_num - iState, len(phi_ext)))
for idx in range(len(phi_ext)):
energies = Ho(N, E_l, E_c, E_j_sum, d, 2 * np.pi * (phi_squid[idx] / phi_o - beta_squid),
2 * np.pi * (phi_ext[idx] / phi_o - beta_ext)) # normalize the flux -> phase here
for level in range(iState + 1, level_num):
trans_energy[level - iState, idx] = energies[level] - energies[iState]
return trans_energy
def coupled_trans_energies(N, E_l, E_c, E_j_sum, d, A_j, A_c, B_coeff, beta_squid, beta_ext, level_num, current, iState,
Nr, wr, g):
B_field = current * B_coeff * 1e-4 # in T, this depends on a seperate measurement
phi_squid = B_field * A_j # these are flux, not normalized
phi_ext = B_field * A_c
trans_energy = np.zeros((level_num - iState, len(phi_ext)))
for idx in range(len(phi_ext)):
energies = coupled_H(N, E_l, E_c, E_j_sum, d, 2 * np.pi * (phi_squid[idx] / phi_o - beta_squid),
2 * np.pi * (phi_ext[idx] / phi_o - beta_ext), Nr, wr,
g) # normalize the flux -> phase here
for level in range(iState + 1, level_num):
trans_energy[level - iState, idx] = energies[level] - energies[iState]
return trans_energy
########################################################################
# Fitting for bottom spectrum
N=50
E_l=0.722729827116
E_c=0.552669197076
E_j_sum=17.61374383
A_j=4.76321410213e-12
A_c=1.50075181762e-10
d=0.125005274368
beta_squid=0.129912406349
beta_ext=0.356925557542
B_coeff = 60
level_num = 5
current = np.linspace(0.0412, 0.0421, 100)
iState = 0
spectrum = trans_energies(N, E_l, E_c, E_j_sum, d, A_j, A_c, B_coeff, beta_squid, beta_ext, level_num, current, iState)
for idx in range(iState, 3):
line = plt.plot(current * 1e3, spectrum[idx, :]) # transition from state (iState)
plt.setp(line, linewidth=3.0, linestyle='--', color="black")#, alpha=0.2)
# line = plt.plot(current, spectrum[idx,:]+10.304) # transition from state (iState)
# plt.setp(line,linewidth=2.0, linestyle ='--', color = "black", alpha=0.5)
# line = plt.plot(current, -spectrum[idx,:]+10.304) # transition from state (iState)
# plt.setp(line,linewidth=2.0, linestyle ='--', color = "black", alpha=0.5)
# iState = 1
# spectrum = trans_energies(N, E_l, E_c, E_j_sum, d, A_j, A_c, B_coeff, beta_squid, beta_ext, level_num, current, iState)
# for idx in range(iState,level_num):
# line = plt.plot(current*1e3, spectrum[idx-iState,:]) # transition from state (iState)
# plt.setp(line,linewidth=1.0, linestyle ='--', color = "red", alpha=0.5)
# line = plt.plot(current, spectrum[idx-iState,:]+10.304) # transition from state (iState)
# plt.setp(line,linewidth=2.0, linestyle ='--', color = "red", alpha=0.5)
# line = plt.plot(current, -spectrum[idx-iState,:]+10.304) # transition from state (iState)
# plt.setp(line,linewidth=2.0, linestyle ='-.', color = "red", alpha=0.5)
'''
#Coupled Transition energy calculation and fitting for top spectrum here
N = 40
Nr = 10
E_l=0.735773762652
E_c=0.537375025825
E_j_sum=22.3
A_j=3.83424869313e-12
A_c=1.46689233147e-10
d=0.185865262485
beta_squid=-2.58488114861e-05
beta_ext=-0.0251115059548
B_coeff = 60
g=0.0845608058905
wr = 10.304
current = np.linspace(0.02,0.03,100)
level_num = 7
iState = 0
spectrum = coupled_trans_energies(N, E_l, E_c, E_j_sum, d, A_j, A_c, B_coeff, beta_squid, beta_ext, level_num, current, iState, Nr, wr, g)
for idx in range(iState,level_num):
line = plt.plot(current*1e3, spectrum[idx,:]) # transition from state (iState)
plt.setp(line,linewidth=1.0, linestyle ='-', color = "black", alpha=0.5)
line = plt.plot(current*1e3, spectrum[idx,:]/2) # transition from state (iState)
plt.setp(line,linewidth=1.0, linestyle ='-.', color = "black", alpha=0.9)
iState = 1
spectrum = coupled_trans_energies(N, E_l, E_c, E_j_sum, d, A_j, A_c, B_coeff, beta_squid, beta_ext, level_num, current, iState, Nr, wr, g)
for idx in range(iState,level_num):
line = plt.plot(current*1e3, spectrum[idx-iState,:]) # transition from state (iState)
plt.setp(line,linewidth=1.0, linestyle ='--', color = "red", alpha=0.5)
'''
plt.ylim([3.1, 4.0])
plt.xlim([41.3,42.1])
plt.xticks([41.3,41.7, 42.1])
plt.yticks([3.2,3.6,4])
plt.tick_params(labelsize=26)
# plt.colorbar()
plt.show()