'Program SURFACE.BAS (c) 2007 by J. C. Sprott <sprott@physics.wisc.edu>
'Produces 24 frames of a surface described by z=z(x,y)
'Compile this program with PowerBASIC for DOS

defext a-z      	'Use extended (80-bit) precision throughout
dim p%(15)
surface$="Torus"	'Shape to animate (Saddle, Cube, Cyl, Parab, Hyperb, Torus)
nd&=10000      		'Number of initial points to discard
nf%=24				'Number of frames
thn%=0				'Starting angle (0 to nf% -1)
screen 7
call colorscale		'Reset color palette for red on cyan
c$=command$
more:
xmin=1e37: ymin=xmin: zmin=xmin
xmax=-xmin: ymax=-ymin: zmax=-zmin
n&=0
randomize 0			'Use the same random numbers for each frame
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
        cls
    end if
    select case lcase$(surface$)
    	case "saddle": x=1-2*rnd: y=1-2*rnd: z=x*x-y*y
        case "cube"
        	select case rnd
            	case <1/3: x=sgn(rnd-.5): y=1-2*rnd: z=1-2*rnd
                case >2/3: x=1-2*rnd: y=sgn(rnd-.5): z=1-2*rnd
                case else: x=1-2*rnd: y=1-2*rnd: z=sgn(rnd-.5)
            end select
        case "cyl": x=1-2*rnd: theta=6.2832*rnd: y=cos(theta): z=sin(theta)
        case "parab": theta=6.2832*rnd: y=rnd*sin(theta): z=rnd*cos(theta): x=y*y+z*z
    	case "hyperb": x=1-2*rnd: theta=6.2832*rnd: y=sqr(1+x*x)*cos(theta): z=sqr(1+x*x)*sin(theta)
    	case "torus": theta=6.2832*rnd: phi=6.2832*rnd: x=cos(theta)*(1+.5*cos(phi)): y=sin(theta)*(1+.5*cos(phi)): z=.5*sin(phi)
	end select
	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(surface$+serialno$(thn%,2)+".bmp"): exit loop
loop
incr thn%
print thn%: call morse(str$(thn%))
if thn%<24 then goto more
print"Done";: call whoop: call whoop
thn%=0
end

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

SUB colorscale ()
    'Sets the RGB palette to red on cyan 16-color VGA screen mode
    palette 0,0
    out &h3C8, 0
    out &h3C9, 0
    out &h3C9, 30
    out &h3C9, 30
    palette 1,1
    out &h3C8, 1
    out &h3C9, 0
    out &h3C9, 17
    out &h3C9, 17
    FOR pal% = 2 TO 15
        PALETTE pal%, pal%
        c% = INT(4.2 * pal%)
        OUT &h3C8, pal%
        OUT &h3C9, c%
        OUT &h3C9, 0 'c%
        OUT &h3C9, 0 'c%
    NEXT pal%
END SUB