ass3: Draft the script for task 12

ass3-12-a-weak-junction.py

@@ -0,0 +1,38 @@
+#!/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        = 10e-3 #Ohm
+    I_J      = 1e-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, 1000, 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_RF in [0., .5e-3, 2.e-3]:
+        phi = odeint(phi_dot, phi_0, t, (I_DC, I_RF))
+        I_DC_bar = np.mean(phi[N_points//2:]/t[N_points//2:])
+        V_DC_bar = hbar/(2*e)*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)