Das deutsche Python-Forum

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()