I used the Runge-Kutta method for simulation. You can try different number of bodies, initial velocities and time interval between two positions.
Code:
Code: Select all
'===============================================================================
' Simulation of a symmetric n-body problem
' NBody_symmetric.bas
' Created on May, 02 2024
'===============================================================================
Dim Shared As Integer n = 2, k
Dim Shared As Double v0 = 0.76, xymax = 3.0, dt = 1e-5, _
colour(), x0(), y0(), vx0(), vy0(), x1(), y1(), vx1(), vy1()
Dim As String key
Function Edit(ByVal text As String) As String
'Simple line editor to edit a predefined text (only BACKSPACE And ENTER are supported)
Dim As Integer row, col
Dim As String s, key
row = CsrLin
col = Pos
s = text
Do
Locate row, col
Print Space(Len(s) + 1)
Locate row, col
Print s;
Do
Sleep 10
key = InKey
Loop Until Len(key) = 1
Select Case ASC(key)
Case 13
Exit Do
Case 8
s = Left(s, Len(s) - 1)
Case Else
If Asc(key) > 30 Then s = s & key
End Select
Loop
Print 'new line
Return s
End Function
Sub Console()
'Console window with information and data input
Screen 0
Width 80, 25
Shell "Title " & "Simulation of a symmetric n-body problem"
Print "Initial conditions: a number of bodies are located at equal distances on "
Print "a circle and move along a tangent to this circle with the same initial "
Print "velocity. In the simulation, the radius of the circle, all masses and the "
Print "gravity constant are set to 1. The user can select the number of bodies "
Print "(2 to 7), the initial velocity, the time interval between two calculated "
Print "positions of each body and the mathematical size of the graphic window. "
Print "The program starts with a set of parameters for 3 bodies. The initial "
Print "velocity should be increased when the number of bodies is increased."
Print
Print "Edit parameters here. You can use the BACKSPACE and the ENTER key."
Print "Number of bodies (2 to 7): ";
n = Val(Edit(Str(n + 1))) - 1
Print "Initial velocity: ";
v0 = Val(Edit(Str(v0)))
Print "Time interval: ";
dt = Val(Edit(Str(dt)))
Print "Size of graphic window +/-";
xymax = Val(Edit(Str(xymax)))
ReDim colour(n), x0(n), y0(n), vx0(n), vy0(n), x1(n), y1(n), vx1(n), vy1(n)
End Sub
Function Func(Byval i As Integer, Byval t As Double, u() As Double) As Double
'Function used by the Sub DiffSysn. For each body with index k (= 0 to n), the
'following system of differential equations is solved:
'dx / dt = vx
'dy / dt = vy
'dvx / dt = -Sum( m(j) * (x(k) - x(j)) / r(k, j)^3 ), j = 0 to n, j <> k
'dvy / dt = -Sum( m(j) * (y(k) - y(j)) / r(k, j)^3 ), j = 0 to n, j <> k
'r(j, k) = sqr((x - x0(j))^2 + (y - y0(j))^2
Dim As Double x, y, vx, vy, r, sum
Dim As Integer j
x = u(1)
y = u(2)
vx = u(3)
vy = u(4)
Select Case i
Case 1
Return vx
Case 2
Return vy
Case 3
sum = 0
For j = 0 To n
If j <> k Then
r = Sqr((x - x0(j))^2 + (y - y0(j))^2)
sum = sum - (x - x0(j)) / r^3
End If
Next
Return sum
Case 4
sum = 0
For j = 0 To n
If j <> k Then
r = Sqr((x - x0(j))^2 + (y - y0(j))^2)
sum = sum - (y - y0(j)) / r^3
End If
Next
Return sum
End Select
End Function
Sub Diffsysn(ByVal x0 As Double, y0() As Double, ByVal dx As Double, Byref x1 As Double, y1() As Double)
'Stepwise solution of a system of first-order differential equations
'dy(i)/dx = Func(i, x, y(1), ... y(n)), i = 1 to n
'(Runge-Kutta method)
'Parameters:
'- x0, y0(1) to y0(n) = initial values
'- dx = stepwith
'- x1 = x0 + dx
'- y1(1) to y1(n) = results at x1
'Required function:
'Function Func(i As Integer, x As Double, y() As Double)
'Select Case i
'Case 1: Func = ... (Functions of x, y(1) to y(n))
'Case 2: Func = ...
'....
'End Select
'End Function
Dim As Integer n, i
n = Ubound(y0)
Dim As Double k1(1 To n), k2(1 To n), k3(1 To n), k4(1 To n), z(1 To n)
x1 = x0 + dx
For i = 1 To n
k1(i) = dx * Func(i, x0, y0())
Next i
For i = 1 To n
z(i) = y0(i) + k1(i) / 2
Next i
For i = 1 To n
k2(i) = dx * Func(i, x0 + dx / 2, z())
Next i
For i = 1 To n
z(i) = y0(i) + k2(i) / 2
Next i
For i = 1 To n
k3(i) = dx * Func(i, x0 + dx / 2, z())
Next i
For i = 1 To n
z(i) = y0(i) + k3(i)
Next i
For i = 1 To n
k4(i) = dx * Func(i, x1, z())
Next i
For i = 1 To n
y1(i) = y0(i) + (k1(i) + 2 * k2(i) + 2 * k3(i) + k4(i)) / 6
Next i
End Sub
Sub Simulation()
Dim As Integer i, xp(n), yp(n)
Dim As Any Ptr img(n)
Dim As Double pi = 4 * Atn(1), phi, u0(4), u1(4), t0 = 0, t1
'Open graphic window :
Screenres 600,650
Windowtitle "Simulation"
Width 600\8, 650\16
Locate Hiword(Width), 10
Print "Stop: Press any key";
'Image buffers for bodies:
For i = 0 To n
colour(i) = 9 + i
img(i) = Imagecreate(20, 20, 0)
Circle img(i), (10, 10), 9, colour(i),,,,F
Next
'Define initial conditions:
For i = 0 To n
phi = pi / 2 - 2 * i * pi / (n + 1)
x0(i) = Cos(phi)
y0(i) = Sin(phi)
vx0(i) = v0 * Cos(phi - pi / 2)
vy0(i) = v0 * Sin(phi - pi / 2)
Next
View (0, 0) - (599, 599),, 3
Window (- xymax, - xymax) - ( xymax, xymax)
Do
'Draw bodies:
For k = 0 To n
xp(k) = Pmap(x0(k), 0) - 10
yp(k) = Pmap(y0(k), 1) - 5
Put (Pmap(xp(k), 2), Pmap(yp(k), 3)), img(k), XOR
Next
'Calculate next positions:
For k = 0 To n
u0(1) = x0(k)
u0(2) = y0(k)
u0(3) = vx0(k)
u0(4) = vy0(k)
DiffSysn(t0, u0(), dt, t1, u1())
x1(k) = u1(1)
y1(k) = u1(2)
vx1(k) = u1(3)
vy1(k) = u1(4)
Next
'Delete bodies and draw actual position as a point:
For k = 0 To n
Put (Pmap(xp(k), 2), Pmap(yp(k), 3)), img(k), Xor
PSet (x0(k), y0(k)), colour(k)
Next
'The new positions are used as the next starting points:
t0 = t1
For k = 0 To n
x0(k) = x1(k)
y0(k) = y1(k)
vx0(k) = vx1(k)
vy0(k) = vy1(k)
Next
Loop Until Inkey <> ""
For k = 0 To n
ImageDestroy img(k)
Next
End Sub
'Main:
Do
Console()
Simulation()
Locate Hiword(Width), 10
Print Space(50);
Locate HiWord(Width), 10
Input; "Again y/n: ", key
Loop Until LCase(key) = "n"
End
