from matplotlib import pyplot as plt
import numpy as np
import sys

def pausa():
    return input("Premi ENTER per continuare...")

# periodo di rivoluzione in unità dell'anno terrestre 
T = np.array([0.241, 0.615, 1, 1.88, 11.9, 29.4, 83.7, 163.7])
# 'distanza' dal Sole (-> semiasse maggiore) in funzione dell'U.A.
d = np.array([0.387, 0.723,  1 , 1.52, 5.20, 9.57, 19.17, 30.18])
n = len(T)

# limitiamo il nr di pianeti da analizzare (min 2, max 8)
risp = input("Numero di pianeti (2-8; default 8): ")
nl = 8 
if len(risp) >0 :
    nl = int(risp)

if  1 < nl <= 8:
    n = nl
    T = T[:n]
    d = d[:n]
    print(" => Consideriamo solo i primi %d pianeti"%n)
else:
    print(" %d non è valido"%n1)
    sys.exit(1)


# primo plot, d vs T
figura, ax = plt.subplots()
ln1, = ax.plot(T, d, 'o', color="blue")
ax.grid()

ax.set_title("Pianeti del sistema solare: 'distanza' vs periodo")
ax.set_xlabel("T (y)")
ax.set_ylabel("d (R_T)")
figura.show()

pausa()

# linearizzazione, confidando nella Terza Legge di Keplero
figura, ax = plt.subplots()

ln1, = ax.plot(T**2, d**3, 'o', color="blue")
ax.grid()

ax.set_title("Pianeti del sistema solare: -- linearizzazione mediante Keplero")
ax.set_xlabel("T^2 (y^2)")
ax.set_ylabel("d^3 (R_T^3)")

ascisse = np.linspace(0, max(T*T))
linea = ascisse # coeff angolare = 1
ax.plot(ascisse, linea, color="cyan")

figura.show()

pausa()

# linearizzazione mediante scale logaritmiche ('doppio-log', o 'log-log')
figura, ax = plt.subplots()

ln1, = ax.plot(T, d, 'o', color="blue")
ax.grid()

ax.set_xscale('log')
ax.set_yscale('log')

ax.set_title("Pianeti del sistema solare -- scale logaritmiche")
ax.set_xlabel("T^2 (y^2)")
ax.set_ylabel("d^3 (R_T^3)")

# sovrapponiamo quanto atteso dalla Terza Legge di Keplero
ascisse = np.linspace(0.5*T[0], 1.5*T[-1], 100)
ax.plot(ascisse, ascisse**(2/3), linestyle="dashed", linewidth=1.5, color="cyan")
figura.show()

print(" -> la linea dà l'andamento secondo Keplero, ovvero d prop. a T^(3/2)") 

pausa()

# aggiorniamo la figura esistente...
ax.plot(T[0], d[0], 'o',  fillstyle="none", color="magenta", markersize = 16)
ax.plot(T[-1], d[-1], 'o', fillstyle="none", color="magenta", markersize = 16)

m1 = ( np.log(d[-1]) - np.log(d[0]) ) / ( np.log(T[-1]) - np.log(T[0]) )
print("Potenza empirica dell'andamento di d in funzione di T: %.4f"%m1)
print(" --> Ricavata facendo uso del primo e dell'ultimo punto")
print("     (inutile fare la retta che passa meglio fra i punti, dati gli scarti esigui)")

ascisse = np.linspace(0.5*T[0], 1.5*T[-1], 100)
ax.plot(ascisse, ascisse**m1, linestyle="dotted", linewidth=1.5, color="magenta")

figura.canvas.draw()
input("Premi return per uscire")

# volendo, si può anche trovare la slope della migliore
# retta fra i punti, ma gli scostamenti dei punti dalle
#  due rette trovate sono talmente esigui che non ne vale la pena