'Program FIG16.BAS calculates Hurst exponent, autocorrelation function,
'power spectral density, and return map
'by J. C. Sprott <http://sprott.physics.wisc.edu/>
'Compile with the PowerBASIC Console Compiler <http://www.powerbasic.com/>
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&=1e9       'maximum time to follow trajectory
kmax&=2e3       'maximum delay for autocorrelation function
imax&=20        'number of frequency points (powers of 2)
pi=4##*ATN(1##)
CLS
CONSOLE NAME "Fig 16"
CONSOLE SET LOC 648,0
GRAPHIC WINDOW "Fig 16", 0, 0, 640, 480 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, REDRAW
GRAPHIC BOX(0,0)-(640,480)
ticks&=LOG10(tmax&)
FOR i&=1 TO ticks&-1
	xp&=640*i&/ticks&
	GRAPHIC LINE(xp&,0)-(xp&,20)
	GRAPHIC LINE(xp&,460)-(xp&,480)
	GRAPHIC LINE(0,480-xp&)-(20,480-xp&)
	GRAPHIC LINE(620,480-xp&)-(640,480-xp&)
NEXT i&
CONSOLE SET FOCUS
DIM x(n&), dxdt(n&), y(640), dx(kmax&), c(kmax&), sm(imax&), yp(imax&)
x(1)=0.2: x(2)=0: x(3)=0
FOR t&=1 TO tmax&
	FOR i&=1 TO 1/h
		CALL rk4(x(),dxdt(),n&,h,x())
	NEXT i&
	t0&=t& MOD (kmax&+1)
	dx(t0&)=dxdt(1)
	IF t&>kmax& THEN
		FOR k&=0 TO kmax&
			tt&=(t0&-k&+kmax&+1) MOD (kmax&+1)
			c(k&)=c(k&)+dx(t0&)*dx(tt&)
		NEXT k&
	END IF
	r2=0: FOR i&=1 TO n&: r2=r2+x(i&)*x(i&): NEXT i&
	r=MAX(r,SQR(r2))
	xp&=640*LOG10(t&)/ticks&
	yp&=480-640*LOG10(r)/ticks&
	IF r>rbest THEN
		rbest=r
		IF t&=1 THEN GRAPHIC SET PIXEL(xp&,yp&) ELSE GRAPHIC LINE-(xp&,yp&)
		IF xp&>xpold& THEN
			GRAPHIC REDRAW
			LOCATE 1,1
			PRINT t&
			xpold&=xp&
		END IF
		FOR i&=xp& TO 640: y(i&)=yp&: NEXT i&
		IF INKEY$=CHR$(27) THEN EXIT FOR
	END IF
NEXT t&
GRAPHIC LINE-(640,yp&)
CALL linreg(y(),a,b)
GRAPHIC LINE(0,b)-(640,a*640+b)
GRAPHIC SET POS(25,25)
GRAPHIC PRINT"R ="+STR$(10^((480-b)*(ticks&/640)))+"t^"+STR$(-a)
GRAPHIC REDRAW
GRAPHIC SAVE"fig16a.bmp"
BEEP
SLEEP 5000
GRAPHIC CLEAR
GRAPHIC BOX(0,0)-(640,480)
GRAPHIC LINE(0,240)-(640,240)
FOR i&=1 TO 9
	xp&=640*i&/10
	yp&=480*i&/10
	GRAPHIC LINE(xp&,0)-(xp&,20)
	GRAPHIC LINE(xp&,460)-(xp&,480)
	GRAPHIC LINE(0,yp&)-(20,yp&)
	GRAPHIC LINE(620,yp&)-(640,yp&)
NEXT i&
GRAPHIC SET PIXEL(0,0)
FOR k&=1 TO 100
	xp&=640*k&/100
	yp&=240-240*c(k&)/c(0)
	GRAPHIC LINE-(xp&,yp&)
NEXT tau&
GRAPHIC REDRAW
GRAPHIC SAVE"fig16b.bmp"
SLEEP 5000
GRAPHIC CLEAR
GRAPHIC BOX(0,0)-(640,480)
FOR i&=1 TO imax&-1
	xp&=640*i&/imax&
	yp&=640*i&/imax&
	GRAPHIC LINE(xp&,0)-(xp&,20)
	GRAPHIC LINE(xp&,460)-(xp&,480)
	GRAPHIC LINE(0,yp&)-(20,yp&)
	GRAPHIC LINE(620,yp&)-(640,yp&)
NEXT i&
FOR i&=0 TO imax&
	m&=2^i&
	sm(i&)=1+c(kmax&)/c(0)
	FOR k&=1 TO kmax&-1
		sm(i&)=sm(i&)+2*c(k&)*COS(0.5*pi*m&*k&/kmax&)/c(0)
	NEXT k&
	xp&=640*i&/imax&
	yp&=640*LOG2(sm(0)/sm(i&))/imax&
	yp(i&)=yp&
	IF i&=0 THEN GRAPHIC SET PIXEL(xp&,yp&) ELSE GRAPHIC LINE-(xp&,yp&)
NEXT i&
CALL linreg(yp(),a,b)
GRAPHIC LINE(0,b)-(640,a*imax&+b)
GRAPHIC SET POS(25,440)
GRAPHIC PRINT"a ="+STR$(a*imax&/640)
GRAPHIC REDRAW
GRAPHIC SAVE"fig16c.bmp"
GRAPHIC CLEAR
GRAPHIC BOX(0,0)-(640,480)
FOR i&=1 TO 7
	xp&=640*i&/8
	yp&=640*i&/8
	GRAPHIC LINE(xp&,0)-(xp&,20)
	GRAPHIC LINE(xp&,460)-(xp&,480)
	GRAPHIC LINE(0,yp&)-(20,yp&)
	GRAPHIC LINE(620,yp&)-(640,yp&)
NEXT i&
GRAPHIC LINE(0,240)-(640,240)
GRAPHIC LINE(320,0)-(320,480)
FOR t&=1 TO 5e4
	xp&=320+240*dxdt(1)
	FOR i&=1 TO 1/h
		CALL rk4(x(),dxdt(),n&,h,x())
	NEXT i&
	yp&=240-240*dxdt(1)
	GRAPHIC SET PIXEL(xp&,yp&)
NEXT i&
GRAPHIC REDRAW
GRAPHIC SAVE"fig16d.bmp"
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

SUB linreg (y(), a, b)
	'Fits a function y(x) to a linear form y = ax + b by least squares
	'x is assumed to be in equal (unit) intervals
	xmin& = LBOUND(y()): xmax& = UBOUND(y())
	n& = xmax& - xmin& + 1
	j& = 0: k = 0: l& = 0: r2 = 0
	FOR x& = xmin& TO xmax&
		y = y(x&)
		j& = j& + x&
		k = k + y
		l& = l& + x& * x&
		r2 = r2 + x& * y
	NEXT x&
	a = (n& * r2 - k * j&) / (n& * l& - j& * j&)
	b = (k - a * j&) / n&
END SUB