Example: Section 8.4 (page 281), The Real 3 Body Problem

[Michael Clease]

[code=thinbasic]
' Example: Section 8.4 (page 281), The Real 3 Body Problem

' From Stan Blank's Book:
' "Python Programming in OpenGL
' "A Graphical Approach to Programming

' Converted by Michael Clease (Stolen timing ideas from Petr..hehe)
' Last modified: February 27, 2010

Uses "TBGL"

' Handle for our window
Global hWnd As DWord
Global Width As Long
Global Height As Long
Global vx1, vy1, vz1, x1, y1, z1, ax1, ay1, az1 As Double
Global vx2, vy2, vz2, x2, y2, z2, ax2, ay2, az2 As Double
Global vx3, vy3, vz3, x3, y3, z3, ax3, ay3, az3, dt, G As Double
Global m1,m2,m3, rad1 ,rad2, rad3 As Double
Global r12,r13,r23,r312,r23,r313,r323 As Double

Dim Title As String Value "A little MORE Gravity! 2 Bodies Simulation - 'Esc' to quit"
width = 500
height = 500
'
' Create and show window
hWnd = TBGL_CreateWindowEx(Title, Width, Height, 32, %TBGL_WS_WINDOWED )
TBGL_ShowWindow

Init()'Anim,x,y,dx,dy,axrng)
' Define "Callback" to be fired + that it should be fired each 10ms
TBGL_BindPeriodicFunction(hWnd, "PlotFunc", 10)

' -- Once the command below is executed, further script execution
' -- is halted and only periodic calling of the bound function is performed
TBGL_ProcessPeriodicFunction(hWnd)

TBGL_DestroyWindow

Sub Init()
x1 = 1.0
y1 = 1.0
z1 = 0.0
x2 = -1.0
y2 = -1.0
z2 = 0.0
x3 = 0.50
y3 = -1.0
z3 = 0.0
'# Initial values For velocity components In x,y,z space
vx1 = 0.0
vy1 = 0.0
vz1 = 0.0
vx2 = 0.0
vy2 = 0.0
vz2 = 0.0
vx3 = 0.0
vy3 = 0.0
vz3 = 0.0
'# Initial acceleration components
ax1 = 0.0
ay1 = 0.0
az1 = 0.0
ax2 = 0.0
ay2 = 0.0
az2 = 0.0
ax3 = 0.0
ay3 = 0.0
az3 = 0.0
'# Initial star masses
m1 = 0.7
m2 = 0.4
m3 = 0.5
'# Gravitational Constant
G = 1.0
'# radius of stars used In the plotFunc Function
rad1 = 0.2*m1
rad2 = 0.2*m2
rad3 = 0.2*m3
'# Calculate r**3 denominators For 3 Body Gravitation
'# More complex because the motion of EACH star depends
'# On where the other two stars are located!
r12 = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)
r312 = r12*Sqr(r12) + 0.01
r13 = (x1-x3)*(x1-x3) + (y1-y3)*(y1-y3) + (z1-z3)*(z1-z3)
r313 = r13*Sqr(r13) + 0.01
r23 = (x2-x3)*(x2-x3) + (y2-y3)*(y2-y3) + (z2-z3)*(z2-z3)
r323 = r23*Sqr(r23) + 0.01
'# Calculate the initial accelerations
'# MUCH more complex than 2 Body dynamics
'# Because each star must use the combined forces
'# due To gravity of the other 2 stars.
'# THIS is why there are TWO ax1, etc statements.
ax1 += -G*(x1-x2)*m2/r312
ax1 += -G*(x1-x3)*m3/r313
ay1 += -G*(y1-y2)*m2/r312
ay1 += -G*(y1-y3)*m3/r313
az1 += -G*(z1-z2)*m2/r312
az1 += -G*(z1-z3)*m3/r313
ax2 += -G*(x2-x1)*m1/r312
ax2 += -G*(x2-x3)*m3/r323
ay2 += -G*(y2-y1)*m1/r312
ay2 += -G*(y2-y3)*m3/r323
az2 += -G*(z2-z1)*m1/r312
az2 += -G*(z2-z3)*m3/r323
ax3 += -G*(x3-x2)*m2/r323
ax3 += -G*(x3-x1)*m1/r313
ay3 += -G*(y3-y2)*m2/r323
ay3 += -G*(y3-y1)*m1/r313
az3 += -G*(z3-z2)*m2/r323
az3 += -G*(z3-z1)*m1/r313
'# THIS Value keeps a smooth orbit On my workstation
'# Smaller values slow down the orbit, higher values speed things
'# up
dt = 0.01
TBGL_BackColor(0,0,0)
' Resets status of all keys
TBGL_ResetKeyState()
End Sub

Sub PlotFunc()
' -- Which window is calling?
Local hWnd As DWord = TBGL_CallingWindow
' TBGL_Camera( x, y, z, x2, y2, z2 )
ORBITS()
TBGL_ClearFrame '(%TBGL_CLEAR_Color)

TBGL_Camera(0.0, 0.0, 5.0, 0.0, 0.0, 0.0)

TBGL_PushMatrix
TBGL_Translate(x1,y1,z1)
TBGL_Color(245, 150, 30)
TBGL_Sphere(Rad1)
TBGL_PopMatrix

'# plot the second star (m2) position

TBGL_PushMatrix
TBGL_Translate(x2,y2,z2)
TBGL_Color(245, 230, 100)
TBGL_Sphere(Rad2)
TBGL_PopMatrix

'# plot the second star (m3) position

TBGL_PushMatrix
TBGL_Translate(x3,y3,z3)
TBGL_Color(245, 230, 200)
TBGL_Sphere(Rad3)
TBGL_PopMatrix
TBGL_DrawFrame

' ESCAPE key to disable callback
If TBGL_GetWindowKeyState(hWnd, %VK_ESCAPE) Then
TBGL_UnBindPeriodicFunction( hWnd )
Exit Sub
End If

End Sub

Sub orbits()
'# More complex due to 3 Body instead of simply 2 Body
'# interactions. THIS is the first half of the velocity
'# Calculations. Known As Leap Frog!
vx1 += 0.5*ax1*dt
vy1 += 0.5*ay1*dt
vz1 += 0.5*az1*dt
vx2 += 0.5*ax2*dt
vy2 += 0.5*ay2*dt
vz2 += 0.5*az2*dt
vx3 += 0.5*ax3*dt
vy3 += 0.5*ay3*dt
vz3 += 0.5*az3*dt
'# Calculate New positions
x1 += vx1*dt
y1 += vy1*dt
z1 += vz1*dt
x2 += vx2*dt
y2 += vy2*dt
z2 += vz2*dt
x3 += vx3*dt
y3 += vy3*dt
z3 += vz3*dt
'# Reset acceleration components To zero.
'# THIS is important!
ax1 = 0.0
ay1 = 0.0
az1 = 0.0
ax2 = 0.0
ay2 = 0.0
az2 = 0.0
ax3 = 0.0
ay3 = 0.0
az3 = 0.0
'# Recalculate r**3 denominators
r12 = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)
r312 = r12*Sqr(r12) + 0.01
r13 = (x1-x3)*(x1-x3) + (y1-y3)*(y1-y3) + (z1-z3)*(z1-z3)
r313 = r13*Sqr(r13) + 0.01
r23 = (x2-x3)*(x2-x3) + (y2-y3)*(y2-y3) + (z2-z3)*(z2-z3)
r323 = r23*Sqr(r23) + 0.01
'# Calculate acceleration components from each body.
'# We Add Or accumulate the acceleration components provided
'# by each of the other two stars To arrive At ONE resultant
'# Acceleration. We avoid self-gravity!
ax1 += -G*(x1-x2)*m2/r312
ax1 += -G*(x1-x3)*m3/r313
ay1 += -G*(y1-y2)*m2/r312
ay1 += -G*(y1-y3)*m3/r313
az1 += -G*(z1-z2)*m2/r312
az1 += -G*(z1-z3)*m3/r313
ax2 += -G*(x2-x1)*m1/r312
ax2 += -G*(x2-x3)*m3/r323
ay2 += -G*(y2-y1)*m1/r312
ay2 += -G*(y2-y3)*m3/r323
az2 += -G*(z2-z1)*m1/r312
az2 += -G*(z2-z3)*m3/r323
ax3 += -G*(x3-x2)*m2/r323
ax3 += -G*(x3-x1)*m1/r313
ay3 += -G*(y3-y2)*m2/r323
ay3 += -G*(y3-y1)*m1/r313
az3 += -G*(z3-z2)*m2/r323
az3 += -G*(z3-z1)*m1/r313
'# Calculate the second half of the velocity components
vx1 += 0.5*ax1*dt
vy1 += 0.5*ay1*dt
vz1 += 0.5*az1*dt
vx2 += 0.5*ax2*dt
vy2 += 0.5*ay2*dt
vz2 += 0.5*az2*dt
vx3 += 0.5*ax3*dt
vy3 += 0.5*ay3*dt
vz3 += 0.5*az3*dt
End Sub

[/code]