Example: Section 8.3 (page 265), A little MORE Gravity! 2 Bodies

[Michael Clease]

[code=thinbasic]
' Example: Section 8.3 (page 265), A little MORE Gravity! 2 Bodies

' 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, r2, r3, ax1, ay1, az1, dt As Double
Global vx2, vy2, vz2, x2, y2, z2, ax2, ay2, az2, G As Double
Global m1,m2, rad1 ,rad2 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()
'# initial x,y,z positions For both stars
x1 = 1.0
y1 = 0.0
z1 = 0.0
x2 = -1.0
y2 = 0.0
z2 = 0.0
'# initial vx,vy,vz velocities For both stars
vx1 = 0.0
vy1 = -0.128571428
vz1 = 0.0
vx2 = 0.0
vy2 = 0.3
vz2 = 0.0
'# initial masses For both stars
m1 = 0.7
m2 = 0.3
rad1 = 0.1*m1
rad2 = 0.1*m2
'# arbitrary "Big G" gravitational constant
G = 1.0
'# calculate distance And r**3 denominator For universal
'# gravitation
r2 = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)
r3 = r2*Sqr(r2)
'# calculate acceleration components along x,y,z axes
'# First For m1
ax1 = -G*(x1-x2)*m2/r3
ay1 = -G*(y1-y2)*m2/r3
az1 = -G*(z1-z2)*m2/r3
'# now For m2
ax2 = -G*(x2-x1)*m1/r3
ay2 = -G*(y2-y1)*m1/r3
az2 = -G*(z2-z1)*m1/r3
'# THIS Value keeps a smooth orbit On my workstation
'# Smaller values slow down orbit, higher values speed up orbit
dt = 0.04
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, 3.0, 0.0, 0.0, 0.0)


TBGL_PushMatrix
TBGL_Translate(x1,y1,z1)
TBGL_Color(245, 230, 100)
TBGL_Sphere(Rad1)
TBGL_PopMatrix

'# plot the second star (m2) position

TBGL_PushMatrix
TBGL_Translate(x2,y2,z2)
TBGL_Color(245, 150, 100)
TBGL_Sphere(Rad2)
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()
'# calculate front half of velocity vector 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
'# calculate x,y,z positions For both stars
x1 += vx1*dt
y1 += vy1*dt
z1 += vz1*dt
x2 += vx2*dt
y2 += vy2*dt
z2 += vz2*dt
'# calculate the New r**3 denominator For each star position
r2 = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)
r3 = r2*Sqr(r2)
'# calculate the New acceleration components
ax1 = -G*(x1-x2)*m2/r3
ay1 = -G*(y1-y2)*m2/r3
az1 = -G*(z1-z2)*m2/r3
ax2 = -G*(x2-x1)*m1/r3
ay2 = -G*(y2-y1)*m1/r3
az2 = -G*(z2-z1)*m1/r3
'# calculate the back half 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
'#Send calculated x,y,z star positions To the display
'glutPostRedisplay()
End Sub

[/code]

[kryton9]

Again thanks Michael, fun to see these conversions and what is coming up in the book!!