# -*- coding: utf-8 -*-

"""Trajectoire d'un point soumis à un champ de force centrale conservative."""

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button, RadioButtons
from matplotlib import rc


m = 7.36e22  # kg masse de la Lune
G = 6.67e-11 # m^3.kg^-1.s^-2
MT = 6e24# kg masse de la Terre
# Conditions initiales
theta0 = 0  # rad
dtheta0 = 2*np.pi/(27.3*86400)  # rad/s
r0 = 384e6  # m
dr0 = 0  # m/s



def Ep_eff(r, dtheta0) : 
    C = r0**2*dtheta0
    return 0.5*m*C**2/r**2 - G*MT*m/r
     



def Fc(r: float) -> float:
    """Force centrale."""
    # Constante de gravitation
    G = 6.67e-11  # N m²/kg²
    # Masse de l'astre central
    MA = 5.97e24  # kg
    return -G*MA*m/r**2


# Durée de la simulation
tmax = 10*30*86400  # s
# Nombre de points
N = 1000
# Liste des instants
t = np.linspace(0, tmax, N)


def deriv(y, t):
    """Système d'équations du premier ordre à résoudre.
    dX représente la dérivée temporelle de X"""
    theta, dtheta, r, dr = y
    dy = [dtheta, -2*dr*dtheta/r, dr, r*dtheta**2 + Fc(r)/m]
    return dy




def slider_changed (event):
    dessine()
    

def dessine():	


    dtheta0  = amplitude1.val

    # C = r0**2*dtheta0

    r = np.linspace(1e8, 10e9, 1000)
    
    # plt.plot(r, Ep_eff(r, dtheta0))
    # plt.ylim([-4e28, 1e28])

    l1.set_xdata(r)
    l1.set_ydata(Ep_eff(r, dtheta0))
    
    l2.set_xdata([1e8, 1e10])
    
    y=[0, 0]
    l2.set_ydata(y)
    Em = 0.5*m*(r0*dtheta0)**2 - G*MT*m/(r0)
    l3.set_xdata([1e8, 1e10])
    
    y=[Em, Em]
    l3.set_ydata(y)
    
    ax.axes.axis([1e8, 1e10, -0.5e29,0.5e29])
    # ax.grid()
    
    sol = odeint(deriv, [theta0, dtheta0, r0, dr0], t)
    theta = sol[:, 0]
    r = sol[:, 2]
    l4.set_xdata(r*np.cos(theta))
    
    l4.set_ydata(r*np.sin(theta))

    l5.set_xdata([0])
    
    l5.set_ydata([0])


    bx.axes.axis([5e8-1e10, 5e8, -0.5e10, 0.5e10])
    ax.set_ylabel(r'$E_{p, {\rm eff}}$ (J)' , fontsize=20)
    ax.set_xlabel(r'$r $ (m)' , fontsize=20)
    
    plt.draw()


fig = plt.figure(figsize=[16,8])

ax = fig.add_axes([0.1,0.18,0.4,0.8])
l1,l2,l3 =plt.plot([],[],'k',[],[],'-r',[],[],'--k')

bx = fig.add_axes([0.52,0.18,0.4,0.8])
l4, l5 = plt.plot([], [], 'r', [],[], 'ko', markersize = 12)

sld_pos        = [0.1, 0.05, 0.8, 0.05]

sld_amplitude1 = plt.axes(sld_pos, facecolor='grey')
amplitude1     = Slider(sld_amplitude1, r'$\dot{\theta}(t=0) rad.s^{-1}$', 2e-6, 5e-6, valinit=2.66e-6)
amplitude1.on_changed(slider_changed)

dessine()

plt.show()