11-16-2023, 06:28 PM
Code: (Select All)
100 Rem ITERATED FUNCTION SYSTEMS Ver 1.0 (c) 1993 by J. C. Sprott
110 DefDbl A-Z: Dim A(12)
120 Randomize Timer 'Reseed random number generator
130 Screen 12 'Assume VGA graphics
140 GoSub 300 'Set parameters
150 GoSub 400 'Iterate equations
160 GoSub 500 'Display results
170 GoSub 600 'Test results
180 On T% GOTO 130, 140, 150
190 End
300 Rem Set parameters
310 X = .05: Y = .05: XE = X + .00001: YE = Y
320 For I% = 1 To 12: A(I%) = .1 * (Int(25 * Rnd) - 12): Next I%
330 J0 = Abs(A(1) * A(4) - A(2) * A(3))
340 J1 = Abs(A(7) * A(10) - A(8) * A(9))
350 If J0 + J1 = 0 Or J0 > 1 Or J1 > 1 Then GoTo 320 'Not contracting
360 P = J0 / (J0 + J1)
370 T% = 3: LSUM = 0: N = 1: N1 = 0: N2 = 0
380 XMIN = 1000000#: XMAX = -XMIN: YMIN = XMIN: YMAX = XMAX
390 Return
400 Rem Iterate equations
410 If X <> XE Then If Rnd > P Then R% = 6 Else R% = 0
420 XNEW = A(1 + R%) * X + A(2 + R%) * Y + A(5 + R%)
430 YNEW = A(3 + R%) * X + A(4 + R%) * Y + A(6 + R%)
490 Return
500 Rem Display results
510 If N < 100 Or N > 1000 Then GoTo 560
520 If X < XMIN Then XMIN = X
530 If X > XMAX Then XMAX = X
540 If Y < YMIN Then YMIN = Y
550 If Y > YMAX Then YMAX = Y
560 If N = 1000 Then GoSub 800
570 If X > XL And X < XH And Y > YL And Y < YH And N > 1000 Then PSet (X, Y)
590 Return
600 Rem Test results
610 GoSub 700 'Calculate Lyapunov exponent
620 GoSub 900 'Calculate correlation dimension
630 If N > 21000 Then T% = 2 'Candidate IFS found
640 If Abs(XNEW) + Abs(YNEW) > 1000000# Then T% = 2 'Unbounded
650 If N > 998 And (F < 1 Or L > -.2) Then T% = 2 'Uninteresting
660 If Len(InKey$) Then T% = 0 'User key press
670 X = XNEW: Y = YNEW: N = N + 1
690 Return
700 Rem Calculate Lyapunov exponent
710 XSAVE = XNEW: YSAVE = YNEW: X = XE: Y = YE
720 GoSub 400 'Reiterate equations
730 DLX = XNEW - XSAVE: DLY = YNEW - YSAVE: DL2 = DLX * DLX + DLY * DLY
740 DF = 10000000000# * DL2: RS = 1# / Sqr(DF)
750 XE = XSAVE + RS * (XNEW - XSAVE): YE = YSAVE + RS * (YNEW - YSAVE)
760 XNEW = XSAVE: YNEW = YSAVE
770 LSUM = LSUM + Log(DF): L = .721347 * LSUM / N
790 Return
800 Rem Resize the screen (and discard the first thousand iterates)
810 DX = .1 * (XMAX - XMIN): DY = .1 * (YMAX - YMIN)
820 XL = XMIN - DX: XH = XMAX + DX: YL = YMIN - DY: YH = YMAX + DY
830 If XH - XL < .000001 Or YH - YL < .000001 Then GoTo 890
840 Window (XL, YL)-(XH, YH): Cls
850 Line (XL, YL)-(XH, YH), , B
890 Return
900 Rem Calculate fractal dimension
910 If N < 200 Then GoTo 990 'Wait for transient to settle
920 If N = 200 Then D2MAX = (XMAX - XMIN) ^ 2 + (YMAX - YMIN) ^ 2
930 If N = 200 Or Rnd < .02 Then XS = X: YS = Y 'New reference point
940 DX = XNEW - XS: DY = YNEW - YS
950 D2 = DX * DX + DY * DY
960 If D2 < .001 * D2MAX Then N2 = N2 + 1
970 If D2 < .00001 * D2MAX Then N1 = N1 + 1: F = .434294 * Log(N2 / N1)
990 ReturnMore ancient era code off the hard drives, as I copy and clean things up.
Methinks this one could use a delay of some sort in it, but I'll be danged if I could figure out where to place one so that it doesn't just slow the whole program down to uselessness!
