From a5eea861df05363fcc7dd7ca8fc2082a41bd1a7f Mon Sep 17 00:00:00 2001
From: Kees van Kempen <k.vankempen@student.science.ru.nl>
Date: Thu, 17 Mar 2022 17:12:50 +0100
Subject: [PATCH] ass3: Task 12: does not work, but it is something

---
 ass3-12-a-weak-junction.py | 34 +++++++++++++++++++++++++---------
 1 file changed, 25 insertions(+), 9 deletions(-)

diff --git a/ass3-12-a-weak-junction.py b/ass3-12-a-weak-junction.py
index b767016..6d3dab1 100644
--- a/ass3-12-a-weak-junction.py
+++ b/ass3-12-a-weak-junction.py
@@ -7,21 +7,21 @@ 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
+    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)) )
+    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)
+t = np.linspace(0, 628, N_points)
 
 hbar     = 1.0545718e-34 #m^2kg/s
 e        = 1.60217662e-19 #C
@@ -29,10 +29,26 @@ 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 [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))
-        I_DC_bar = np.mean(phi[N_points//2:]/t[N_points//2:])
-        V_DC_bar = hbar/(2*e)*I_DC_bar
+        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)
\ No newline at end of file
+        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()
\ No newline at end of file