Liebe Sternfreunde, das Heft 92 sollte demnächst in eurem Kasten liegen. Ich habe diesmal ein stark astrophysikalisch durchdrungenes Thema als Schwerpunkt. Ich fand es sehr erkenntnisreich, die Simulation von Sternhaufen slebst nachzurechnen. Hier an diese Mail klebe ich das Python-Programm dran. Ihr braucht das Programm
plot.py dazu, es muss im selben Verzeichnis liegen
Die Rechenzeiten sind lang, manches lief bei mir mehrere Tage.
Dagegen gibt es zwei "Medikamente":
1) Ich habe denselben Algorithmus in C++ vorliegen. Den müsst ihr aber auf eurem Computer übersetzen, Gnu-C++ geht. Wer das haben möchte: Mich bitte anschreiben.
2) Das sehr weit entwickelte
Programm nbody6 von Sverre Aarseth. Ich kann den leicht modifizierten Fortran-Code weggeben (einen Fehler gefunden und die Ausgabe angepasst) oder eine ausführbare Datei für Windows. Gnu-Fortran genügt zum übersetzen, zum Benutzen braucht ihr sicher Unterstützung von mir.
Sverre hat sein ganzes Berufsleben daran gearbeitet und war ein Vorreiter freier Software. Er ist übrigens
vor wenigen Tagen verstorben. Ich hatte zuletzt 2022 zuletzt Kontakt zu ihm. Seinerzeit fiel es ihm schon schwer, am Bildschirm etwas zu lesen.
~
Das PythonProgramm findet sich auch auf unserer
Internetseite.
Code:
# nbodyCluster (ohne Doppelsternbehandlung): Analyse der Entwicklung von Sternhaufen
from math import *
from random import *
from plot import *
from sys import *
G = 1 # Gravitationskonstante ist 1 gesetzt (einheitenfreies rechnen)
N = 100 # Anzahl der Körper
Q = 1 # Anfangsradius R0
M = 1 # Gesamtmasse
o = 1000 # Anzahl Ausgabewerte gesamt.
J = 0.03 # maximaler Ruck, J<0.03
g = 1 # Genauigkeitssteuerung an/aus
S = 30 # t* bis Simulationsende
l = 3.0 # Limit für Entweicher. 3 ist in Ordnung
L = l
xNext = 3 # der drittnächste Stern für die Dichteberechnung
E=100 # Abbruch bei xx % Entweichern
rP = 0.05 / (N ** 0.333) # Plummer-Radius, spart Rechenzeit und erschwert das Entstehen von Doppelsternen
d = 0.1 * sqrt(1/N) # deltaT vor der ersten Steuerung. (#8-22)
grafik=1 # nur Text oder auch Grafik? 0/1
r = [[0.0 for i in range(3)] for j in range(N)]
v = [[0.0 for i in range(3)] for j in range(N)]
w = [[0.0 for i in range(3)] for j in range(N)] # die Hilfsgeschwindigkeit v05
a = [[0 for i in range(3)] for j in range(N)]
b = [[0.0 for i in range(3)] for j in range(N)]
m = [0.0 for i in range(N)]
R = [0.0 for i in range(N)]
def sq(x):
return x * x
def baryzentrischR(N):
x = 0
y = 0
z = 0
LL = 0
for i in range(N): # Schwerpunkt berechnen ...
x+=m[i] * r[i][0]
y+=m[i] * r[i][1]
z+=m[i] * r[i][2]
LL+=m[i]
for i in range(N): # ... korrigieren
r[i][0]-=x / LL
r[i][1]-=y / LL
r[i][2]-=z / LL
def baryzentrischV(N):
x = 0
y = 0
z = 0
for i in range(N): # mittl. Geschwindigkeit ber...
x+=v[i][0]
y+=v[i][1]
z+=v[i][2]
for i in range(N): # ... korrigieren
v[i][0]-=x / N
v[i][1]-=y / N
v[i][2]-=z / N
def wPot(N):
P = 0
for i in range(N):
for j in range(i + 1,N):
R = 0
for k in range(3):
R = R + sq(r[i][k] - r[j][k])
R = sqrt(R)
P = P + m[i] * m[j] / R
P = G * P
return (P)
def wKin(N):
K = 0
for i in range(N):
V = sq(v[i][0]) + sq(v[i][1]) + sq(v[i][2])
K+=m[i] / 2 * V
return K
def beschleunigung():
global d
global D
global c
for p in range(N): # alle Beschleunigungen löschen / kopieren
for k in range(3):
b[p][k] = a[p][k]
a[p][k] = 0
B = 0
for p in range(N): # alle Körper
for q in range(p + 1, N): # gegen alle anderen Körper
R = sqrt(sq(r[p][0] - r[q][0]) + sq(r[p][1] - r[q][1]) + sq(r[p][2] - r[q][2]))
# E = R * R * R # E ist die umgekehrte 3, also R3
E = R * (R * R + rP * rP) # mit Plummer-Radius
for k in range(3):
A = G * (r[q][k] - r[p][k]) / E
a[p][k]+=m[q] * A
a[q][k]-=m[p] * A
# beide Beschleunigungen existieren jetzt, größten Ruck ber.
j = sqrt(sq(a[p][0] - b[p][0]) + sq(a[p][1] - b[p][1]) + sq(a[p][2] - b[p][2]))
A = sqrt(sq(a[p][0]) + sq(a[p][1]) + sq(a[p][2]))
if (j / A > B): # größten Ruck in B aufheben
B = j / A
# neuen Zeitschritt ausrechnen
D = d
if B > 0:
D = J / B * d # Vorschlag für den neuen Zeitschritt
# nach Hermite-Verfahren
def startwerte():
# seed(31) # seed / random nur für Tests
A = 6.0 / 5.0 # für Radius
for i in range(N):
R = A * ( random()**0.33333333)
cos_theta = 2.0 * random() - 1.0
sin_theta = sqrt(1 - (cos_theta*cos_theta))
phi = 2 * pi * random()
r[i][0] = R * sin_theta * cos(phi)
r[i][1] = R * sin_theta * sin(phi)
r[i][2] = R * cos_theta
A = 5.0 / 6.0 # für Geschwindigkeit
for i in range(N):
V = A * (random()**0.333333333)
cos_theta = 2.0 * random() - 1
sin_theta = sqrt(1 - (cos_theta*cos_theta))
phi = 2 * pi * random()
v[i][0] = V * sin_theta * cos(phi)
v[i][1] = V * sin_theta * sin(phi)
v[i][2] = V * cos_theta
for i in range(N):
m[i] = 1.0 / N
baryzentrischR(N)
# virialisieren
K = wKin(N)
P = wPot(N)
c = sqrt(2 * K / P)
for i in range(N):
for k in range(3):
v[i][k]/=c # /=c: virialisiert
v[i][k] = 0 # =0:kalter Kollaps. Auskommentieren für virialisierten Start
baryzentrischV(N)
beschleunigung() # die erste Beschleunigung gleich hier
return c
def leapfrog1():
# eine Hilfsgeschwindigkeit in Intervallmitte wird benötigt
for i in range(N):
for k in range(3):
w[i][k] = v[i][k] + a[i][k] * d / 2 # v05 = v + a * dt / 2
# der neue Ort nimmt die mittlere Geschwindigkeit
# Damit ist es ein quadratisches Verfahren
for i in range(N):
for k in range(3):
r[i][k] = r[i][k] + w[i][k] * d # s = s + v05 * dt
# die Geschwindigkeit muss korrigiert werden
beschleunigung()
for i in range(N):
for k in range(3):
v[i][k] = w[i][k] + a[i][k] * d / 2 # v=v05+a*dt/2
def leapfrog():
global t
global d
while t <= T:
leapfrog1()
t = t + d
if (g > 0.1): # neuer Zeitschritt nur, falls das Genauigkeitsflag gesetzt ist
d = D
def abstand(i, j):
return sqrt(sq(r[i][0] - r[j][0]) + sq(r[i][1] - r[j][1]) + sq(r[i][2] - r[j][2]))
def volumen(r):
return 1.333 * pi * r * r * r
def rC(N):
Next=[0.0 for i in range(xNext)]
rho = [[0.0 for i in range(2)] for j in range(N)]
for i in range(N): # Dichte um jeden Stern
for j in range (xNext):
Next[j] = 7e77
for j in range (N): # gegen jeden anderen Stern
dist = abstand(i, j)
for k in range (xNext):
if (i != j) and (dist < Next[k]):
for l in range(xNext - 1, k, -1):
Next[l] = Next[l - 1] # ein Element Platz schaffen
Next[k] = dist
break
rho[i][0] = R[i]
rho[i][1] = (xNext - 1) / volumen(Next[xNext-1])
# zentrale Dichte: die inneren 2%. das ist zweckmäßig bis N~500. Darüber: verkleinern.
anzC = int(0.02 * N + 0.5)
if (anzC < 2):
anzC = 2
rhoC = 0.0
for i in range(anzC):
rhoC += rho[i][1]
rhoC /= anzC
# Radius der viertel zentralen Dichte zurückmelden. Wird durch Doppelsterne verfälscht, hier nicht behoben.
radC=0
anzC=0
for i in range(N):
Rho = rho[i][1]
if (Rho < 0.25 * rhoC):
radC = rho[i][0]
anzC = i + 1
break
return radC, anzC, rhoC
# Doppelsterne: näher als 0.03, hängt nicht von N ab!
def DS(N):
dslimit = 0.03
ds = 0
for i in range(N): # Dichte jeden Stern prüfen
for j in range(i+1,N): # gegen jeden anderen Stern
dist = abstand(i, j)
if (dist < dslimit):
ds=ds+1
return ds/2 # jeder Stern wurde 2x berücksichtigt
def rH(n, t):
return R[int((N - n) / 2)] # nur valide für identische Massen aller Körper. Die Entweicher abziehen.
def ausgabe():
# Grafik
global L
n=0 # Anzahl Entweicher
for i in range(N):
if (grafik>0): Plot(r[i][0] / Q, r[i][1] / Q, 4, "black")
R[i] = sqrt(sq(r[i][0]) + sq(r[i][1]) + sq(r[i][2]))
if (R[i]>L): # L: Limit für Entweicher
n=n+1
R.sort()
radC, anzC, rhoC=rC(N)
radH=rH(n, t)
L = l * radH # Grenze für Entweicher aktualisieren
K = wKin(N)
P = wPot(N)
ds=DS(N)
# 1 2 3 4 5 6 7 8 9
# t* N rH wKin wPot rC nC dbl pC
print(t, N-n, radH, K, P, radC, anzC, ds, rhoC)
abbruch=0
if (n>(E/100.0)*N-0.5):
abbruch=1
return abbruch # Abbruch durch Entweicher
# Hauptprogramm
c = startwerte()
# Koordinatensystem einrichten
# x0 x1 xScale y0 y1 yScale
if (grafik>0): initKoor(-3 ,3 ,0.5, -2.5,2.5,0.5)
s = S / o # Zeitschritt für Mehrfach-Leapfrog
t = 0
clrgraf=0
while t < S:
T = t + s # Tend für leapfrog
if (T > S): T = S
clrgraf=clrgraf+s
if (clrgraf>1):
if (grafik>0): initKoor(-3 ,3 ,0.5, -2.5,2.5,0.5)
clrgraf=0
abbruch=ausgabe()
if (abbruch>0):
t=S+1
leapfrog()
if (grafik>0): update()
if (grafik>0): name = input("Enter!") # wg. Aufruf von Konsole und grep/plot
Registrieren
Anmelden
FAQ
Suche
