Visualizzazione di orbite di vario tipo di un corpo di massa minore attorno ad uno di massa maggiore

from matplotlib import pyplot as grafico
import numpy as np
import math
from scipy import special

print("Benvenuto nel calcolatore di orbite!")
tipo_corpo = input("Specificare se il corpo minore è un satellite: ")
nome = input("Inserire il nome del corpo minore: ")
NOMONE = input("Inserire il nome del corpo maggiore: ")
stella = input("Specificare se il corpo maggiore è una stella: ")
M = float(input("Inserire la massa del corpo maggiore in kg (includere decimali): "))
if stella.lower() == 'stella':
    Lum = float(input("Inserire luminosità solare (includere decimali): "))
    albedo_domanda = input("L'albedo del corpo minore è noto? ")
    if albedo_domanda.lower() == "sì" or albedo_domanda.lower() == "si":
        albedo = float(input("Inserire l'albedo del corpo minore, preferibilmente Bond (includere decimali): "))
R_s = float(input("Inserire il raggio del corpo maggiore in m (includere decimali): "))
R_p = float(input("Inserire il raggio del corpo minore in m (includere decimali): "))
m = float(input("Inserire la massa del corpo minore in kg (includere decimali): "))
T = float(input("Inserire il tempo di rivoluzione del corpo minore in s (includere decimali, in secondi): "))
e = float(input("Inserire l'eccentricità dell'orbita (includere decimali): "))
G = 6.67430e-11
AU = 1.495978707e11
sigma = 5.670374419e-8 #costante di Stefan-Boltzmann

a = np.cbrt(((T**2)*G*M)/(4*(np.pi**2)))
#E = (-G*M*m)/(2*a)
#L = m*(2*np.pi/T)*(a**2)


#calcoli
if stella.lower() == 'stella':
    Rmin = np.sqrt(Lum/1.1)*AU
    Rmax = np.sqrt(Lum/0.53)*AU
    Lim_Roche = R_s*2.44*((M/(4/3*np.pi*R_s**3))/(m/(4/3*np.pi*R_p**3)))**(1/3)
b = a*np.sqrt(1-e**2)
if tipo_corpo.lower() == 'satellite':
    rP = a*(1-e) - R_s
    rA = a*(1+e) - R_s
    v_h1 = np.sqrt((G*M)/(R_s+rP))
    v_h2 = np.sqrt((G*M)/(R_s+rA))
if tipo_corpo.lower() != 'satellite':
    rP = a*(1-e)
    rA = a*(1+e)
    Lim_Roche = R_s*2.44*((M/(4/3*np.pi*R_s**3))/(m/(4/3*np.pi*R_p**3)))**(1/3)
if stella.lower() == 'stella':
    if albedo_domanda.lower() == "sì" or albedo_domanda.lower() == "si":
        L_stella = Lum*3.828e26
        T_p = ((L_stella*(1-albedo))/(16*np.pi*sigma*(rP)**2))**(1/4)
        T_a = ((L_stella*(1-albedo))/(16*np.pi*sigma*(rA)**2))**(1/4)
dF = e*a
rA_AU = rA/AU
rP_AU = rP/AU
v_p = np.sqrt((G*M*(1+e))/(a*(1-e)))
v_a = np.sqrt((G*M*(1-e))/(a*(1+e)))

#Costruiamo l'orbita
angolo = np.linspace(0, 2*np.pi)
x = a*np.cos(angolo)
y = b*np.sin(angolo)
if stella.lower() == 'stella':
    Rmin_x = Rmin*np.cos(angolo) #abitabile min
    Rmin_y = Rmin*np.sin(angolo)
    Rmax_x = Rmax*np.cos(angolo) #abitabile max
    Rmax_y = Rmax*np.sin(angolo)
    #concatenazione per fare percorso chiuso
    x_R = np.concatenate([Rmin_x, Rmax_x[::-1]]) #parte quindi dai valori di Rmax, finisce in quelli di Rmin e inizia però dall'inverso di questi ultimi
    y_R = np.concatenate([Rmin_y, Rmax_y[::-1]])
if tipo_corpo.lower() != "satellite":
    Roche_x = Lim_Roche*np.cos(angolo)
    Roche_y = Lim_Roche*np.sin(angolo)
    Area_occupata_x = R_s*np.cos(angolo)
    Area_occupata_y = R_s*np.sin(angolo)
if tipo_corpo.lower() == 'satellite':
    afelio = rA
    perielio = -rP
    xA = rA*np.cos(angolo)
    xP = rP*np.sin(angolo)
else:
    afelio = rA-dF
    perielio = -rP-dF

grafico.grid()
if tipo_corpo.lower() == 'satellite':
    orbita = grafico.plot(xA, xP, label="Orbita di " + nome + " attorno a " + NOMONE + ".", color = "blue")
else:
    orbita = grafico.plot(x, y, label="Orbita di " + nome + " attorno a " + NOMONE + ".", color = "blue")
posizione_M = grafico.scatter(-dF, 0, label="Posizione di " + NOMONE, s=50, color="yellow", zorder=5) #zorder è per la SOVRAPPOSIZIONE: più alto il valore, più sovrapposto a qualcos'altro
if tipo_corpo.lower() == 'satellite':
    grafico.scatter(afelio, 0, color="red", zorder=5, label="Distanza maggiore dalla superficie: " + str(rA/1000) + " km; velocità = " + str(v_h2*3.6) + " km/h")
    grafico.scatter(perielio, 0, color="purple", zorder=5, label="Distanza minore dalla superficie: " + str(rP/1000) + " km; velocità = " + str(v_h1*3.6) + " km/h")
else:
    if stella.lower() == 'stella':
        if albedo_domanda.lower() == "sì" or albedo_domanda.lower() == "si":
            grafico.scatter(afelio, 0, color="red", zorder=5, label="Distanza maggiore: " + str(rA/1000) + " km = " + str(rA_AU) + " AU; velocità = " + str(v_a*3.6) + " km/h. Temperatura (senza effetto serra): " + str(T_a) + " K = " + str(T_a-273.15) + " °C")
            grafico.scatter(perielio, 0, color="purple", zorder=5, label="Distanza minore: " + str(rP/1000) + " km = " + str(rP_AU) + " AU; velocità = " + str(v_p*3.6) + " km/h. Temperatura (senza effetto serra): " + str(T_p) + " K = " + str(T_p-273.15) + " °C")
        else:
            grafico.scatter(afelio, 0, color="red", zorder=5, label="Distanza maggiore: " + str(rA/1000) + " km = " + str(rA_AU) + " AU; velocità = " + str(v_a*3.6) + " km/h")
            grafico.scatter(perielio, 0, color="purple", zorder=5, label="Distanza minore: " + str(rP/1000) + " km = " + str(rP_AU) + " AU; velocità = " + str(v_p*3.6) + " km/h")
    else:
        grafico.scatter(afelio, 0, color="red", zorder=5, label="Distanza maggiore: " + str(rA/1000) + " km = " + str(rA_AU) + " AU; velocità = " + str(v_a*3.6) + " km/h")
        grafico.scatter(perielio, 0, color="purple", zorder=5, label="Distanza minore: " + str(rP/1000) + " km = " + str(rP_AU) + " AU; velocità = " + str(v_p*3.6) + " km/h")
if stella.lower() == 'stella':
    raggio_int = grafico.plot(Rmin_x-dF, Rmin_y, label="Raggio minimo della zona abitabile di " + NOMONE + " (" + str(Rmin/AU) + " AU)", color = "black")
    raggio_ext = grafico.plot(Rmax_x-dF, Rmax_y, label="Raggio massimo della zona abitabile di " + NOMONE + " (" + str(Rmax/AU) + " AU)", color = "black")
    area = grafico.fill(x_R-dF, y_R, label="Area della zona abitabile di " + NOMONE, color = "lightgreen", alpha=0.4) #alpha rende più opaco
    Roche = grafico.fill(Roche_x-dF, Roche_y, label="Area del limite di Roche di " + NOMONE + " nel caso di " + nome + ". Limite di Roche: " + str(Lim_Roche/1000) + " km = " + str(Lim_Roche/AU) + " AU", color = "red", alpha=0.4)
    area_stella = grafico.fill(Area_occupata_x-dF, Area_occupata_y, label="Area occupata da " + NOMONE + ". Raggio: " + str(R_s/1000) + " km = " + str(R_s/AU) + " AU", color="yellow", alpha = 0.3)
if tipo_corpo.lower() != "satellite" and stella.lower() != "stella":
    Roche_pianeta = grafico.fill(Roche_x-dF, Roche_y, label="Area del limite di Roche di " + NOMONE + " nel caso di " + nome + ". Limite di Roche: " + str(Lim_Roche/1000) + " km = " + str(Lim_Roche/AU) + " AU", color = "red", alpha=0.4)
    area_pianeta = grafico.fill(Area_occupata_x-dF, Area_occupata_y, label="Area occupata da " + NOMONE + ". Raggio: " + str(R_s/1000) + " km = " + str(R_s/AU) + " AU", color="yellow", alpha = 0.3)
grafico.axhline(y=0, color="grey", linestyle="solid")
grafico.axvline(x=0, color="grey", linestyle="solid")
grafico.axis("equal") #EVITA CHE I CERCHI SI DISTORGANO, centra gli assi
grafico.title("Orbita ellittica di un corpo attorno a un altro")
grafico.legend(loc="upper right", fontsize=7.5) #Mostra i nomi delle funzioni presentate e i colori, in alto a destra
grafico.xlabel("Raggio (m)")
grafico.ylabel("Raggio (m)")
grafico.show()