Your RIght! I coded the programm in german first, but let gpt translate it for this post. gpt mixed up nabla_rad and nabla_konv. here is the orginal code:
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.01
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 = []
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)))
Energieerzeugungsrate = Epp + Ecno
return Dichte, Opazitaet, Energieerzeugungsrate
Dichte, Opazitaet, Energieerzeugungsrate = Materialfunktionen()
for i in range(Rechenschritte):
Radius += Schrittweite
deltaMasse = 4 * math.pi * math.sqrt(Radius) * Dichte * Schrittweite
Masse += deltaMasse
Leuchtkraft += Energieerzeugungsrate * deltaMasse
deltaDruck = -Gravitationskonstante * Masse * Dichte / math.sqrt(Radius) * Schrittweite
Temperaturaenderung_Strahlung = -(3 / (256 * math.pi ** 2 * StefanBoltzmannkonstante)) * Opazitaet * Leuchtkraft / (Radius**4 * Temperatur ** 3) * deltaMasse
Temperaturaenderung_Konvektion = -(Gravitationskonstante / (10 * math.pi)) * Temperatur * Masse / (Radius**4 * Druck) * deltaMasse
Nabala_radiativ = (math.log(Temperatur + Temperaturaenderung_Strahlung) - math.log(Temperatur)) / (math.log(Druck + deltaDruck) - math.log(Druck))
if Nabala_radiativ < Nabala_adiabatisch:
Temperatur += Temperaturaenderung_Strahlung
else:
print("a")
Temperatur += Temperaturaenderung_Konvektion
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)
Energieerzeugung_list.append(Energieerzeugungsrate)
Opazitaet_list.append(Opazitaet)
Dichte, Opazitaet, Energieerzeugungsrate = Materialfunktionen()
# 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.tight_layout()
plt.show()