'Program SEARCH.BAS searches for periodic animated attractors
'Compile with QuickBASIC, VisualBASIC for DOS, or PowerBASIC
'(c) 1996 by J. C. Sprott (sprott@jun.physics.wisc.edu)

DECLARE SUB advancexy (x#, y#, z#, n%)
DECLARE SUB display (x#, y#, z#, n%)
DECLARE SUB testsoln (x#, y#, z#, n%)
DECLARE SUB scn2bmp (filename$)
DECLARE SUB setparams (x#, y#, z#)
DECLARE SUB lyapunov (x#, y#, z#, n%, l#)
DECLARE FUNCTION mkhex$ (a$)

DEFDBL A-Z                                  'Use double precision
DIM a(6)                                    'Array of coefficients
RANDOMIZE TIMER                             'Reseed random numbers
SCREEN 12                                   'Assume VGA graphics
n% = 0
WHILE INKEY$ = ""                           'Loop until key is pressed
	IF n% = 0 THEN CALL setparams(x, y, z)
	CALL advancexy(x, y, z, n%)             'Advance the solution
	CALL display(x, y, z, n%)               'Display the results
	CALL testsoln(x, y, z, n%)              'Test the solution
WEND
END

SUB advancexy (x, y, z, n%)
	SHARED a()
	xnew = a(1) + a(2) * x + a(3) * x * x + a(4) * y + a(5) * z + a(6) * SIN(6.283185307# * n% / 16)
	ynew = x
	znew = y
	x = xnew: y = ynew: z = znew
	n% = n% + 1
END SUB

SUB display (x, y, z, n%)
	SHARED code$
	STATIC xmin, xmax, ymin, ymax, zmin, zmax, w, h, pix%
	SELECT CASE n%
		CASE 1                                  'Initialize limits
			xmin = 1000: ymin = xmin: zmin = xmin
			xmax = -xmin: ymax = -ymin: zmax = -zmin
			pix% = 0
		CASE 2 TO 99                            'Skip these
		CASE 100 TO 999                         'Update limits
			IF x < xmin THEN xmin = x
			IF x > xmax THEN xmax = x
			IF y < ymin THEN ymin = y
			IF y > ymax THEN ymax = y
			IF z < zmin THEN zmin = z
			IF z > zmax THEN zmax = z
		CASE 1000                               'Clear the screen
			CLS
            FOR i% = 1 to 3
				LINE(0, 120 * i%) - (639, 120 * i%)
                LINE(160 * i%, 0) - (160 * i%, 479)
            NEXT i%
			IF CSNG(xmax) = CSNG(xmin) THEN xmax = xmin + 1
			IF CSNG(ymax) = CSNG(ymin) THEN ymax = ymin + 1
			dx = (xmax - xmin) / 10: xmin = xmin - dx: xmax = xmax + dx
			dy = (ymax - ymin) / 10: ymin = ymin - dy: ymax = ymax + dy
			w = 160 / (xmax - xmin): h = 120 / (ymin - ymax)
		CASE ELSE                               'Plot data
			IF x < 2 * xmin OR x > 2 * xmax THEN EXIT SUB
			xo% = 160 * (n% MOD 4)
			yo% = 120 * (n% \ 4 MOD 4)
			xp% = w * (x - xmin) + xo%
			yp% = h * (y - ymax) + yo%
			zp% = 1 + 15 * (z - zmin) / (zmax - zmin)
			IF POINT(xp%, yp%) = 0 THEN pix% = pix% + 1
			PSET (xp%, yp%), zp%                 'Illuminate screen pixel
			IF pix% > 10000 THEN PRINT code$: n% = 0: CALL scn2bmp(code$ + ".BMP")
	END SELECT
END SUB

SUB lyapunov (x, y, z, n%, l)
	STATIC xe, ye, ze, lsum
	IF n% = 1 THEN lsum = 0: xe = .000001#: ye = 0: ze = 0
	xsave = x: ysave = y: zsave = z: x = xe: y = ye: z = ze
	n% = n% - 1
	CALL advancexy(x, y, z, n%)                 'Reiterate equations
	dx = x - xsave: dy = y - ysave: dz = z - zsave
	d2 = dx * dx + dy * dy + dz * dz
	df = 100000000000# * d2: rs = 1# / SQR(df)
	xe = xsave + rs * (x - xsave): x = xsave
	ye = ysave + rs * (y - ysave): y = ysave
	ze = zsave + rs * (z - zsave): z = zsave
	lsum = lsum + LOG(df): l = .721347 * lsum / (n% + 1)
END SUB

FUNCTION mkhex$ (a$)
	'Converts a sequence of space-delimited hexadecimal bytes to a string
	'For example, mkhex$("4A 4B 4C") evaluates to "JKL"
	bytes% = (LEN(a$) + 1) / 3
	c$ = ""
	FOR i% = 1 TO bytes%
		b$ = "&H" + MID$(a$, 3 * i% - 2, 3)
		c$ = c$ + CHR$(VAL(b$))
	NEXT
	mkhex$ = c$
END FUNCTION

SUB scn2bmp (filename$)
	'Captures the screen to a file filename$ in bitmaped (BMP) format
	'Palette correct in EGA and VGA color modes 7 - 12 only
	x1% = 0: y1% = 0: x2% = 639: y2% = 479      'Corners of image
	x2% = x1% + 16 * INT(1 + (x2% - x1%) \ 16 + .01) - 1
	BEEP
	sw% = 1 + x2% - x1%         'Width of image in pixels
	sh% = 1 + y2% - y1%         'Height of image in pixels
	sf& = sw% * CLNG(sh%) \ 2   'Bytes in image
	f% = FREEFILE
	OPEN filename$ FOR OUTPUT AS f%
	PRINT #f%, mkhex$("42 4D");         'ASCII "BM"
	PRINT #f%, MKL$(sf& + 118);         'Size of file in bytes
	PRINT #f%, mkhex$("00 00 00 00");   'Must be zero
	PRINT #f%, mkhex$("76 00 00 00 28 00 00 00");
	PRINT #f%, MKL$(sw%);               'Screen width
	PRINT #f%, MKL$(sh%);               'Screen height
	PRINT #f%, mkhex$("01 00 04 00 00 00 00 00");
	PRINT #f%, MKL$(sf&);               'Size of image in bytes
	PRINT #f%, mkhex$("00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00");
	PRINT #f%, mkhex$("00 00 00 00");   'Palette RGB values (0 - 15)
	PRINT #f%, mkhex$("AA 00 00 00");
	PRINT #f%, mkhex$("00 AA 00 00");
	PRINT #f%, mkhex$("AA AA 00 00");
	PRINT #f%, mkhex$("00 00 AA 00");
	PRINT #f%, mkhex$("AA 00 AA 00");
	PRINT #f%, mkhex$("00 55 AA 00");
	PRINT #f%, mkhex$("AA AA AA 00");
	PRINT #f%, mkhex$("FF FF FF 00");
	PRINT #f%, mkhex$("FF 55 55 00");
	PRINT #f%, mkhex$("55 FF 55 00");
	PRINT #f%, mkhex$("FF FF 55 00");
	PRINT #f%, mkhex$("55 55 FF 00");
	PRINT #f%, mkhex$("FF 55 FF 00");
	PRINT #f%, mkhex$("55 FF FF 00");
	PRINT #f%, mkhex$("FF FF FF 00");
	a$ = SPACE$(sw% \ 2)
	FOR y% = y2% TO y1% STEP -1
		x% = x1%
		WHILE x% < x2%
			a% = POINT(x%, y%)
			a% = 16 * a%
			x% = x% + 1
			a% = a% OR POINT(x%, y%)
			x% = x% + 1
			MID$(a$, (x% - x1%) \ 2, 1) = CHR$(a%)
		WEND
		PRINT #f%, a$;
	NEXT y%
	CLOSE f%
END SUB

SUB setparams (x, y, z)
	SHARED a(), code$
	x = 0: y = 0: z = 0
	code$ = ""
	FOR i% = 1 TO 6
		a(i%) = (INT(25 * RND) - 12) / 10#
		IF a(6) = 0 THEN a(6) = 1               'Don't allow zero drive
		code$ = code$ + CHR$(77 + CINT(10 * a(i%)))
	NEXT i%
END SUB

SUB testsoln (x, y, z, n%)
	nmax% = 32000                               'Bailout value
	CALL lyapunov(x, y, z, n%, l)               'Get Lyapunov exponent (l)
	IF n% = nmax% THEN n% = 0                   'Bailout value reached
	IF n% > 100 AND l < .005 THEN n% = 0        'Solution is not chaotic
	IF ABS(x) + ABS(y) + ABS(z) > 100000# THEN n% = 0   'Solution is unbounded
END SUB