ich bin gerade dabei ein Sternenmodell zu programmieren, und bekomme unerwartete Probleme bei dessen Ergebnissen:
Code: Alles auswählen
import matplotlib.pyplot as plt
from astropy import constants as const
import math
# Eingabewerte
Sonnenmassen = 2
Zentraltemperatur = 15710000
Zentraldruck = 2.477E15
Wasserstoffteilchenzahl = 0.75
Stickstoffteilchenzahl = 0.10
Rechenschritte = 100
# Konstanten
Gaskonstante = const.R.cgs.value
Gravitationskonstante = const.G.cgs.value
StefanBoltzmannkonstante = const.sigma_sb.cgs.value
Lichtgeschwindigkeit = const.c.cgs.value
# Zielwerte
MaxMasse = Sonnenmassen * const.M_sun.cgs.value
if Sonnenmassen < 1.2:
MaxRadius = const.R_sun.cgs.value * Sonnenmassen ** 0.7
else:
MaxRadius = const.R_sun.cgs.value * Sonnenmassen ** 0.6
MaxLeuchtkraft = const.L_sun.cgs.value * Sonnenmassen ** 3.5
Schrittweite = MaxRadius / Rechenschritte
MittleresMolekulargewicht = 0.5 * Wasserstoffteilchenzahl + 4 / 3 * (
1 - Wasserstoffteilchenzahl - Stickstoffteilchenzahl) + 2 * Stickstoffteilchenzahl
Nabala_adiabatisch = 0.4
# Werte der Schleife
Radius = 0
Masse = 0
Leuchtkraft = 0
Temperatur = Zentraltemperatur
Druck = Zentraldruck
Nabala_radiativ = 0
# Listen für die Plots
radius_list = []
masse_list = []
leuchtkraft_list = []
temperatur_list = []
druck_list = []
dichte_list = []
Energieerzeugung_list = []
Opazitaet_list = []
nabla_list = []
def Materialfunktionen():
Dichte = (MittleresMolekulargewicht / Gaskonstante) * (Druck / Temperatur - (4 * StefanBoltzmannkonstante * Temperatur ** 3) / (3 * Lichtgeschwindigkeit))
Opazitaet = 2.5E25 * Dichte / Temperatur ** 3.5
Epp = 5.4879E9 * Wasserstoffteilchenzahl ** 2 * Dichte * Temperatur ** (-2 / 3) * math.exp(
-3380 * Temperatur ** (-1 / 3))
Ecno = (8E31 * Wasserstoffteilchenzahl * Stickstoffteilchenzahl * Dichte * Temperatur ** (-2 / 3) *
math.exp(1.88E8 * math.sqrt(1.5 + Wasserstoffteilchenzahl / 2) * math.sqrt(Dichte) * Temperatur ** (
-1.5) - 15228 * Temperatur ** (-1 / 3)))
Energieerzeugung = Epp + Ecno
print(Energieerzeugung)
return Dichte, Opazitaet, Energieerzeugung
for i in range(Rechenschritte):
Dichte, Opazitaet, Energieerzeugung = Materialfunktionen()
Radius += Schrittweite
deltaMasse = 4 * math.pi * Radius**2 * Dichte * Schrittweite
Masse += deltaMasse
deltaDruck = -(Gravitationskonstante * Masse * Dichte / Radius ** 2) * Schrittweite
Leuchtkraft += deltaMasse * Energieerzeugung
Temperaturaenderung_Strahlung = -(3 / (64 * math.pi * StefanBoltzmannkonstante)) * ((Opazitaet * Dichte)/Radius ** 2) * ((Leuchtkraft) / Temperatur ** 3) * Schrittweite
Temperaturaenderung_Konvektion = (-2/5) * (Temperatur/Druck) * deltaDruck
Nabala_radiativ = (math.log(Temperatur + Temperaturaenderung_Strahlung) - math.log(Temperatur)) / (math.log(Druck + deltaDruck) - math.log(Druck))
if abs(Nabala_radiativ) < abs(Nabala_adiabatisch):
Temperatur += Temperaturaenderung_Strahlung
else:
Temperatur += Temperaturaenderung_Konvektion * Schrittweite
Druck += deltaDruck
# Hinzufügen der Werte zur Liste
radius_list.append(Radius)
masse_list.append(Masse)
leuchtkraft_list.append(Leuchtkraft)
temperatur_list.append(Temperatur)
druck_list.append(Druck)
dichte_list.append(Dichte)
Opazitaet_list.append(Opazitaet)
Energieerzeugung_list.append(Energieerzeugung)
nabla_list.append(Nabala_radiativ)
# Plotting
plt.figure(figsize=(12, 16)) # Ändere die Größe des Plots, um Platz für die zusätzlichen Plots zu machen
plt.subplot(4, 2, 1)
plt.plot([r / MaxRadius for r in radius_list], [m / MaxMasse for m in masse_list], color='blue')
plt.xlabel('Radius (%)')
plt.ylabel('Masse (%)')
plt.subplot(4, 2, 2)
plt.plot([r / MaxRadius for r in radius_list], [l / MaxLeuchtkraft for l in leuchtkraft_list], color='orange')
plt.xlabel('Radius (%)')
plt.ylabel('Leuchtkraft (%)')
plt.subplot(4, 2, 3)
plt.plot([r / MaxRadius for r in radius_list], temperatur_list, color='green')
plt.xlabel('Radius (%)')
plt.ylabel('Temperatur')
plt.subplot(4, 2, 4)
plt.plot([r / MaxRadius for r in radius_list], druck_list, color='red')
plt.xlabel('Radius (%)')
plt.ylabel('Druck')
plt.subplot(4, 2, 5)
plt.plot([r / MaxRadius for r in radius_list], dichte_list, color='purple')
plt.xlabel('Radius (%)')
plt.ylabel('Dichte')
plt.subplot(4, 2, 6)
plt.plot([r / MaxRadius for r in radius_list], Energieerzeugung_list, color='brown')
plt.xlabel('Radius (%)')
plt.ylabel('Energieerzeugung')
plt.subplot(4, 2, 7)
plt.plot([r / MaxRadius for r in radius_list], Opazitaet_list, color='magenta')
plt.xlabel('Radius (%)')
plt.ylabel('Opazität')
plt.subplot(4, 2, 8)
plt.plot([r / MaxRadius for r in radius_list], nabla_list, color='magenta')
plt.xlabel('Radius (%)')
plt.ylabel('nabla')
plt.tight_layout()
plt.show()
https://imgur.com/AozqtnG
https://imgur.com/9pzSmBQ
https://imgur.com/XEI59yz
https://imgur.com/bNP6Xrl
[https://imgur.com/JdCc7LH
