import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integr

grandeurs=["L","C","RP","RN","R","Usat","V0","derivV0","Tmax"]
unites=["H","F","Ω","Ω","Ω","V","V","V/s","s"]

simu1=[0.04,1e-7,1050,708,1416,15,0.01,0.1,0.005]
simu2=[0.04,1e-7,708,708,1416,15,0.01,0.1,0.005]
simu3=[0.04,1e-7,750,708,1416,15,0.1,0.1,0.02]

def Euler(ED,Y0,T):
    Y=np.zeros((len(T),2))
    Y[0]=Y0
    for i in range(len(T)-1):
        Y[i+1]=Y[i]+(T[i+1]-T[i])*ED(Y[i],T[i])
    return Y

def simulation(valeurs,nombrePas=2000):
    print("***Paramètres imposés***")
    for i in range(len(valeurs)):
        print(grandeurs[i],"=",valeurs[i],unites[i])
    print("***Paramètres calculés***")
    L,C,RP,RN,R,Usat,V0,derivV0,Tmax=valeurs
    S=RN/(R+RN)*Usat
    omega=1/np.sqrt(L*C)
    m=np.sqrt(L/C)*(1/RP-1/RN)/2
    rho=np.sqrt(L/C)*(1/RP+1/R)/2
    f=omega/2/np.pi
    print("S=",S,"V")
    print("m=",m)
    print("omega=",omega,"rad/s")
    print("rho=",rho)
    print("f=",f,"Hz")
    T=np.linspace(0,Tmax,nombrePas)
    def ED(vectV,t):
        if np.abs(vectV[0])>S:
            coeff=rho
        else:
            coeff=m
        rep=[vectV[1],-omega**2*vectV[0]-2*omega*coeff*vectV[1]]
        return np.array(rep)
    Y0=[V0,derivV0]
    #choisir la méthode voulue (mettre # pour désactiver l'autre)
    rep=integr.odeint(ED,Y0, T)
    #rep=Euler(ED,Y0,T)
    plt.figure("X(t)")
    plt.plot(T,rep[:,0])
    plt.figure("V(X)")
    plt.plot(rep[:,0],rep[:,1])
    plt.show()