54 lines
1.7 KiB
Python
54 lines
1.7 KiB
Python
#!/bin/env python3
|
|
|
|
import numpy as np
|
|
from scipy.integrate import odeint
|
|
from matplotlib import pyplot as plt
|
|
import pandas as pd
|
|
|
|
#phi_dot = lambda phi, t, I_DC, I_RF, I_J, omega_RF, hbar, e, R: 2*e*R/hbar*( I_DC + I_RF*np.cos(omega_RF*t) - I_J(np.sin(phi)) )
|
|
def phi_dot(phi, t, I_DC, I_RF):
|
|
R = 10.e-3 #Ohm
|
|
I_J = 1.e-3 #A
|
|
omega_RF = 2*np.pi*.96e9 #rad/s
|
|
hbar = 1.0545718e-34 #m^2kg/s
|
|
e = 1.60217662e-19 #C
|
|
|
|
return 2*e*R/hbar*( I_DC + I_RF*np.cos(omega_RF*t) - I_J*(np.sin(phi)) )
|
|
#return R*( I_DC + I_RF*np.sin(omega_RF*t) - I_J*(np.sin(phi)) )
|
|
|
|
# We need an initial value to phi
|
|
phi_0 = 0
|
|
|
|
# Let's try it for 10 periods
|
|
N_points = 1000
|
|
t = np.linspace(0, 628, N_points)
|
|
|
|
hbar = 1.0545718e-34 #m^2kg/s
|
|
e = 1.60217662e-19 #C
|
|
|
|
df = pd.DataFrame(columns=['I_DC','I_RF','V_DC_bar'])
|
|
|
|
#phi = odeint(phi_dot, phi_0, t, (.5e-3, .5e-3))[:, 0]
|
|
#for I_DC in [1e-4, .5e-3, 1.e-3, 1.5e-3, 2.e-3, 2.5e-3]:
|
|
for I_DC in np.arange(0, 1e-3, 1e-5):
|
|
for I_RF in [0., .5e-3, 2.e-3]:
|
|
phi = odeint(phi_dot, phi_0, t, (I_DC, I_RF))
|
|
N_asymp = N_points//2
|
|
I_DC_bar = np.mean(phi[N_asymp:]/t[N_asymp:])
|
|
#V_DC_bar = hbar/(2*e)*I_DC_bar
|
|
V_DC_bar = I_DC_bar
|
|
print("For I_DC =", I_DC, "\t I_RF = ", I_RF, "\twe find V_DC_bar =", V_DC_bar)
|
|
df = df.append({'I_DC': I_DC, 'I_RF': I_RF, 'V_DC_bar': V_DC_bar}, ignore_index = True)
|
|
|
|
## Plotting the thing
|
|
plt.figure()
|
|
|
|
plt.xlabel("$\\overline{V_{DC}}$")
|
|
plt.ylabel("$I_{DC}$")
|
|
|
|
for I_RF in df.I_RF.unique():
|
|
x, y = df[df.I_RF == I_RF][["V_DC_bar", "I_DC"]].to_numpy().T
|
|
plt.plot(x[10:], y[10:], label="$I_{RF} = " + str(I_RF) + "$")
|
|
|
|
plt.legend()
|
|
plt.show() |