PureBasic Forums - English

http://www.purebasic.com
https://www.purebasic.fr/english/

NBody UI (Thanks Jack for the code)

https://www.purebasic.fr/english/viewtopic.php?t=79164

Page 1 of 1

NBody UI (Thanks Jack for the code)

Posted: Sun May 15, 2022 12:06 pm
by pdwyer
Jack's code from here:
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

;====================================================================

Re: NBody UI (Thanks Jack for the code)

Posted: Sun May 15, 2022 12:19 pm
by mk-soft
Nice,

but ist better to use window timer. My machine is to fast and save cpu power

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() 
  
  Protected EventID.i
  Protected Dayloop.i
  Protected EventGadget.i
  Protected do
  
  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())
              AddWindowTimer(#Frm_Main, 1, 10)
          EndSelect
          
        Case #PB_Event_Timer
          advance(#NBODIES, bodies(), 0.1)
          UpdateDisplay(dayloop) 
          
      EndSelect
    Until EventID = #PB_Event_CloseWindow
  EndIf
  
EndProcedure

;====================================================================

Re: NBody UI (Thanks Jack for the code)

Posted: Sun May 15, 2022 12:46 pm
by jack
nice :)
thanks for sharing

Re: NBody UI (Thanks Jack for the code)

Posted: Sun May 15, 2022 1:29 pm
by pdwyer
If you want to see older trails, if you set:

Code: Select all

VectorSourceColor(RGBA(0,0, 0, 40))
to

Code: Select all

VectorSourceColor(RGBA(255,255, 255, 40))
you can watch the change. (eg after tweaking a mass or x position so that the orbits cross and then wait for a bit till they planets nearly collide all sorts of fun things can happen

But with the black background it's hard to see how slight orbits are affected.

@Mk-soft. thanks for the timer idea, I'll re-use that I think
All times are UTC+01:00
Page 1 of 1

Powered by phpBB® Forum Software © phpBB Limited