import numpy as np
import matplotlib.pyplot as plt

##
f = np.array([...])
UR = np.array([...])
deltaT = np.array([...])

##

R = ...
C = ...
L = ...
r = ...
r_L =  ...

##
omega = 2*np.pi*f
I = UR/R
phi = -omega*deltaT
##
NbPoints = 10000
omega0 = 1/np.sqrt(L*C)
Q = 1/(R+r+r_L)*np.sqrt(L/C)
omega_th = np.linspace(np.min(omega), np.max(omega), NbPoints)
I_th = 1/(1+1j*Q*(omega_th/omega0 - omega0/omega_th))*2/(R+r+r_L)
##

plt.close('all')
plt.figure('Résonance en intensité', figsize = [16/2.54, 20/2.54])
plt.axes([0.18, 0.55, 0.8, 0.44])
plt.plot(omega_th, np.abs(I_th), 'k')
plt.plot(omega, I, 'rs', markeredgecolor = 'k', ms = 6 )
plt.ylabel(r'$I \,{(\rm A)}$', fontsize = 12)
plt.xlim([np.min(omega_th), np.max(omega_th)])

plt.axes([0.18, 0.1, 0.8, 0.44])
plt.plot(omega_th, np.angle(I_th), 'k')
plt.plot(omega, phi, 'rs', markeredgecolor = 'k', ms = 6 )
plt.ylabel(r'$\arg(\varphi)\,{(\rm rad)}$', fontsize = 12)
plt.xlabel(r'$\omega\,{\rm (rad.s^{-1})}$')
plt.xlim([np.min(omega_th), np.max(omega_th)])
plt.show()