'Program FIG10.BAS calculates quasiperiodic regions in a 3-D stereogram
'by J. C. Sprott <http://sprott.physics.wisc.edu/>
'Compile with the PowerBASIC Console Compiler <http://www.powerbasic.com/>
DECLARE FUNCTION gauss AS EXT
DEFEXT a-z      'do all calculations in extended (80-bit) precision

FUNCTION PBMAIN()
n&=3            'number of nonlinear odes
h=0.05          'stepsize of integration
tmax=1e3        'time to spend at each initial condition
pi=4#*ATN(1##)
s=pi            'scale factor on graph
CONSOLE NAME "Fig 10"
CONSOLE SET LOC 808,0
GRAPHIC WINDOW "Fig 10", 0, 0, 800, 400 TO hWin1&
GRAPHIC ATTACH hWin1&, 0
GRAPHIC WIDTH 1
GRAPHIC COLOR RGB(0,0,0), RGB(255,255,255)
GRAPHIC CLEAR
GRAPHIC ATTACH hWin1&, 0
GRAPHIC BOX(0,0)-(800,400)
GRAPHIC LINE(400,0)-(400,800)
GRAPHIC BITMAP LOAD "fig10.bmp", 800, 400 TO hBmp???
GRAPHIC COPY hBmp???, hWin1&
DIM xo(n&), x(n&), xb(n&), dxdt(n&),zb&(800,400)
RANDOMIZE TIMER
xb(2)=pi/2
eps=.01
WHILE INKEY$<>CHR$(27)  'try different random initial conditions
	FOR i&=1 TO n&      'explore the neighborhood of the last solution
		xo(i&)=xb(i&)+eps*gauss
		WHILE xo(i&)>pi: xo(i&)=xo(i&)-2##*pi: WEND
		WHILE xo(i&)<-pi: xo(i&)=xo(i&)+2##*pi: WEND
		x(i&)=xo(i&)
	NEXT i&
	eps=1.01##*eps
	IF eps>2## THEN eps=2##
	xp&=200*xo(1)
	zp&=200*xo(3)
	xp1&=200+(xp&+.05*zp&)/s
	xp2&=600+(xp&-.05*zp&)/s
	yp&=200-200*xo(2)/s
	IF xp1&<=1 THEN ITERATE
	IF xp2&>=799 THEN ITERATE
	IF yp&<=1 THEN ITERATE
	IF yp&>=399 THEN ITERATE
	zb1&=zb&(xp1&,yp&)           'this is the z-buffer
	zb2&=zb&(xp2&,yp&)
	IF zp&<zb1& AND zp&<zb2& AND zb1&<>0 THEN ITERATE
	zb&(xp1&,yp&)=zp&
	zb&(xp2&,yp&)=zp&
	FOR t=0 TO tmax/h
		CALL rk4(x(),dxdt(),n&,h,x())
	NEXT t
	d2=0: FOR i&=1 TO n&: d2=d2+x(i&)^2: NEXT i&
	x2=0: FOR i&=1 TO n&: x2=MAX(x2,x(i&)^2): NEXT i&
	IF d2>0.16*tmax^2 THEN      'orbit escaped
		IF x(1)^2>x(2)^2 AND x(1)^2>x(3)^2 THEN c&=RGB(255,0,0): IF x(1)<0 THEN c&=RGB(127,127,0)
		IF x(2)^2>x(3)^2 AND x(2)^2>x(1)^2 THEN c&=RGB(0,255,0): IF x(2)<0 THEN c&=RGB(0,255,255)
		IF x(3)^2>x(1)^2 AND x(3)^2>x(2)^2 THEN c&=RGB(0,0,255): IF x(3)<0 THEN c&=RGB(255,0,255)
		PRINT CLNG(d2),CLNG(d2-x2),ROUND(eps,3)
		IF zp&>zb1& THEN GRAPHIC SET PIXEL(xp1&,yp&),c&
		IF zp&>zb2& THEN GRAPHIC SET PIXEL(xp2&,yp&),c&
		FOR i&=1 TO n&: xb(i&)=xo(i&): NEXT i&
		eps=.01##
	END IF
WEND
GRAPHIC SAVE"fig10.bmp"
BEEP
WAITKEY$
END FUNCTION

SUB DERIVS (x(), dxdt(), n&)
	'Returns the time derivatives dxdt(i&) of x(i&)
	FOR i&=1 TO n&
		i1&=1+i& MOD n&
		dxdt(i&)=SIN(x(i1&))
	NEXT i&
END SUB

SUB RK4 (X(), DXDT(), N&, H, XNEW())
	'Fourth-order Runge-Kutta integrator
	DIM XT(N&), DXT(N&), DXM(N&)
	HH = H * .5##
	H6 = H / 6##
	CALL DERIVS(X(), DXDT(), N&)
	FOR I& = 1 TO N&
		XT(I&) = X(I&) + HH * DXDT(I&)
	NEXT I&
	CALL DERIVS(XT(), DXT(), N&)
	FOR I& = 1 TO N&
		XT(I&) = X(I&) + HH * DXT(I&)
	NEXT I&
	CALL DERIVS(XT(), DXM(), N&)
	FOR I& = 1 TO N&
		XT(I&) = X(I&) + H * DXM(I&)
		DXM(I&) = DXT(I&) + DXM(I&)
	NEXT I&
	CALL DERIVS(XT(), DXT(), N&)
	FOR I& = 1 TO N&
		XNEW(I&) = X(I&) + H6 * (DXDT(I&) + DXT(I&) + 2## * DXM(I&))
	NEXT I&
	ERASE DXM, DXT, XT
END SUB

FUNCTION gauss
	'Returns a normally (Gaussian) distributed random deviate
	'with mean of zero and standard deviation of 1.0
	DO
		V1 = 2 * RND - 1
		V2 = 2 * RND - 1
		R = V1 * V1 + V2 * V2
	LOOP WHILE R >= 1 OR R = 0
	GAUSS = V2 * SQR(-2 * LOG(R) / R)
END FUNCTION