REM STRANGE ATTRACTOR ICON SEARCH (c) 1995 by J. C. Sprott
REM Run under QBASIC, QuickBASIC, or VB/DOS
DEFDBL A-Z 'Use double precision
DIM A(12), PAL&(255)
NMAX = 61000 'Maximum number of iterations
TWOPI = 8 * ATN(1) 'A useful constant (2 pi)
RANDOMIZE TIMER 'Reseed random number generator
SCREEN 13 'Set graphics mode
SW% = 320 'Screen width
SH% = 200 'Screen height
NC% = 254 'Number of colors
WHILE T% = 0
GOSUB Init 'Initialize
WHILE T% = 1
GOSUB Iteqn 'Iterate equations
GOSUB Display 'Display results
GOSUB Test 'Test results
WEND
WEND
PALETTE: CLS
END
Init: 'Initialize
X = .05: Y = .05: Z = .05 'Initial condition
XE = X + .000001: YE = Y: ZE = Z
T% = 1: LSUM = 0: N = 0: NL = 0
XL = 1000000!: XH = -XL: YL = XL: YH = XH: ZL = XL: ZH = XH
NS% = 2 + INT(8 * RND) 'Number of sectors (2 to 9)
OF = .5 + 1.5 * RND 'Overlap factor (.5 to 2)
FOR I% = 1 TO 12 'Get random coefficients
A(I%) = (INT(25 * RND) - 12) / 10
NEXT I%
RETURN
Iteqn: 'Iterate equations
XNEW = A(1) + X * (A(2) + A(3) * X + A(4) * Y) + Y * (A(5) + A(6) * Y)
YNEW = A(7) + X * (A(8) + A(9) * X + A(10) * Y) + Y * (A(11) + A(12) * Y)
ZNEW = X * X + Y * Y
N = N + 1
RETURN
Display: 'Display results
IF N > 100 AND N < 1000 THEN 'Get scale limits for graph
IF X < XL THEN XL = X ELSE IF X > XH THEN XH = X
IF Y < YL THEN YL = Y ELSE IF Y > YH THEN YH = Y
IF Z < ZL THEN ZL = Z ELSE IF Z > ZH THEN ZH = Z
END IF
IF N = 1000 THEN 'Set palette and clear screen
RANB = RND: CB% = 1 + INT(3 * RND) 'Blue phase and period
RANG = RND: CG% = 1 + INT(3 * RND) 'Green phase and period
RANR = RND: CR% = 1 + INT(3 * RND) 'Red phase and period
BC% = INT(2 + NC% * RND) 'Choose random background color
FOR I% = 2 TO NC% + 1 'Redefine palette colors
B% = INT(32 + 32 * SIN(CB% * TWOPI * (I% / NC% + RANB)))
G% = INT(32 + 32 * SIN(CG% * TWOPI * (I% / NC% + RANG)))
R% = INT(32 + 32 * SIN(CR% * TWOPI * (I% / NC% + RANR)))
PAL&(I%) = 65536 * B% + 256 * G% + R%
IF I% = BC% THEN 'Set background and shadow colors
PAL&(0) = 65536 * INT(B% / 2) + 256 * INT(G% / 2) + INT(R% / 2)
PAL&(1) = 65536 * INT(B% / 3) + 256 * INT(G% / 3) + INT(R% / 3)
END IF
NEXT I%
CLS : PALETTE USING PAL&(0)
XZ = .05 * SW% / (ZH - ZL)
YZ = .05 * SH% / (ZH - ZL)
END IF
IF N > 1000 THEN 'Plot point on screen
IF XH > XL THEN XP = SW% * (X - XL) / (XH - XL)
IF YH > YL THEN YP = SH% * (YH - Y) / (YH - YL)
S% = INT(NS% * RND) 'Choose sector randomly
TH = TWOPI * (OF * YP / SH% + S%) / NS%
IF (NS% MOD 2) = 0 AND (S% MOD 2) = 0 THEN TH = TWOPI / NS% - TH
YP = .5 * SH% * (1 + XP * COS(TH) / SW%)
XP = .5 * SW% * (1 + XP * SIN(TH) / SW%)
C% = 2 + INT(NC% * (Z - ZL) / (ZH - ZL) + NC%) MOD NC%
IF C% > POINT(XP, YP) THEN PSET (XP, YP), C%
XP = XP + XZ * (Z - ZL): YP = YP + YZ * (Z - ZL)
IF POINT(XP, YP) = 0 THEN PSET (XP, YP), 1
END IF
RETURN
Test: 'Test results
IF ABS(XNEW) + ABS(YNEW) + ABS(ZNEW) > 1000000! THEN T% = 0 'Unbounded
XSAVE = XNEW: YSAVE = YNEW: ZSAVE = ZNEW
X = XE: Y = YE: Z = ZE: N = N - 1
GOSUB Iteqn 'Reiterate equations
DLX = XNEW - XSAVE: DLY = YNEW - YSAVE: DLZ = ZNEW - ZSAVE
DL2 = DLX * DLX + DLY * DLY + DLZ * DLZ
IF CSNG(DL2) > 0 THEN 'Don't divide by zero
DF = 1000000000000# * DL2
RS = 1 / SQR(DF)
XE = XSAVE + RS * (XNEW - XSAVE): XNEW = XSAVE
YE = YSAVE + RS * (YNEW - YSAVE): YNEW = YSAVE
ZE = ZSAVE + RS * (ZNEW - ZSAVE): ZNEW = ZSAVE
LSUM = LSUM + LOG(DF): NL = NL + 1
L = .721347 * LSUM / NL 'This is the Lyapunov exponent
END IF
IF N > 100 AND L < .005 THEN T% = 0 'Not chaotic
IF N > NMAX THEN T% = 0 'Strange attractor found
IF LEN(INKEY$) THEN T% = 2 'Exit on keypress
X = XNEW: Y = YNEW: Z = ZNEW 'Update value of variables
RETURN