# %% Imports & Constants
import numpy as np
import matplotlib.pyplot as plt


## Intégration par la méthode d'Euler explicite - formulation matricielle
def euler_explicite(F, Y0, t0, tf, n):
    """
    Résout numériquement le système d'équations différentielles Y' = F(t, Y) avec la méthode d'Euler explicite.

    Paramètres:
    F : fonction
        La fonction définissant le système d'équations différentielles.
    Y0 : array_like
        La condition initiale Y(t0).
    t0 : float
        Le temps initial.
    tf : float
        Le temps final.
    n : int
        Le nombre de pas de temps.

    Retourne:
    t_values : array_like
        Les valeurs de temps.
    Y_values : array_like
        Les valeurs approchées de Y aux instants correspondants.
    """
    dt = (tf - t0) / n  # Taille du pas de temps
    t_values = np.linspace(t0, tf, n + 1)
    m = len(Y0)  # Dimension du système
    Y_values = np.zeros((n + 1, m))
    Y_values[0, :] = Y0

    for i in range(n):
        Y_values[i + 1, :] = Y_values[i, :] + dt * F(t_values[i], Y_values[i, :])

    return t_values, Y_values


# %% Application: résolution du pendule simple avec frottement
if __name__ == "__main__":
    # Paramètres du pendule
    g = 9.81  # Accélération gravitationnelle (m/s^2)
    R = 0.1  # Longueur du pendule (m)
    m = 0.5  # Masse du pendule (kg)
    f = 1.0  # Coefficient de frottement (N·m·s)

    # Fonction définissant le système d'équations différentielles
    def F(t, Y):
        # à compléter
        return np.array([dtheta_dt, domega_dt])

    # Paramètres de la simulation
    t0 = 0.0  # Temps initial (s)
    tf = 8.0  # Temps final (s)
    n = 1000  # Nombre de pas de temps

    # Conditions initiales
    Theta0 = [0, np.pi]  # Angle initial (radians)
    Omega0 = [25, 0]  # Vitesse angulaire initiale (rad/s)

    for theta0, omega0 in zip(Theta0, Omega0):
        Y0 = np.array([theta0, omega0])

        # Résolution du système d'équations différentielles
        t_values, Y_values = euler_explicite()  # arguments à compléter

        # Extraction des résultats
        theta_values = Y_values[:, 0]
        omega_values = Y_values[:, 1]

        # Affichage des résultats
        plt.figure()
        plt.plot(
            t_values,
            theta_values,
            label=rf"Angle $\theta(t)$, $\theta_0={theta0/np.pi:.1f}\times\pi$ rad",
        )
        plt.plot(
            t_values,
            omega_values,
            label=rf"Vitesse angulaire $\omega(t)$, $\omega_0={omega0:.1f}$ rad/s",
        )
        plt.xlabel("Temps (s)")
        plt.title("Résolution du pendule simple avec Euler explicite")
        plt.xlim(t0, tf)
        plt.legend()
        plt.grid()
        plt.show()

# %%