# -*- coding: utf-8 -*-
# Created on Tue Mar 28 20:05:00 2023
# @author: mhebding

# Bibliothèques

import numpy as np
import matplotlib.pyplot as plt

## Sujet 2 ##
# Attracteur de Lorentz

def euler(ED, t0, tf, y0, n):
    '''
    Parameters
    ----------
    ED : une fonction de l equa diff de y
    t0 : instant initial
    tf : instant final
    y0 : condition initiale de y
    n : un nombre de points

    Returns
    -------
    liste_t : liste des instants
    liste_y : liste des y(t)

    '''
    h = (tf - t0) / n  # le pas ∆t noté h
    y, t = y0, t0 # initialisation des conditions initiales
    liste_y = [y0]  # liste des y(t)
    liste_t = [t0]  # liste des t
    for k in range(n):  # n iterations
        y = y + h * ED(y, t) # COEUR DE LA METHODE
        t = t + h  # instant suivant
        liste_t.append(t)  # on ajoute la valeur de t dans la liste
        liste_y.append(y)  # on ajoute la valeur de y dans la liste
    return liste_t, liste_y  # on renvoie les deux listes

def ED_Lorentz(M, t):
    x, y, z = M
    dM = np.array([
    10*(y-x),
    28*x-y-x*z,
    x*y-8/3*z
    ])
    return dM # renvoie la matrice dM = dx, dy, dz (voir l equa diff)

M0 = np.array([1, 1, 1])  # conditions initiales
t, M = euler(ED_Lorentz, 0, 50, M0, 10000) # resolution
M = np.array(M) # pour assurer le type np.array
x, y, z = M[:, 0], M[:, 1], M[:, 2]

# Graphes
#plt.plot(t,x)
#plt.plot(t,y)
#plt.plot(t,z)

from mpl_toolkits import mplot3d
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot3D(x, y, z)
plt.show()