Advertisement
Gondus

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

May 15th, 2024
484
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 7.28 KB | Science | 0 0
  1. from matplotlib import pyplot as grafico
  2. import numpy as np
  3. import math
  4. from scipy import special
  5.  
  6. print("Benvenuto nel calcolatore di orbite!")
  7. tipo_corpo = input("Specificare se il corpo minore è un satellite: ")
  8. nome = input("Inserire il nome del corpo minore: ")
  9. NOMONE = input("Inserire il nome del corpo maggiore: ")
  10. stella = input("Specificare se il corpo maggiore è una stella: ")
  11. M = float(input("Inserire la massa del corpo maggiore in kg (includere decimali): "))
  12. if stella.lower() == 'stella':
  13.     Lum = float(input("Inserire luminosità solare (includere decimali): "))
  14.     albedo_domanda = input("L'albedo del corpo minore è noto? ")
  15.     if albedo_domanda.lower() == "sì" or albedo_domanda.lower() == "si":
  16.         albedo = float(input("Inserire l'albedo del corpo minore, preferibilmente Bond (includere decimali): "))
  17. R_s = float(input("Inserire il raggio del corpo maggiore in m (includere decimali): "))
  18. R_p = float(input("Inserire il raggio del corpo minore in m (includere decimali): "))
  19. m = float(input("Inserire la massa del corpo minore in kg (includere decimali): "))
  20. T = float(input("Inserire il tempo di rivoluzione del corpo minore in s (includere decimali, in secondi): "))
  21. e = float(input("Inserire l'eccentricità dell'orbita (includere decimali): "))
  22. G = 6.67430e-11
  23. AU = 1.495978707e11
  24. sigma = 5.670374419e-8 #costante di Stefan-Boltzmann
  25.  
  26. a = np.cbrt(((T**2)*G*M)/(4*(np.pi**2)))
  27. #E = (-G*M*m)/(2*a)
  28. #L = m*(2*np.pi/T)*(a**2)
  29.  
  30.  
  31. #calcoli
  32. if stella.lower() == 'stella':
  33.     Rmin = np.sqrt(Lum/1.1)*AU
  34.     Rmax = np.sqrt(Lum/0.53)*AU
  35.     Lim_Roche = R_s*2.44*((M/(4/3*np.pi*R_s**3))/(m/(4/3*np.pi*R_p**3)))**(1/3)
  36. b = a*np.sqrt(1-e**2)
  37. if tipo_corpo.lower() == 'satellite':
  38.     rP = a*(1-e) - R_s
  39.     rA = a*(1+e) - R_s
  40.     v_h1 = np.sqrt((G*M)/(R_s+rP))
  41.     v_h2 = np.sqrt((G*M)/(R_s+rA))
  42. if tipo_corpo.lower() != 'satellite':
  43.     rP = a*(1-e)
  44.     rA = a*(1+e)
  45.     Lim_Roche = R_s*2.44*((M/(4/3*np.pi*R_s**3))/(m/(4/3*np.pi*R_p**3)))**(1/3)
  46. if stella.lower() == 'stella':
  47.     if albedo_domanda.lower() == "sì" or albedo_domanda.lower() == "si":
  48.         L_stella = Lum*3.828e26
  49.         T_p = ((L_stella*(1-albedo))/(16*np.pi*sigma*(rP)**2))**(1/4)
  50.         T_a = ((L_stella*(1-albedo))/(16*np.pi*sigma*(rA)**2))**(1/4)
  51. dF = e*a
  52. rA_AU = rA/AU
  53. rP_AU = rP/AU
  54. v_p = np.sqrt((G*M*(1+e))/(a*(1-e)))
  55. v_a = np.sqrt((G*M*(1-e))/(a*(1+e)))
  56.  
  57. #Costruiamo l'orbita
  58. angolo = np.linspace(0, 2*np.pi)
  59. x = a*np.cos(angolo)
  60. y = b*np.sin(angolo)
  61. if stella.lower() == 'stella':
  62.     Rmin_x = Rmin*np.cos(angolo) #abitabile min
  63.     Rmin_y = Rmin*np.sin(angolo)
  64.     Rmax_x = Rmax*np.cos(angolo) #abitabile max
  65.     Rmax_y = Rmax*np.sin(angolo)
  66.     #concatenazione per fare percorso chiuso
  67.     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
  68.     y_R = np.concatenate([Rmin_y, Rmax_y[::-1]])
  69. if tipo_corpo.lower() != "satellite":
  70.     Roche_x = Lim_Roche*np.cos(angolo)
  71.     Roche_y = Lim_Roche*np.sin(angolo)
  72.     Area_occupata_x = R_s*np.cos(angolo)
  73.     Area_occupata_y = R_s*np.sin(angolo)
  74. if tipo_corpo.lower() == 'satellite':
  75.     afelio = rA
  76.     perielio = -rP
  77.     xA = rA*np.cos(angolo)
  78.     xP = rP*np.sin(angolo)
  79. else:
  80.     afelio = rA-dF
  81.     perielio = -rP-dF
  82.  
  83. grafico.grid()
  84. if tipo_corpo.lower() == 'satellite':
  85.     orbita = grafico.plot(xA, xP, label="Orbita di " + nome + " attorno a " + NOMONE + ".", color = "blue")
  86. else:
  87.     orbita = grafico.plot(x, y, label="Orbita di " + nome + " attorno a " + NOMONE + ".", color = "blue")
  88. 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
  89. if tipo_corpo.lower() == 'satellite':
  90.     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")
  91.     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")
  92. else:
  93.     if stella.lower() == 'stella':
  94.         if albedo_domanda.lower() == "sì" or albedo_domanda.lower() == "si":
  95.             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")
  96.             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")
  97.         else:
  98.             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")
  99.             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")
  100.     else:
  101.         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")
  102.         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")
  103. if stella.lower() == 'stella':
  104.     raggio_int = grafico.plot(Rmin_x-dF, Rmin_y, label="Raggio minimo della zona abitabile di " + NOMONE + " (" + str(Rmin/AU) + " AU)", color = "black")
  105.     raggio_ext = grafico.plot(Rmax_x-dF, Rmax_y, label="Raggio massimo della zona abitabile di " + NOMONE + " (" + str(Rmax/AU) + " AU)", color = "black")
  106.     area = grafico.fill(x_R-dF, y_R, label="Area della zona abitabile di " + NOMONE, color = "lightgreen", alpha=0.4) #alpha rende più opaco
  107.     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)
  108.     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)
  109. if tipo_corpo.lower() != "satellite" and stella.lower() != "stella":
  110.     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)
  111.     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)
  112. grafico.axhline(y=0, color="grey", linestyle="solid")
  113. grafico.axvline(x=0, color="grey", linestyle="solid")
  114. grafico.axis("equal") #EVITA CHE I CERCHI SI DISTORGANO, centra gli assi
  115. grafico.title("Orbita ellittica di un corpo attorno a un altro")
  116. grafico.legend(loc="upper right", fontsize=7.5) #Mostra i nomi delle funzioni presentate e i colori, in alto a destra
  117. grafico.xlabel("Raggio (m)")
  118. grafico.ylabel("Raggio (m)")
  119. grafico.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement