'
' *****************************************************************
'
'                       FROM CALCULUS TO CHAOS
'      
'                    An Introduction to Dynamics
'
'                 
'                                by
'
'                           David Acheson
'
' ******************************************************************
'
'        The program which follows is part of the above book,
'            published in 1997 by Oxford University Press.
'                   ISBN 0-19-850077-7 (paperback)
'             
'
'                   Copyright ¸ David Acheson 1997
'
'
'
'
'                      Program 9:  THREEBP.BAS
'
'
 DEFDBL A-H, K-M, O-Z: DEFINT I-J, N
 n = 12: OPTION BASE 1
 DIM x(n), xc(n), f(n), c1(n), c2(n), c3(n), c4(n)

REM ****** Setting up graphics ******
 
  CLS : SCREEN 9
  PAINT (1, 1), 9
   xm = .75: ym = xm
  VIEW (180, 17)-(595, 330), 0, 13
  WINDOW (-xm, -ym)-(xm, ym)
  LINE (-xm, 0)-(xm, 0), 8: LINE (0, -ym)-(0, ym), 8
  LOCATE 13, 76: PRINT xm

REM ****** Step-by-step method ******

  m1 = .5#: m2 = .5#: m3 = .5#
  LOCATE 13, 1: INPUT "x3,y3"; x3, y3

  t = 0#
   
     REM ** initial x-coordinates **
   xc(1) = -.5#: xc(2) = .5#: xc(3) = x3
     REM ** initial y-coordinates **
   xc(4) = 0#: xc(5) = 0#: xc(6) = y3
     REM ** initial x-velocities **
   xc(7) = 0: xc(8) = 0#: xc(9) = 0#
     REM ** initial y-velocities **
   xc(10) = -.3#: xc(11) = .3#: xc(12) = -.3#
     
  h = .003#

  DO
    GOSUB Runge
    PSET (xc(1), xc(4)), 12
    PSET (xc(2), xc(5)), 9
    PSET (xc(3), xc(6)), 15
    t = t + h
      
     REM ** Energy conserved? **
     kin1 = .5# * m1 * (xc(7) ^ 2# + xc(10) ^ 2#)
     kin2 = .5# * m2 * (xc(8) ^ 2# + xc(11) ^ 2#)
     kin3 = .5# * m3 * (xc(9) ^ 2# + xc(12) ^ 2#)
     pot = -(m1 * m2 / r12 + m2 * m3 / r23 + m3 * m1 / r31)
     energy = kin1 + kin2 + kin3 + pot
     LOCATE 21, 1: PRINT "Energy"
     LOCATE 22, 1: PRINT energy
       
    LOCATE 17, 7: PRINT "             "
    LOCATE 17, 1: PRINT "Time ="; CSNG(t)
    
  LOOP UNTIL INKEY$ = "q"

END

REM ****** Subroutines ******

Equations:
    
 d21 = x(2) - x(1): d32 = x(3) - x(2): d13 = x(1) - x(3)
 d54 = x(5) - x(4): d65 = x(6) - x(5): d46 = x(4) - x(6)
   
 r12 = (d21 ^ 2# + d54 ^ 2#) ^ .5#
 r23 = (d32 ^ 2# + d65 ^ 2#) ^ .5#
 r31 = (d13 ^ 2# + d46 ^ 2#) ^ .5#
 
 p12 = r12 ^ 3#: p23 = r23 ^ 3#: p31 = r31 ^ 3#
   
 f(1) = x(7): f(2) = x(8): f(3) = x(9)
 f(4) = x(10): f(5) = x(11): f(6) = x(12)
   
 f(7) = m2 * d21 / p12 - m3 * d13 / p31
 f(8) = m3 * d32 / p23 - m1 * d21 / p12
 f(9) = m1 * d13 / p31 - m2 * d32 / p23
 f(10) = m2 * d54 / p12 - m3 * d46 / p31
 f(11) = m3 * d65 / p23 - m1 * d54 / p12
 f(12) = m1 * d46 / p31 - m2 * d65 / p23
   
 RETURN

Runge:
     
    FOR i = 1 TO n: x(i) = xc(i): NEXT
    GOSUB Equations
    FOR i = 1 TO n: c1(i) = h * f(i): NEXT
   
    FOR i = 1 TO n: x(i) = xc(i) + c1(i) / 2#: NEXT
    GOSUB Equations
    FOR i = 1 TO n: c2(i) = h * f(i): NEXT
   
    FOR i = 1 TO n: x(i) = xc(i) + c2(i) / 2#: NEXT
    GOSUB Equations
    FOR i = 1 TO n: c3(i) = h * f(i): NEXT
   
    FOR i = 1 TO n: x(i) = xc(i) + c3(i): NEXT
    GOSUB Equations
    FOR i = 1 TO n: c4(i) = h * f(i): NEXT
     
    FOR i = 1 TO n
    xc(i) = xc(i) + (c1(i) + 2# * c2(i) + 2# * c3(i) + c4(i)) / 6#
    NEXT
   
    RETURN