'Program FIG5.BAS calculates standard deviation and kurtosis
'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
GLOBAL b AS EXT

FUNCTION PBMAIN()
n&=3            'number of nonlinear odes
h=0.01          'stepsize of integration
tmax=5e5        'time to spend at each b value
CONSOLE NAME "Fig 5"
CONSOLE SET LOC 648,0
GRAPHIC WINDOW "Fig 5", 0, 0, 640, 960 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(40,0)-(640,440)
GRAPHIC BOX(40,480)-(640,920)
GRAPHIC STYLE 0
FOR i&=1 TO 10
	xp&=40+600*i&/11
	GRAPHIC LINE(xp&,0)-(xp&,20)
	GRAPHIC LINE(xp&,420)-(xp&,440)
	GRAPHIC LINE(xp&,480)-(xp&,500)
	GRAPHIC LINE(xp&,900)-(xp&,920)
NEXT i&
FOR i&=1 TO 9
	yp&=440*i&/10
	GRAPHIC LINE(40,yp&)-(60,yp&)
	GRAPHIC LINE(620,yp&)-(640,yp&)
NEXT i&
FOR i&=1 TO 9
	yp&=480+440*i&/10
	GRAPHIC LINE(40,yp&)-(60,yp&)
	GRAPHIC LINE(620,yp&)-(640,yp&)
NEXT i&
DIM x(n&), dxdt(n&)
FOR i&=1 TO n&: x(i&)=RND: NEXT i&
bmin=0.11##
bmax=0##
FOR ib&=1 TO 600
	b=bmin+(bmax-bmin)*ib&/600
	v=0: k=0
	FOR t=0 TO tmax/h
		CALL rk4(x(),dxdt(),n&,h,x())
		d2=0
		FOR i&=1 TO n&: d2=d2+x(i&)*x(i&): NEXT i&
		v=v+d2
		var=v/t
		k=k+d2*d2/var^2
	NEXT t
	sd=SQR(var)
	kur=k/t-3
	PRINT b,sd,kur
	yp&=440-440*sd/20
	IF ib&=1 THEN yp1&=yp&
	GRAPHIC LINE(39+ib&,yp1&)-(40+ib&,yp&)
	yp1&=yp&
	yp&=480-440*kur/2
	IF ib&=1 THEN yp2&=yp&
	GRAPHIC LINE(39+ib&,yp2&)-(40+ib&,yp&)
	yp2&=yp&
NEXT ib&
GRAPHIC SAVE"fig5.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&)=-b*x(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