Orbite V2

from matplotlib import pyplot as plt
import numpy as np
import math
import scipy
from scipy import special
import scipy.constants

#interazione gravitazionale tra 2 corpi, usando la formula di Newton
nome_corpo_maggiore = input("Inserisci il nome del corpo di massa maggiore: ")
massa_maggiore = float(input(f"Inserisci la massa di {nome_corpo_maggiore} in kg: "))
raggio_maggiore = float(input(f"Inserisci il raggio di {nome_corpo_maggiore} in m: "))

nome_corpo_minore = input("Inserisci il nome del corpo di massa minore: ")
massa_minore = float(input(f"Inserisci la massa di {nome_corpo_minore} in kg: "))
raggio_minore = float(input(f"Inserisci il raggio di {nome_corpo_minore} in m: "))

T = float(input(f"Inserisci il tempo di rivoluzione di {nome_corpo_minore} in s: "))
e = float(input(f"Inserisci l'eccentricità dell'orbita di {nome_corpo_minore}: "))

while True:
    condizione_stella = input(f"{nome_corpo_maggiore} è una stella? ")
    if condizione_stella.lower() == "sì":
        stella_switch = True
        lum = float(input(f"Inserisci la luminosità di {nome_corpo_maggiore} rispetto al Sole: "))
        lum_switch = True
        condizione_albedo = input(f"L'albedo di {nome_corpo_minore} è noto? ")
        if condizione_albedo.lower() == "sì":
            albedo_switch = True
            albedo = float(input(f"Inserisci l'albedo di {nome_corpo_minore}: "))
            break
        elif condizione_albedo.lower() == "no":
            albedo_switch = False
            break
        else:
            print("Input non valido")
    elif condizione_stella.lower() == "no":
        stella_switch = False
        lum_switch = False
        albedo_switch = False
        break

while True:
    condizione_satellite = input(f"{nome_corpo_minore} è un satellite artificale? ")
    if condizione_satellite.lower() == "sì":
        satellite_switch = True
        break
    elif condizione_satellite.lower() == "no":
        satellite_switch = False
        break

while True:
    condizione_rotazione = input(f"Il periodo di rotazione di {nome_corpo_maggiore} è noto? ")
    if condizione_rotazione.lower() == "sì":
        rot = float(input(f"Inserisci il periodo di rotazione di {nome_corpo_maggiore} in s: "))
        rotazione_switch = True
        break
    elif condizione_rotazione.lower() == "no":
        rotazione_switch = False
        break

#costanti
G = scipy.constants.gravitational_constant
AU = scipy.constants.astronomical_unit
sigma = scipy.constants.Stefan_Boltzmann
ly = scipy.constants.light_year

#calcoli
a = np.cbrt(((T**2)*G*massa_maggiore)/(4*(np.pi**2))) #semiasse maggiore
b = a*np.sqrt(1-e**2)

raggio_perielio = a*(1-e)
raggio_afelio = a*(1+e)
Lim_Roche = raggio_maggiore*2.44*((massa_maggiore/(4/3*np.pi*raggio_maggiore**3))/(massa_minore/(4/3*np.pi*raggio_minore**3)))**(1/3) #satelliti liquidi

velocità_perielio = np.sqrt(((G*massa_maggiore)/a)*((1+e)/(1-e)))
velocità_afelio = np.sqrt(((G*massa_maggiore)/a)*((1-e)/(1+e)))

if stella_switch is True:
    if lum_switch is True:
        #da https://web.archive.org/web/20230314210750/https://www.planetarybiology.com/calculating_habitable_zone.html
        raggio_abitabile_minimo = np.sqrt(lum/1.1)*AU #in AU il calcolo
        raggio_abitabile_massimo = np.sqrt(lum/0.53)*AU #in AU il calcolo
        if albedo_switch is True:
            energia_stellare = lum*3.828e26
            Temp_perielio = ((energia_stellare*(1-albedo))/(16*np.pi*sigma*(raggio_perielio)**2))**(1/4)
            Temp_afelio = ((energia_stellare*(1-albedo))/(16*np.pi*sigma*(raggio_afelio)**2))**(1/4)

if rotazione_switch is True:
    numero_orbite_giorno = rot/T


#calcoli per il plot
distanza_centro_fuoco = e*a
angolo = np.linspace(0, 2*np.pi, 100000)
x_orbita = a*np.cos(angolo) #costruzione orbita sull'asse x
y_orbita = b*np.sin(angolo) #costruzione orbita sull'asse y
distanza_corpomaggiore_centro = -distanza_centro_fuoco #per semplicità nella visualizzazione, si assume quindi che il perielio sia a sx nel plot e l'afelio a dx nel plot

if stella_switch is True:
    abitabile_minimo_x = raggio_abitabile_minimo*np.cos(angolo) #raggio abitabile min
    abitabile_minimo_y = raggio_abitabile_minimo*np.sin(angolo)
    abitabile_massimo_x = raggio_abitabile_massimo*np.cos(angolo) #raggio abitabile max
    abitabile_massimo_y = raggio_abitabile_massimo*np.sin(angolo)
    #concatenazione per fare percorso chiuso
    #si crea un anello, prima con dati in senso antiorario (min_x,y) e poi in senso orario (inverso di max_x,y)
    x_R = np.concatenate([abitabile_minimo_x, abitabile_massimo_x[::-1]])
    y_R = np.concatenate([abitabile_minimo_y, abitabile_massimo_y[::-1]])

if satellite_switch is True:
    raggio_x_corpomaggiore = raggio_maggiore*np.cos(angolo)
    raggio_y_corpomaggiore = raggio_maggiore*np.sin(angolo)

#costruzione orbita
if rotazione_switch is False:
    plt.plot(x_orbita,
            y_orbita,
            label=f"Orbita di {nome_corpo_minore} attorno a {nome_corpo_maggiore}\nTempo di rivoluzione: {(((T/60)/60)/24):.3f} d = {((((T/60)/60)/24)/365):.3f} y\nEccentricità: {e}")
elif rotazione_switch is True:
    plt.plot(x_orbita,
            y_orbita,
            label=f"Orbita di {nome_corpo_minore} attorno a {nome_corpo_maggiore}\nTempo di rivoluzione: {(((T/60)/60)/24):.3f} d = {((((T/60)/60)/24)/365):.3f} y\nNumero di orbite attorno a {nome_corpo_maggiore} in un giorno: {(numero_orbite_giorno):.3f}\n(Periodo di rotazione di {nome_corpo_maggiore}: {(((rot/60)/60)/24):.3f} d = {((((rot/60)/60)/24)/365):.3f} y)\nEccentricità: {e}")
plt.scatter(distanza_corpomaggiore_centro,
            0,
            label=f"Posizione di {nome_corpo_maggiore}",
            zorder = 5,
            color="yellow")

if albedo_switch is False and satellite_switch is False:
    plt.scatter(raggio_afelio+distanza_corpomaggiore_centro,
                0,
                label=f"Afelio dell'orbita di {nome_corpo_minore}\nAfelio: {((raggio_afelio)/1000):.3f} km = {((raggio_afelio)/AU):.3f} AU = {((raggio_afelio)/ly):.3f} ly\nVelocità: {((velocità_afelio)*3.6):.3f} km/h",
                zorder = 5,
                color="green")
    plt.scatter(-(raggio_perielio-distanza_corpomaggiore_centro),
                0,
                label=f"Perielio dell'orbita di {nome_corpo_minore}\nPerielio: {((raggio_perielio)/1000):.3f} km = {((raggio_perielio)/AU):.3f} AU = {((raggio_perielio)/ly):.3f} ly\nVelocità: {((velocità_perielio)*3.6):.3f} km/h",
                zorder = 5,
                color="red")

elif albedo_switch is True and satellite_switch is False:
    plt.scatter(raggio_afelio+distanza_corpomaggiore_centro,
                0,
                label=f"Afelio dell'orbita di {nome_corpo_minore}\nAfelio: {((raggio_afelio)/1000):.3f} km = {((raggio_afelio)/AU):.3f} AU = {((raggio_afelio)/ly):.3f} ly\nVelocità: {((velocità_afelio)*3.6):.3f} km/h\nTemperatura: {((Temp_afelio)-273.15):.3f}°C",
                zorder = 5,
                color="green")
    plt.scatter(-(raggio_perielio-distanza_corpomaggiore_centro),
                0,
                label=f"Perielio dell'orbita di {nome_corpo_minore}\nPerielio: {((raggio_perielio)/1000):.3f} km = {((raggio_perielio)/AU):.3f} AU = {((raggio_perielio)/ly):.3f} ly\nVelocità: {((velocità_perielio)*3.6):.3f} km/h\nTemperatura: {((Temp_perielio)-273.15):.3f}°C",
                zorder = 5,
                color="red")

if stella_switch is True:
    plt.plot(abitabile_minimo_x+distanza_corpomaggiore_centro,
            abitabile_minimo_y,
            label = f"Limite minore dell'area abitabile di {nome_corpo_maggiore}\nRaggio = {(raggio_abitabile_minimo/AU):.3f} AU",
            zorder = 5,
            color = "black"
            )
    plt.plot(abitabile_massimo_x+distanza_corpomaggiore_centro,
            abitabile_massimo_y,
            label = f"Limite maggiore dell'area abitabile di {nome_corpo_maggiore}\nRaggio = {(raggio_abitabile_massimo/AU):.3f} AU",
            zorder = 5,
            color = "black"
            )
    plt.fill(x_R+distanza_corpomaggiore_centro,
            y_R,
            label = f"Area abitabile di {nome_corpo_maggiore}",
            color="lightgreen",
            alpha=0.4)

if satellite_switch is True:
    plt.scatter(raggio_afelio+distanza_corpomaggiore_centro,
                0,
                label=f"Apogeo dell'orbita di {nome_corpo_minore}\nApogeo: {((raggio_afelio-raggio_maggiore)/1000):.3f} km = {((raggio_afelio)/AU):.3f} AU\nVelocità: {((velocità_afelio)*3.6):.3f} km/h",
                zorder = 5,
                color="green")
    plt.scatter(-(raggio_perielio-distanza_corpomaggiore_centro),
                0,
                label=f"Perigeo dell'orbita di {nome_corpo_minore}\nPerigeo: {((raggio_perielio-raggio_maggiore)/1000):.3f} km = {((raggio_perielio)/AU):.3f} AU\nVelocità: {((velocità_perielio)*3.6):.3f} km/h",
                zorder = 5,
                color="red")
    plt.plot(raggio_x_corpomaggiore+distanza_corpomaggiore_centro,
            raggio_y_corpomaggiore,
            color="black",
            zorder=5)
    plt.fill(raggio_x_corpomaggiore+distanza_corpomaggiore_centro,
            raggio_y_corpomaggiore,
            color="red",
            label=f"Area occupata da {nome_corpo_maggiore}\nRaggio: {(raggio_maggiore):.3f} m",
            alpha=0.4)

plt.grid()
plt.axis("equal") #evita distorsione cerchi
plt.legend(loc="upper right", fontsize=7.5)
plt.xlabel("Raggio (m)")
plt.ylabel("Raggio (m)")
plt.title(f"Orbita di {nome_corpo_minore} attorno a {nome_corpo_maggiore}")
plt.show()