# -*- coding: utf-8 -*-

# Programme de tracé de la courbe u = f(w/w0)
# Dominique Lefebvre pour TangenteX.com
# 25 janvier 2015
#

# importation des librairies
from numpy import array,arange,pi,cos,zeros
from scipy.integrate import odeint 
import matplotlib.pyplot as plt
     
# définition des constantes du circuit RLC
Omega0 = 5.  # Pulsation propre du circuit  

# définition des paramètres du forçage
A = 1.       # amplitude
phi = 0.     # phase du forçage


# fonction de définition du circuit RLC
def RLC(y,t):
    u,u_dot = y
    u_dot_dot = -K*u_dot - Omega0**2*u + A*cos(OmegaF*t + phi) 
    return [u_dot,u_dot_dot]   

# saisie du coefficent d'amortissement
K = float(input("Coefficient d'amortissement : "))

# conditions initiales
u0 = 1      
du0 = 0
C0 = array([u0,du0])

# définition du vecteur temps de l'expérience
t0 = 0.0
tmax = 10*pi
pastemps = 0.01
time = arange(t0, tmax, pastemps)

# définition du vecteur de variation de w/w0
vrmin = 0.1
vrmax = 2.
vrpas = 0.01
nbpas = (2.0 - 0.1)/vrpas

# calcul de la courbe
cvr = zeros(nbpas)
vr = zeros(nbpas)
i = vrmin
j = 0
while i < vrmax:
    OmegaF = Omega0*i
    y = odeint(RLC,C0,time)
    vr[j] = i
    cvr[j] = y[:,0].max()
    j = j + 1
    i = i + vrpas
        
# affichage de l'évolution de la tension aux bornes de C
plt.figure()
plt.xlim(vr.min(),vr.max())
plt.ylim(cvr.min(),cvr.max())
plt.plot(vr, cvr)
plt.title('Resonance - Amortissement = ' + str(K))
plt.xlabel('w/w0')
plt.ylabel('Amplitude')
plt.show()