"""
CHARGE D'UN CIRCUIT RC PAR LA METHODE D'EULER

"""

## Importation des bibliothèques utiles
## ------------------------------------

import numpy as np                  # pour la manipulation des tableaux
import matplotlib.pyplot as plt     # pour les représentations graphiques


## Implémentation de la méthode d'Euler explicite
## ----------------------------------------------

def euler(F, y0, t0, tf, dt):
  """ 
  Résout le problème de Cauchy y'(t)=F(y(t),t) avec y(0)=y0 par la méthode 
  d'Euler explicite.
  Arguments d'entrée : 
    - F : fonction donnant y' (fonction de 2 variables y et t) ;
    - y0 : condition initiale sur y (flottant) ;
    - t0 et tf : bornes de l'intervalle de résolution (flottants) ;
    - dt : pas de discrétisation utilisé pour la résolution (flottant).
  Variables de sortie :
    - t : vecteur contenant l'ensemble des instants tk (array numpy) ;
    - y : vecteur contenant l'ensemble des valeurs approchées yk (array numpy).
  """
  # Initialisation des variables de sortie
  t = np.arange(t0, tf+dt, dt)   # la dernière valeur est dans [tf,tf+dt[
  N = len(t)                     
  y = np.zeros(N)                # initialisation du array y
  y[0] = y0                      # prise en compte de la CI

  # Boucle permettant le calcul des yk par récurrence  
  for k in range(0,N-1):
    y[k+1] = y[k] + F(y[k],t[k])*dt
    
  return t, y # retourne un tableau contenant les temps tk
              # et un tableau contenant les valeurs approchées yk


## Application à la charge d'un circuit RC
## ---------------------------------------

# Définition des constantes utiles
Em = 1      # amplitude du RSF (en volt)
w=1E7       # pulsation du RSF (en rad.s-1)
R = 1000    # en ohm
C = 1E-9    # en farad      
u0 = -2      # condition initiale en volt
tau = R*C   # en seconde

# bornes de l'intervalle de résolution
t0, tf = 0, 10*tau      
t_ex = np.linspace(t0, tf, 200) # vecteur temps


# Définition des fonctions 

def circuitRC_RSF(u, t):
    """
    Fonction explicitant du/dt en fonction de u et t
    pour un circuit RC série soumis à un échelon de tension.
    On a F(u(t),t) = du/dt = -u/tau + (Em/tau) cos(wt)
    """
    return (Em*np.cos(w*t)-u)/tau
        
def u_exacte(t):
    """
    Expression analytique de la solution exacte (forme adimensionnée, correspondant à la
    condition initiale u(0) = u0).
    """ 
    A=Em/(1+(w*tau)**2)
    B=w*tau*A
    return (u0-A)*np.exp(-t/tau)+A*np.cos(w*t)+B*np.sin(w*t)

u_ex = u_exacte(t_ex)             # calcul de la solution exacte

## Résolution  - Représentation graphique
## ---------------------------------------

plt.figure(figsize=(14,8))          # création d'une fenêtre graphique

plt.plot(1e6*t_ex, u_ex, '-', label = "Solution exacte")   # tracé de la solution exacte
ymin=ymax=u0
for N in [1, 0.5, 0.2, 0.1, 0.01]:           # pour différentes valeurs du pas
  dt=N*tau
  t, u = euler(circuitRC_RSF, u0, t0, tf, dt)   # résolution par la méthode d'Euler
  ymin=min(ymin,min(u))
  ymax=max(ymax,max(u))
  plt.plot(1e6*t, u, '--', label = r"Solution correspondant au pas $\delta t$ = "+str(N)+"τ")    
                                    # puis représentation graphique
plt.ticklabel_format(style='plain')
#plt.title("Simulation par la méthode d'Euler explicite, pour $e(t)=E_m \\cos(\\omega t)$ avec :\n$E$ = "+str(Em)+" V ; $\\omega$ = "+str(w)+" rad $s^{-1}$ ; $u(t_0)$ = "+str(u0)+\
#          " V ; R = "+str(R)+" $\Omega$ ; C = "+str(C)+" F")
plt.title(f"Simulation par la méthode d'Euler explicite, pour $e(t)=E_m \\cos(\\omega t)$ avec \
          :\n$E_m$ = {Em} V ; $\\omega$ = {w:.1e} rad $s^{-1}$ ; $u({t0})$ = {u0} V ; R = {R} $\Omega$ ; C = {C} F")

       
plt.xlim(1e6*t0,1e6*tf), plt.xlabel(r"$t$ ($\mu s$)")   # habillage de l'axe des abscisses
plt.ylim(1.05*ymin,1.05*ymax), plt.ylabel(r"$u$ (V)")     # habillage de l'axe des ordonnées
plt.legend(loc = 'lower right')     # affichage de la légende
plt.grid()

plt.show()

## Comparaison avec odeint
## -------------------------------------

from scipy.integrate import odeint

rdt=0.2     # rapport de dt sur tau
dt=rdt*tau  # pas de temps utilisé pour la comparaison


#Résolution par Euler explicite
te,ue=euler(circuitRC_RSF, u0, t0, tf, dt)

#Résolution par odeint
to=[t0+i*dt for i in range(int((tf-t0)/dt))]
uo=odeint(circuitRC_RSF,u0,te)

plt.figure(figsize=(14,8))          # création d'une fenêtre graphique
plt.plot(1e6*t_ex, u_ex, '-', label = "Solution exacte")   
plt.plot(1e6*te,ue,'--',label="Solution par Euler explicite")
plt.plot(1e6*te,uo,'-.',label="Solution par odeint")

plt.ticklabel_format(style='plain')
plt.title(f"Simulation par la méthode d'Euler explicite, pour $e(t)=E_m \\cos(\\omega t)$  \
          \n avec un pas de temps d'intégration dt = {rdt} $\\tau$")

ymin,ymax=min(min(u_ex),min(ue),min(uo)),max(max(u_ex),max(ue),max(uo)) 
      
plt.xlim(1e6*t0,1e6*tf), plt.xlabel(r"$t$ ($\mu s$)")   # habillage de l'axe des abscisses
plt.ylim(1.05*ymin,1.05*ymax), plt.ylabel(r"$u$ (V)")     # habillage de l'axe des ordonnées
plt.legend(loc = 'lower right')     # affichage de la légende
plt.grid()

plt.show()