Dear sentaro and all, my most advanced project is the simulation of an open cluster with the calculator. If you let a bunch of stars do what they want, you see interesting effects. Which occur in reality too
- development of a dense core and a halo
- some of the stars escape
- you get double stars, which generate energy = velocity for the other stars
Double stars are a remarkable numerical problem, because I have a time step control. The simulation of double stars requires small time steps.
I changed my comments into English, but I let the changed version not run. May be you have to eliminate tiny problems lie a missing ' or so.
The program recommends 50 stars, but I worked with larger numbers too.
The first such calculations were done by Sebastian Hoerner starting around 1951 with very few bodies. Sverre Aarseth dedicated all of his professional life to that problems. His Fortran programs are very advanced and freely available.
Here is a wiki site to the problem.
https://en.wikipedia.../N-body_problem
~
If you run the program you get a graph with the traces of the stars. In the top line is given:
time (ahead the .)number of escaper (after the .)Wkin/Wpot, that is the relation kinetic to potential energy
'ProgramMode:RUN
'nbody, Open and galactic clusters. Leapfrog method.
' Stop for looking at results: line 236 R/S
' w0=|r| w1=|v| r0-2=Koor v0-2:velocity
' t = sqrt(rrr/Gm) t*.
't(Mio.a)=sqrt(6.4rrr/m) r/Ly m/mS
Goto M 'Main
Lbl i'init
ClrGraph
0->X
0->Y
1->G
2->l 'Limit for escapers
10->N
"N"?N
1->Q
1->M
0.136*Sqrt(N)->W'Relax.
30->S
"final time"?S
0.03->J'use J=0.03 for balance between accuracy and calculation time
"J"?J
1->g
0.1*Sqrt(1/N)->d
0.05/(N^.333)->a'Plummer radius reduces creation of double stars
{N,3}->Dim Mat r
{N,3}->Dim Mat v
{N,3}->Dim Mat w
{N,3}->Dim Mat a
{N,3}->Dim Mat b
{N,1}->Dim Mat m
'Startwerte
0->i
While i<N
0->R
0->V
For 0->j To 2
Ran#->Z
2Q*(Z-0.5)->r[i,j]
R+r[i,j]^<2>->R
Ran#->Z
2*(Z-0.5)->v[i,j]
V+v[i,j]^<2>->V
Next
SqrtR->R
SqrtV->V
If (R<Q) And (V<1)
M/N->m[i]
i+1->i
IfEnd
WhileEnd
Gosub b 'baryz R
'virialisieren
Gosub k
Gosub p
Sqrt(2*K/P)->c
For 0->i To N-1
For 0->k To 2
v[i,k]/c->v[i,k] 'virial.
0->v[i,k] 'cc
Next
Next
Gosub B'baryz V
Gosub a'first acceleration here
Return
Lbl g'Grafik
For 0->i To N-1
Blue Plot r[i,0]/Q,r[i,1]/Q
Next
Return
Lbl a'acceleration
For 0->p To N-1
For 0->k To 2
a[p,k]->b[p,k]
0->a[p,k]
Next
Next
0->B
For 0->p To N-1
For p+1->q To N-1
Sqrt((r[p,0]-r[q,0])^<2>+(r[p,1]-r[q,1])^<2>+(r[p,2]-r[q,2])^<2>)->R
R(RR+aa)->E
For 0->k To 2
G*(r[q,k]-r[p,k])/E->A
a[p,k]+m[q]*A->a[p,k]
a[q,k]-m[p]*A->a[q,k]
Next'k
Next'q
Sqrt((a[p,0]-b[p,0])^<2>+(a[p,1]-b[p,1])^<2>+(a[p,2]-b[p,2])^<2>)->j
Sqrt(a[p,0]^<2>+a[p,1]^<2>+a[p,2]^<2>)->A
If j/A>B
j/A->B
IfEnd
Next'p
d->D
If B>0
J/B*d->D
IfEnd
Return
Lbl p'Wpot
0->P
For 0->i To N-1
For i+1->j To N-1
0->R
For 0->k To 2
(r[i,k]-r[j,k])^<2>+R->R
Next
Sqrt(R)->R
P+m[i]*m[j]/R->P
Next
Next
Return
Lbl k'Wkin
0->K
For 0->i To N-1
v[i,0]^<2>+v[i,1]^<2>+v[i,2]^<2>->V
K+m[i]/2*V->K
Next
Return
Lbl b'baryzentrisch R
0->x:0->y:0->z:0->L
For 0->i To N-1
{i,x,r[i,0]}
m[i]*r[i,0]+x->x
m[i]*r[i,1]+y->y
m[i]*r[i,2]+z->z
m[i]+L->L
Next
For 0->i To N-1
r[i,0]-x/L->r[i,0]
r[i,1]-y/L->r[i,1]
r[i,2]-z/L->r[i,2]
Next
Return
Lbl B'baryzentrisch V
0->x:0->y:0->z
For 0->i To N-1
v[i,0]+x->x
v[i,1]+y->y
v[i,2]+z->z
Next
For 0->i To N-1
v[i,0]-x/N->v[i,0]
v[i,1]-y/N->v[i,1]
v[i,2]-z/N->v[i,2]
Next
Return
Lbl l'Leapfrog1
For 0->i To N-1
For 0->k To 2
v[i,k]+a[i,k]*d/2->w[i,k]
Next
Next
For 0->i To N-1
For 0->k To 2
r[i,k]+w[i,k]*d->r[i,k]
Next
Next
Gosub a
For 0->i To N-1
For 0->k To 2
w[i,k]+a[i,k]*d/2->v[i,k]
Next
Next
Return
Lbl L'Leapfrog
While t<=T
Gosub l
t+d->t
If g>0.1
D->d
IfEnd
WhileEnd
Return
Lbl G'Init graf
ClrGraph:CoordOff:AxesOff:GridOff:LabelOff:S-L-Normal:Rad:'#_SketchThin _
'_AxesOn_
' x0 x\1 dx y0 y\1 dy
ViewWindow (-)7.2,7.2,1,(-)4,4,1
Return
'Print.XY / Grafik
Lbl P
Text 1,9," "
Text 1,9,X
Text 1,179,Y
Return
'Stop
Lbl S
Plot :Text 1,1,">"Disps
While Getkey<>31:WhileEnd'exit=47
Text 1,1," "
Plot :
Return
Lbl M
Gosub G
Gosub i
1000->o
"number of outputs"?o
Gosub g
0->X
Gosub P
S/o->s'values for output
0->t
While t<S
t+s->T
If (T>S)
S->T
IfEnd
Gosub L
Gosub g
t->X
Gosub P
0->n ' escaper
For 0->i To N-1
r[i,0]^<2>+r[i,1]^<2>+r[i,2]^<2>->R
v[i,0]^<2>+v[i,1]^<2>+v[i,2]^<2>->V
SqrtR->R
SqrtV->V
R->w[i,0]
Vc->w[i,1]
0.6->w[i,2] ' (O)=place (G)=velocity
R>l=>n+1->n
Next
t->X
Gosub k
Gosub p
100*n+K/P->Y
Gosub P 'R//S for looking at the results in w0,1
n>N/10-0.5=>Stop
WhileEnd
Gosub S
Stop
Edited by piu58, 15 December 2023 - 09:47 AM.
