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