Results 1 to 3 of 3

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

  1. #1
    thinBasic MVPs
    Join Date
    May 2007
    Location
    UK
    Posts
    1,427
    Rep Power
    159

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

    [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]
    Home Desktop : Windows 7 - Intel Pentium (D) - 3.0 Ghz - 2GB - Geforce 6800GS
    Home Laptop : WinXP Pro SP3 - Intel Centrino Duo - 1.73 Ghz - 2 GB - Intel GMA 950
    Home Laptop : Windows 10 - Intel(R) Core(TM) i5-4210U CPU @ 1.70GHz, 2401 Mhz, 2 Core(s), 4 Logical Processor(s) - 4 GB - Intel HD 4400
    Work Desktop : Windows 10 - Intel I7 - 4 Ghz - 8GB - Quadro Fx 370

  2. #2
    thinBasic MVPs kryton9's Avatar
    Join Date
    Nov 2006
    Location
    Naples, Florida & Duluth, Georgia
    Age
    67
    Posts
    3,869
    Rep Power
    404

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

    This is a really sweet example. I always wanted to see how all of this was calculated. I am still back early on in the book. I am going in order, but it is so nice to see what I will be reading about. I am really excited!!

    This is great, getting a refresher course in math and having fun and enjoying it at the same time, that is hard to beat. I gave up gaming night to have fun with this tonight, so you know how much I am enjoying it
    Acer Notebook: Win 10 Home 64 Bit, Core i7-4702MQ @ 2.2Ghz, 12 GB RAM, nVidia GTX 760M and Intel HD 4600
    Raspberry Pi 3: Raspbian OS use for Home Samba Server and Test HTTP Server

  3. #3
    Member sblank's Avatar
    Join Date
    Feb 2010
    Location
    Wayne City, Illinois
    Posts
    80
    Rep Power
    23

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

    Hi Michael!

    Very, very nice work!

    Thanks!

    Cheers,

    Stan

Similar Threads

  1. Example: Section 5.4 (page 41), Combined with Exercise 1 (page 44)
    By kryton9 in forum ThinBASIC programming in OpenGL/TBGL
    Replies: 2
    Last Post: 26-02-2010, 20:46

Members who have read this thread: 1

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •