'Program SIMANIM.BAS (c) 1996 by J. C. Sprott <sprott@juno.physics.wisc.edu>
'Produces 24 frames of the 19 simplest know examples of chaotic flows
'Ref: http://sprott.physics.wisc.edu/paper212.htm
'Compile this program with PowerBASIC

defext a-z      	'Use extended (80-bit) precision throughout
dim a(30) as shared ext
dim p%(15)
nd&=100000      	'Number of initial iterations to discard
nf%=24				'Number of frames
thn%=0				'Starting angle (0 to nf% -1)
screen 7
for i%=0 to 15: read p%(i%): next i%
data 17,24,0,32,48,56,38,46,6,14,39,20,28,62,55,63	'palette data
palette using p%(0)
c$=command$
call getcode(c$,code$)    		'Get code for case to analyze
call coeff(code$,a())           'Get equation coefficients in array a(30)
more:
h=.01           				'Iteration step size
x=.05: y=.05: z=.05				'Initial conditions
if c$="A" then x=0: y=5: z=0
xmin=1e37: ymin=xmin: zmin=xmin
xmax=-xmin: ymax=-ymin: zmax=-zmin
n&=0
do								'Loop until user presses the <Esc> key
	if inkey$=chr$(27) then cls: end
	if n&>nd&/5 and n&<nd& then call minmax(x,y,z,xmin,xmax,ymin,ymax,zmin,zmax)
	if n&=nd& then              'Set up the screen
    	xmm=xmax-xmin: xmax=xmax+.1*xmm: xmin=xmin-.1*xmm
        ymm=ymax-ymin: ymax=ymax+.1*ymm: ymin=ymin-.1*ymm
		zmm=zmax-zmin: zmax=zmax+.1*zmm: zmin=zmin-.1*zmm
		xmm=xmax-xmin
		ymm=ymax-ymin
		zmm=zmax-zmin
        if xmm>1.6*ymm then ymm=xmm/1.6
        if zmm>1.6*ymm then ymm=zmm/1.6
        xav=(xmax+xmin)/2
        yav=(ymax+ymin)/2
        zav=(zmax+zmin)/2
        h=h/10
        cls
	end if
	call rk4(x,y,z,h)           'Advance (x,y,z) by time step h
	if n&>nd& then              'Assume transient has settled
    	x1=(x-xav)/ymm			'Scales attractor to maximum y
        y1=(y-yav)/ymm
        z1=(z-zav)/ymm
		th=thn%*6.2832/nf%
		zp=18*(.5+z1*cos(th)+x1*sin(th))
        if zp<2 then zp=2
        if zp>15 then zp=15
        xp=160+200*(x1*cos(th)-z1*sin(th))/(1.5-zp/30)
		yp=100-200*y1/(1.5-zp/30)
		if point(xp,yp)<zp then pset(xp,yp),zp
        xp=xp+zp: yp=yp+zp
        if point(xp,yp)=0 then pset(xp,yp),1
	end if
	n&=n&+1
    if n&=1990000 then call scn2bmp("Case"+c$+serialno$(thn%,2)+".bmp"): exit loop
loop
incr thn%
print thn%: call morse(str$(thn%))
if thn%<nf% then goto more
print"Done";: call whoop: call whoop
thn%=0
end

sub rk4(x,y,z,h)        'Advance (x,y,z) with fourth order Runge-Kutta
k1x=h*fnx(x,y,z)
k1y=h*fny(x,y,z)
k1z=h*fnz(x,y,z)
k2x=h*fnx(x+.5*k1x,y+.5*k1y,z+.5*k1z)
k2y=h*fny(x+.5*k1x,y+.5*k1y,z+.5*k1z)
k2z=h*fnz(x+.5*k1x,y+.5*k1y,z+.5*k1z)
k3x=h*fnx(x+.5*k2x,y+.5*k2y,z+.5*k2z)
k3y=h*fny(x+.5*k2x,y+.5*k2y,z+.5*k2z)
k3z=h*fnz(x+.5*k2x,y+.5*k2y,z+.5*k2z)
k4x=h*fnx(x+k3x,y+k3y,z+k3z)
k4y=h*fny(x+k3x,y+k3y,z+k3z)
k4z=h*fnz(x+k3x,y+k3y,z+k3z)
x=x+(k1x+2*(k2x+k3x)+k4x)/6
y=y+(k1y+2*(k2y+k3y)+k4y)/6
z=z+(k1z+2*(k2z+k3z)+k4z)/6
end sub

def fnx(x,y,z)=a(1)+a(2)*x+a(3)*x*x+a(4)*x*y+a(5)*x*z+a(6)*y+a(7)*y*y+a(8)*y*z+a(9)*z+a(10)*z*z
def fny(x,y,z)=a(11)+a(12)*x+a(13)*x*x+a(14)*x*y+a(15)*x*z+a(16)*y+a(17)*y*y+a(18)*y*z+a(19)*z+a(20)*z*z
def fnz(x,y,z)=a(21)+a(22)*x+a(23)*x*x+a(24)*x*y+a(25)*x*z+a(26)*y+a(27)*y*y+a(28)*y*z+a(29)*z+a(30)*z*z

sub getcode(c$,code$)           'Get code for case to analyze
if len(c$)=31 then code$=c$: exit sub
again:
if c$="" then
	print"Choose case to analyze (A-S): ";
	while c$="": c$=inkey$: wend
    cls
end if
c$=ucase$(c$)
select case c$
	case"A": code$="QMMMMMWMMMMMCMMMMMWMMWMMMMMCMMM"
	case"B": code$="QMMMMMMMWMMMWMMMCMMMMWMMCMMMMMM"
	case"C": code$="QMMMMMMMWMMMWMMMCMMMMWMCMMMMMMM"
	case"D": code$="QMMMMMCMMMMMWMMMMMMWMMMMMWMkMMM"
	case"E": code$="QMMMMMMMWMMMMWMMCMMMMW%MMMMMMMM"
	case"F": code$="QMMMMMWMMWMMCMMMRMMMMMMWMMMMMCM"
	case"G": code$="QMQMMMMMMWMMMMMWCMMMMMCMMMWMMMM"
	case"H": code$="QMMMMMCMMMWMWMMMRMMMMMWMMMMMMCM"
	case"I": code$="QMMMMMKMMMMMWMMMMMMWMMWMMMMWMCM"
	case"J": code$="QMMMMMMMMaMMMMMM9MMWMMCMMMWWMMM"
	case"K": code$="QMMMWMMMMCMMWMMMCMMMMMWMMMMMMPM"
	case"L": code$="QMMMMMWMMtMMMVMMCMMMMWCMMMMMMMM"
	case"M": code$="QMMMMMMMMCMMMCMMCMMMM^^MMMWMMMM"
	case"N": code$="QMMMMM9MMMMMWMMMMMMMWWMMMMWMM9M"
	case"O": code$="QMMMMMWMMMMMWMMMMMMCMMWMMWhMMMM"
	case"P": code$="QMMMMMhMMWMMCMMMMWMMMMWMMMWMMMM"
	case"Q": code$="QMMMMMMMMCMMWMMMCMMMMMlMMMMWMRM"
	case"R": code$="QVMMMMCMMMMQMMMMMMMWMMMMWMMMMCM"
	case"S": code$="QMCMMM%MMMMMWMMMMMMMWWWMMMMMMMM"
	case else: c$="": beep: goto again
end select
end sub

sub coeff(code$,a(1))           'Return values of 30 equation coefficients
for i%=1 to 30
	a(i%)=(asc(mid$(code$,i%+1))-77)/10
next i%
end sub

sub minmax(x,y,z,xmin,xmax,ymin,ymax,zmin,zmax) 'Find min and max of (x,y,z)
xmin=min(x,xmin)
xmax=max(x,xmax)
ymin=min(y,ymin)
ymax=max(y,ymax)
zmin=min(z,zmin)
zmax=max(z,zmax)
end sub

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% = 319: y2% = 199		'Corners of image
    x2% = x1% + 16 * ceil((x2% - x1%) \ 16 + .01) - 1
    sw% = 1 + x2% - x1%			'Width of image in pixels
	sh% = 1 + y2% - y1%	 		'Height of image in pixels
    sf& = sw% * 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$("77 77 00 00");	'Palette BGR values (0 - 15)
	PRINT #f%, mkhex$("44 44 00 00");
	PRINT #f%, mkhex$("00 00 22 00");
	PRINT #f%, mkhex$("00 00 33 00");
	PRINT #f%, mkhex$("00 00 44 00");
	PRINT #f%, mkhex$("00 00 55 00");
	PRINT #f%, mkhex$("00 00 66 00");
	PRINT #f%, mkhex$("00 00 77 00");
	PRINT #f%, mkhex$("00 00 88 00");
	PRINT #f%, mkhex$("00 00 99 00");
	PRINT #f%, mkhex$("00 00 AA 00");
	PRINT #f%, mkhex$("00 00 BB 00");
	PRINT #f%, mkhex$("00 00 CC 00");
	PRINT #f%, mkhex$("00 00 DD 00");
	PRINT #f%, mkhex$("00 00 EE 00");
	PRINT #f%, mkhex$("00 00 FF 00");
    a$ = SPACE$(sw% \ 2)
    FOR y% = y2% TO y1% STEP -1
        x% = x1%
		WHILE x% < x2%
            a? = POINT(x%, y%)
            SHIFT LEFT a?, 4
            INCR x%
            a? = a? OR POINT(x%, y%)
            INCR x%
        	MID$(a$, (x% - x1%) \ 2, 1) = MKBYT$(a?)
    	WEND
        PRINT #f%, a$;
	NEXT y%
	CLOSE f%
	BEEP
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

FUNCTION serialno$ (i%, j%)
	'Produces a serial number string of length j% corresponding to i%
    a$ = STR$(i%)
    a$ = LTRIM$(a$)
    a$ = STRING$(j%, "0") + a$
    serialno$ = RIGHT$(a$, j%)
END FUNCTION

SUB whoop
	'Makes a "whooping" sound
	FOR i% = 1 TO 250
		SOUND 37 * SQR(i%), .04
	NEXT i%
END SUB

SUB morse (a$)
	'Sends the string a$ in Morse code
    DIM m$(255)
    f% = 750            'Frequency in Hz
    s% = 25             'Speed in words/minute
    m$(44) = "__..__"
    m$(46) = "._._._"
    m$(47) = "_.._."
    m$(48) = "_____"
    m$(49) = ".____"
    m$(50) = "..___"
    m$(51) = "...__"
    m$(52) = "...._"
    m$(53) = "....."
    m$(54) = "_...."
    m$(55) = "__..."
    m$(56) = "___.."
    m$(57) = "____."
    m$(63) = "..__.."
    m$(65) = "._"
    m$(66) = "_..."
    m$(67) = "_._."
    m$(68) = "_.."
    m$(69) = "."
    m$(70) = ".._."
    m$(71) = "__."
    m$(72) = "...."
    m$(73) = ".."
    m$(74) = ".___"
    m$(75) = "_._"
    m$(76) = "._.."
    m$(77) = "__"
    m$(78) = "_."
    m$(79) = "___"
    m$(80) = ".__."
    m$(81) = "__._"
    m$(82) = "._."
    m$(83) = "..."
    m$(84) = "_"
    m$(85) = ".._"
    m$(86) = "..._"
    m$(87) = ".__"
    m$(88) = "_.._"
    m$(89) = "_.__"
    m$(90) = "__.."
    FOR i% = 1 TO LEN(a$)
	b$ = MID$(a$, i%, 1)
	a% = ASC(UCASE$(b$))
	m$ = m$(a%)
	IF m$ = "" THEN
		SOUND 32767, 153 / s%
	ELSE
		SOUND 32767, 65 / s%
		FOR j% = 1 TO LEN(m$)
			SOUND 32767, 22 / s%
			IF MID$(m$, j%, 1) = "." THEN
			SOUND f%, 22 / s%
				ELSE
					SOUND f%, 65 / s%
		    END IF
		NEXT j%
	END IF
    NEXT i%
END SUB