https://www.purebasic.fr/english/viewtopic.php?t=77316
I was going to post there but it's the bug section
Jack,
Thank you for posting this code, I've been looking for nbody code and was thinking about whether to port some of the C code on the web as I was having trouble understanding how to implement from the math.
As a little test I put your code in a UI to see what happens when playing with the starting positions and mass etc.
I stripped out the 3rd dimention (z / vz) parts since I didn't need them and thought it might lighten the processing a bit.
The planet/sun sizes are by necessity arbitrary so there is something to see, I was thinking of making them proportional to mass but the planets disappear.
works well! planet movements look good when playing with the params. I'll try to add some features to the UI so I'm not always tweaking code and re-running
Code: Select all
;https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/nbody-gcc-5.html
; The Computer Language Benchmarks Game
; http://benchmarksgame.alioth.debian.org/
; contributed by Christoph Bauer
;
; https://benchmarksgame-team.pages.debian.net/benchmarksgame/license.html
;
; translated To PureBasic by Jack
; modified a bit To avoid using pointer To Array
; ; UI & 2D tweak by Paul
EnableExplicit
;=============================================================================
;-Constants
;=============================================================================
Enumeration
#Frm_Main
#Canvas
#cmd_Run
#MnuExit
EndEnumeration
#NBODIES = 5
#PI = 3.141592653589793
#solar_mass = 39.478417604357434 ;(4 * #pi) * #pi
#days_per_year = 365.24
Structure planet
x.d
y.d
;z.d
vx.d
vy.d
;vz.d
mass.d
size.l
col.l
EndStructure
Declare.q main(n.l)
Declare ShowFormMain()
Declare UpdateDisplay(frame.i)
;=============================================================================
;-Globals
;=============================================================================
Global gPlanetSize = 8
Global gWindowW.i = 1100
Global gWindowH.i = 950
Global gCanvasW.i = gWindowW - 110
Global gCanvasH.i = gWindowH - 20
Global gPopulation.i = 30
Global Dim bodies.planet(4)
main(1)
;====================================================================
Procedure advance (nbodies.l, Array bodies.planet(1), dt.d)
Define.l i, j
Define.d dx, dy, distance, mag ;dz,
For i = 0 To nbodies - 1
For j = i + 1 To nbodies - 1
dx = bodies(i)\x - bodies(j)\x
dy = bodies(i)\y - bodies(j)\y
;dz = bodies(i)\z - bodies(j)\z
distance = Sqr(((dx * dx) + (dy * dy)) ) ;+ (dz * dz)
mag = dt / ((distance * distance) * distance)
bodies(i)\vx = bodies(i)\vx - (dx * bodies(j)\mass) * mag
bodies(i)\vy = bodies(i)\vy - (dy * bodies(j)\mass) * mag
;bodies(i)\vz = bodies(i)\vz - (dz * bodies(j)\mass) * mag
bodies(j)\vx = bodies(j)\vx + (dx * bodies(i)\mass) * mag
bodies(j)\vy = bodies(j)\vy + (dy * bodies(i)\mass) * mag
;bodies(j)\vz = bodies(j)\vz + (dz * bodies(i)\mass) * mag
Next
Next
For i = 0 To nbodies - 1
bodies(i)\x = bodies(i)\x + dt * bodies(i)\vx
bodies(i)\y = bodies(i)\y + dt * bodies(i)\vy
; bodies(i)\z = bodies(i)\z + dt * bodies(i)\vz
Next
EndProcedure
;====================================================================
Procedure.d energy (nbodies.l, Array bodies.planet(1))
Define.l i, j
Define.d e, dx, dy, distance, mag ;dz,
e = 0.0
For i = 0 To nbodies - 1
e = e + (0.5 * bodies(i)\mass) * (((bodies(i)\vx * bodies(i)\vx) + (bodies(i)\vy * bodies(i)\vy)) ) ;+ (bodies(i)\vz * bodies(i)\vz)
For j = i + 1 To nbodies - 1
dx = bodies(i)\x - bodies(j)\x
dy = bodies(i)\y - bodies(j)\y
; dz = bodies(i)\z - bodies(j)\z
distance = Sqr(((dx * dx) + (dy * dy)) ) ;+ (dz * dz)
e = e - (bodies(i)\mass * bodies(j)\mass) / distance
Next
Next
ProcedureReturn e
EndProcedure
;====================================================================
Procedure offset_momentum (nbody.l, Array bodies.planet(1))
Define.d px, py;, pz
Define.l I
For I = 0 To nbody - 1
px = px + bodies(I)\vx * bodies(I)\mass
py = py + bodies(I)\vy * bodies(I)\mass
; pz = pz + bodies(I)\vz * bodies(I)\mass
Next
bodies(0)\vx = (-px) / ((4 * #PI) * #PI)
bodies(0)\vy = (-py) / ((4 * #PI) * #PI)
; bodies(0)\vz = (-pz) / ((4 * #PI) * #PI)
EndProcedure
;====================================================================
Procedure.q main (n.l)
bodies(0)\x = 0: bodies(0)\y = 0: bodies(0)\vx = 0: bodies(0)\vy = 0: bodies(0)\mass = (4 * #PI) * #PI ; bodies(0)\z = 0: bodies(0)\vz = 0:
bodies(1)\x = 4.84143144246472090e+00: bodies(1)\y = -1.16032004402742839e+00: bodies(1)\vx = 1.66007664274403694e-03 * #days_per_year: bodies(1)\vy = 7.69901118419740425e-03 * #days_per_year: bodies(1)\mass = 9.54791938424326609e-04 * ((4 * #PI) * #PI) ;bodies(1)\vz = (-6.90460016972063023e-05) * #days_per_year: bodies(1)\z = -1.03622044471123109e-01:
bodies(2)\x = 8.34336671824457987e+00: bodies(2)\y = 4.12479856412430479e+00: bodies(2)\vx = (-2.76742510726862411e-03) * #days_per_year: bodies(2)\vy = 4.99852801234917238e-03 * #days_per_year: bodies(2)\mass = 2.85885980666130812e-04 * ((4 * #PI) * #PI) ;bodies(2)\vz = 2.30417297573763929e-05 * #days_per_year: bodies(2)\z = -4.03523417114321381e-01:
bodies(3)\x = 1.28943695621391310e+01: bodies(3)\y = -1.51111514016986312e+01: bodies(3)\vx = 2.96460137564761618e-03 * #days_per_year: bodies(3)\vy = 2.37847173959480950e-03 * #days_per_year: bodies(3)\mass = 4.36624404335156298e-05 * ((4 * #PI) * #PI) ;bodies(3)\z = -2.23307578892655734e-01: bodies(3)\vz = (-2.96589568540237556e-05) * #days_per_year:
bodies(4)\x = 1.53796971148509165e+01: bodies(4)\y = -2.59193146099879641e+01: bodies(4)\vx = 2.68067772490389322e-03 * #days_per_year: bodies(4)\vy = 1.62824170038242295e-03 * #days_per_year: bodies(4)\mass = 5.15138902046611451e-05 * ((4 * #PI) * #PI) ;bodies(4)\z = 1.79258772950371181e-01: bodies(4)\vz = (-9.51592254519715870e-05) * #days_per_year:
bodies(0)\size = 20
bodies(1)\size = 5
bodies(2)\size = 2
bodies(3)\size = 3
bodies(4)\size = 4
bodies(0)\col = RGBA(250,225,165,255)
bodies(1)\col = RGBA(200,155,155,200)
bodies(2)\col = RGBA(100,155,255,200)
bodies(3)\col = RGBA(200,105,105,200)
bodies(4)\col = RGBA(100,255,055,200)
ShowFormMain()
EndProcedure
;=============================================================================
Procedure UpdateDisplay(frame.i)
StartVectorDrawing(CanvasVectorOutput(#Canvas))
;Blank out canvas
AddPathBox(0, 0, gCanvasW, gCanvasH)
VectorSourceColor(RGBA(0,0, 0, 40))
FillPath()
Define PlanetLoop.i
For PlanetLoop = 0 To #NBODIES -1
AddPathCircle(470 + (bodies(PlanetLoop)\x * 15), 470 + (bodies(PlanetLoop)\y * 15), bodies(PlanetLoop)\size) ;bodies(PlanetLoop)\mass
VectorSourceColor(bodies(PlanetLoop)\col)
FillPath()
Next
StopVectorDrawing()
EndProcedure
;====================================================================
Procedure ShowFormMain()
Define EventID.i
Define Dayloop.i
Define EventGadget.i
If OpenWindow(#Frm_Main, 200,100, gWindowW, gWindowH, "N-Body", #PB_Window_SystemMenu|#PB_Window_SizeGadget|#PB_Window_MinimizeGadget|#PB_Window_TitleBar|#PB_Window_MaximizeGadget)
CanvasGadget(#Canvas, 100, 10, gCanvasW, gCanvasH, #PB_Canvas_Border)
ButtonGadget(#cmd_Run, 10,10, 50, 18, "Run")
Repeat
EventID = WaitWindowEvent()
Select eventID
Case #PB_Event_Gadget
EventGadget = EventGadget()
Select EventGadget
Case #Canvas
Case #cmd_Run
offset_momentum(#NBODIES, bodies())
For Dayloop = 1 To 1500000
advance(#NBODIES, bodies(), 0.05)
UpdateDisplay(dayloop)
;nasty exit jump
If Mod(Dayloop,20) = 0
If WindowEvent() = #PB_Event_CloseWindow
Break 2
EndIf
EndIf
Next
Default
EndSelect
EndSelect
Until EventID = #PB_Event_CloseWindow
EndIf
EndProcedure
;====================================================================