#!/bin/env python3

import numpy as np
from scipy.integrate import odeint
from matplotlib import pyplot as plt
import pandas as pd

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)) )
    # I attempted to solve it without the constants, as I suspected overflows
    # were occurring. The next line did not improve the result.
    #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 a lot of periods
N_points = 10000
t = np.linspace(0, 100, 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'])

# For testing:
#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]:
        # The individual solutions for phi do seem sane, at least, the ones
        # I inspected.
        phi = odeint(phi_dot, phi_0, t, (I_DC, I_RF))
        # I initially thought to average over the tail to look at the asymptotic behaviour.
        #N_asymp = N_points//2
        #V_DC_bar = hbar/(2*e)*np.mean(phi[N_asymp:]/t[N_asymp:])
        # Then I choose to just take the last point to see if that gave better results.
        V_DC_bar = hbar/(2*e)*phi[-1]/t[-1]
        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, y, label="$I_{RF} = " + str(I_RF) + "$")

plt.legend()
plt.show()