'Program FIG18.BAS calculates symbolic dynamics
'by J. C. Sprott <http://sprott.physics.wisc.edu/>
'Compile with the PowerBASIC Console Compiler <http://www.powerbasic.com/>
DEFEXT a-z      'do all calculations in extended (80-bit) precision

FUNCTION PBMAIN()
n&=3            'number of nonlinear odes
h=0.003         'stepsize of integration
pi=4#*ATN(1##)
CONSOLE NAME "Fig 18"
CONSOLE SET LOC 646,0
CURSOR OFF
CLS
GRAPHIC WINDOW "Fig 18", 0, 0, 640, 640 TO hWin1&
GRAPHIC ATTACH hWin1&, 0
GRAPHIC WIDTH 1
GRAPHIC COLOR RGB(0,0,0), RGB(255,255,255)
GRAPHIC CLEAR
GRAPHIC ATTACH hWin1&, 0
GRAPHIC BOX(0,0)-(640,640)
DIM x(n&), dxdt(n&), xnew(n&), dxdtold(n&), p&(2^n&,2^n&), m&(2^n&)
OPEN"fig18.dat" FOR OUTPUT AS #1
RANDOMIZE TIMER
x(1)=0.1-0.2*RND: x(2)=0.1-0.2*RND: x(3)=0.1-0.2*RND
'x(1)=0.1: x(2)=0.2: x(3)=0.3           'first chamber start
'x(1)=0: x(2)=0: x(3)=pi/2##            'periodic orbit
'x(1)=0: x(2)=pi/2##: x(3)=pi           'another periodic orbit
WHILE INKEY$<>CHR$(27)
    CALL rk4(x(),dxdt(),n&,h,xnew())
    a$="": s&=0
    FOR i&=1 TO 3
        IF xnew(i&)>=2##*pi THEN xnew(i&)=xnew(i&)-2##*pi
        IF xnew(i&)<0 THEN xnew(i&)=xnew(i&)+2##*pi
        IF dxdt(i&)*dxdtold(i&)<0 THEN  'switched chambers
            FOR j&=1 TO n&
                s&=s&-(dxdt(j&)<0)*2^(j&-1)
            NEXT j&
            a$=CHR$(65+s&)
        END IF
        dxdtold(i&)=dxdt(i&)
        x(i&)=xnew(i&)
    NEXT i&
    IF a$<>"" THEN
        'a$=chr$(int(rnd*2^n&)+65)         'randomness test
        j&=ASC(a$)-64
        i&=ASC(aold$)-64
        'print a$;
        PRINT#1,a$;
        IF i&>0 THEN
            INCR m&
            INCR m&(i&)
            INCR p&(i&,j&)
            FOR jj&=1 TO 2^n&
                LOCATE 2*jj&-1,6*i&-5
                PRINT CINT(1000*p&(i&,jj&)/m&(i&));
            NEXT jj&
            LOCATE 2*2^n&+3,6*i&-5
            PRINT CINT(1000*m&(i&)/m&);
        END IF
        aold$=a$
        x=x/4: y=y/2
        SELECT CASE a$
            CASE"B": y=y+1/2
            CASE"C": x=x+1/4
            CASE"D": x=x+1/4: y=y+1/2
            CASE"E": x=x+2/4
            CASE"F": x=x+2/4: y=y+1/2
            CASE"G": x=x+3/4
            CASE"H": x=x+3/4: y=y+1/2
        END SELECT
        IF m&>20 THEN GRAPHIC SET PIXEL(640*x,640*y)
    END IF
WEND
GRAPHIC SAVE"fig18.bmp"
WAITKEY$
END FUNCTION

SUB DERIVS (x(), dxdt(), n&)
    'Returns the time derivatives dxdt(i&) of x(i&)
    FOR i&=1 TO n&
        i1&=1+i& MOD n&
        dxdt(i&)=SIN(x(i1&))
    NEXT i&
END SUB

SUB RK4 (X(), DXDT(), N&, H, XNEW())
    'Fourth-order Runge-Kutta integrator
    DIM XT(N&), DXT(N&), DXM(N&)
    HH = H * .5##
    H6 = H / 6##
    CALL DERIVS(X(), DXDT(), N&)
    FOR I& = 1 TO N&
        XT(I&) = X(I&) + HH * DXDT(I&)
    NEXT I&
    CALL DERIVS(XT(), DXT(), N&)
    FOR I& = 1 TO N&
        XT(I&) = X(I&) + HH * DXT(I&)
    NEXT I&
    CALL DERIVS(XT(), DXM(), N&)
    FOR I& = 1 TO N&
        XT(I&) = X(I&) + H * DXM(I&)
        DXM(I&) = DXT(I&) + DXM(I&)
    NEXT I&
    CALL DERIVS(XT(), DXT(), N&)
    FOR I& = 1 TO N&
        XNEW(I&) = X(I&) + H6 * (DXDT(I&) + DXT(I&) + 2## * DXM(I&))
    NEXT I&
    ERASE DXM, DXT, XT
END SUB