Startseite
| Programme:
Wellenoptik
| Algorithmen,
Lektionen, Analysen
|
Unsere Vorläufer
| Kontakt, Datenschutz
Startseite
| Programme:
Wellenoptik
| Algorithmen,
Lektionen, Analysen
|
Unsere Vorläufer
| Kontakt, Datenschutz
Das Programm simuliert die dynamische Entwicklung eines Sternhaufens und ermittelt zyklisch dessen Kenngrößen
# 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
'''
Variablenliste
a Hilfsfeld A Beschleunigungen
b Hilfsfeld B größter Ruck
c Virialisierungsfaktor C %
d deltaT D Vorschlag für deltaT
e anzDS E R3
f F
g Flag Genauigkeitsstrg. G Gravkonstante=1
h Anz. Relaxationsztn bis End
i Schleife
j Ruck J erlaubte Grenze für Ruck
k Schleife K wKin
l limit Entweicher L korrigiertes Limit für Entweicher
m Masse M Gesamtmasse
n % N Anzahl Körper
o Anzahl Ausgaben
p Schleife P wPot
q Schleife Q R0
r Ortsvektoren R Abstand vom Origo
s Teilzeit(Ausgabe) S Endzeit
t simulierte Zeit T Endzeit für die LF-Funktion
u
v Geschwindigkeit(Feld) V v²
w % W Relaxationszeit
x
y
z % Zufallszahl
'''