import numpy as np
import matplotlib.pyplot as plt
from pylab import *
from math import *
from scipy.integrate import odeint

## I- Script 1 : force centrale K=-1 et n=-1

# Définition de l'équation différentielle
def système (X,t):  # Attention à l'inversion de t et X
    Vx = X[0]
    Vy = X[1]
    x = X[2]
    y = X[3]
    X_point = [K*x*(x*x+y*y)**(-(n+1)/2),K*y*(x*x+y*y)**(-(n+1)/2),Vx,Vy]
    return X_point

# Données du système
K = -1
n = -1

# Résolution
t = np.linspace(0,10,1000)
X = odeint(système, [0,1,1,0], t)
Vx = X[:,0]
Vy = X[:,1]
x = X[:,2]
y = X[:,3]
print(Vx) # renvoie les valeurs successives de Vx
print(Vy) # renvoie les valeurs successives de Vy
print(x) # renvoie les valeurs successives de x
print(y) # renvoie les valeurs successives de y

# Représentations graphiques
plt.figure(1)
plt.plot(x, y)
plt.title('Représentation de la trajectoire y en fonction de x')
plt.ylabel('y')
plt.xlabel('x')
plt.grid()
plt.show()


## II- Script 2 : problème de Kepler, à la 1ère vitesse cosmique

# Définition de l'équation différentielle
def système (X,t):  # Attention à l'inversion de t et X
    r = X[0]
    rpoint = X[1]
    theta = X[2]
    thetapoint = X[3]
    X_point = [rpoint,r*thetapoint**2-G*Mt/(r**2),thetapoint,-2*rpoint*thetapoint/r]
    return X_point

# Données du système
G = 6.67*10**(-11)
Mt = 6*10**(24)
Rt = 6400000

# Résolution
t = np.linspace(0,100000,100000)
X = odeint(système, [Rt,0,0,(1/Rt)*np.sqrt(1*G*Mt/Rt)], t)
r = X[:,0]
x=r*np.cos(X[:,2])
y=r*np.sin(X[:,2])

# Représentations graphiques
plt.figure(2)
plt.plot(x,y)
plt.title('Représentation de la trajectoire y en fonction de x')
plt.ylabel('y en m')
plt.xlabel('x en m')
axis("equal")
plt.grid()
plt.show()


## III- Script 3 : problème de Kepler, v supérieure à la 1ère vitesse cosmique

# Définition de l'équation différentielle
def système (X,t):  # Attention à l'inversion de t et X
    r = X[0]
    rpoint = X[1]
    theta = X[2]
    thetapoint = X[3]
    X_point = [rpoint,r*thetapoint**2-G*Mt/(r**2),thetapoint,-2*rpoint*thetapoint/r]
    return X_point

# Données du système
G = 6.67*10**(-11)
Mt = 6*10**(24)
Rt = 6400000

# Résolution
t = np.linspace(0,100000,100000)
X = odeint(système, [Rt,0,0,(1/Rt)*np.sqrt(1.5*G*Mt/Rt)], t)
r = X[:,0]
x=r*np.cos(X[:,2])
y=r*np.sin(X[:,2])

# Représentations graphiques
plt.figure(2)
plt.plot(x,y)
plt.title('Représentation de la trajectoire y en fonction de x')
plt.ylabel('y en m')
plt.xlabel('x en m')
axis("equal")
plt.grid()
plt.show()

## IV- Script 4 : problème de Kepler, v égale à la 2nde vitesse cosmique

# Définition de l'équation différentielle
def système (X,t):  # Attention à l'inversion de t et X
    r = X[0]
    rpoint = X[1]
    theta = X[2]
    thetapoint = X[3]
    X_point = [rpoint,r*thetapoint**2-G*Mt/(r**2),thetapoint,-2*rpoint*thetapoint/r]
    return X_point

# Données du système
G = 6.67*10**(-11)
Mt = 6*10**(24)
Rt = 6400000

# Résolution
t = np.linspace(0,100000,100000)
X = odeint(système, [Rt,0,0,(1/Rt)*np.sqrt(2*G*Mt/Rt)], t)
r = X[:,0]
x=r*np.cos(X[:,2])
y=r*np.sin(X[:,2])

# Représentations graphiques
plt.figure(2)
plt.plot(x,y)
plt.title('Représentation de la trajectoire y en fonction de x')
plt.ylabel('y en m')
plt.xlabel('x en m')
axis("equal")
plt.grid()
plt.show()

## V- Script 5 : problème de Kepler, v supérieure à la 2nde vitesse cosmique

# Définition de l'équation différentielle
def système (X,t):  # Attention à l'inversion de t et X
    r = X[0]
    rpoint = X[1]
    theta = X[2]
    thetapoint = X[3]
    X_point = [rpoint,r*thetapoint**2-G*Mt/(r**2),thetapoint,-2*rpoint*thetapoint/r]
    return X_point

# Données du système
G = 6.67*10**(-11)
Mt = 6*10**(24)
Rt = 6400000

# Résolution
t = np.linspace(0,100000,100000)
X = odeint(système, [Rt,0,0,(1/Rt)*np.sqrt(3*G*Mt/Rt)], t)
r = X[:,0]
x=r*np.cos(X[:,2])
y=r*np.sin(X[:,2])

# Représentations graphiques
plt.figure(2)
plt.plot(x,y)
plt.title('Représentation de la trajectoire y en fonction de x')
plt.ylabel('y en m')
plt.xlabel('x en m')
axis("equal")
plt.grid()
plt.show()