NBody UI (Thanks Jack for the code)

Everything else that doesn't fall into one of the other PB categories.
User avatar
pdwyer
Addict
Addict
Posts: 2813
Joined: Tue May 08, 2007 1:27 pm
Location: Chiba, Japan

NBody UI (Thanks Jack for the code)

Post 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

;====================================================================
Paul Dwyer

“In nature, it’s not the strongest nor the most intelligent who survives. It’s the most adaptable to change” - Charles Darwin
“If you can't explain it to a six-year old you really don't understand it yourself.” - Albert Einstein
User avatar
mk-soft
Always Here
Always Here
Posts: 5313
Joined: Fri May 12, 2006 6:51 pm
Location: Germany

Re: NBody UI (Thanks Jack for the code)

Post 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

;====================================================================
My Projects ThreadToGUI / OOP-BaseClass / EventDesigner V3
PB v3.30 / v5.75 - OS Mac Mini OSX 10.xx - VM Window Pro / Linux Ubuntu
Downloads on my Webspace / OneDrive
jack
Addict
Addict
Posts: 1336
Joined: Fri Apr 25, 2003 11:10 pm

Re: NBody UI (Thanks Jack for the code)

Post by jack »

nice :)
thanks for sharing
User avatar
pdwyer
Addict
Addict
Posts: 2813
Joined: Tue May 08, 2007 1:27 pm
Location: Chiba, Japan

Re: NBody UI (Thanks Jack for the code)

Post 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
Paul Dwyer

“In nature, it’s not the strongest nor the most intelligent who survives. It’s the most adaptable to change” - Charles Darwin
“If you can't explain it to a six-year old you really don't understand it yourself.” - Albert Einstein
Post Reply