'
' *****************************************************************
'
'                       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 7:  PENDANIM.BAS
'
'
DEFDBL A-H, K-M, O-Z: DEFINT I-J, N
n = 2: 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
  VIEW (180, 17)-(595, 330), 0, 14
  WINDOW (-2, -2)-(2, 2)
   LOCATE 22, 1: PRINT "Time"
  
REM ****** Step-by-step method ******

DO
 
    k = .1#
    t = 0#
    LOCATE 14, 1: INPUT "angle"; xc(1)
    LOCATE 15, 1: INPUT "angvel"; xc(2)
    h = .05#: animsteps = 4
     
      DO
          CLS
           b1 = SIN(xc(1)): b2 = -COS(xc(1))
           LINE (0, 0)-(b1, b2), 4
           CIRCLE (b1, b2), .05, 9: PAINT (b1, b2), 9
            FOR j = 1 TO animsteps
              GOSUB Runge
              t = t + h
            NEXT
          LOCATE 22, 6: PRINT CSNG(t)
      LOOP UNTIL INKEY$ = "q"

LOOP

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

Equations:
   
 f(1) = x(2)
 f(2) = -k * x(2) - SIN(x(1))
 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