'Program WOLFMAP.BAS calculates the spectrum of LEs for the Henon map
'Ported from Wolf's Fortran code <http://www.cooper.edu/~wolf/chaos/chaos.htm>
'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()
CLS
n&=2            'number of variables in nonlinear map
nn&=n&*(n&+1)   'total number of variables (nonlinear + linear)
DIM x(nn&),xnew(nn&),v(nn&),ltot(n&),znorm(n&),gsc(n&)
irate&=10       'integration steps per reorthonormalization
io&=1e5         'number of iterations between printouts
FOR i&=1 TO n&  'initial conditions for nonlinear maps
	v(i&)=0     'must be within the basin of attraction
NEXT i&
t=0##
FOR i&=n&+1 TO nn&  'initial conditions for linearized maps
	v(i&)=0##   'Don't mess with these; they are problem independent!
NEXT i&
FOR i&=1 TO n&
	v((n&+1)*i&)=1##
	ltot(i&)=0##
NEXT i&
DO
	FOR j&=1 TO irate&
		FOR i&=1 TO nn&
			x(i&)=v(i&)
		NEXT i&
		CALL DERIVS(x(), xnew(), n&)
		FOR i&=1 TO nn&
			v(i&)=xnew(i&)
		NEXT i&
		INCR t
	NEXT j&
	'construct new orthonormal basis by gram-schmidt:
	znorm(1)=0##    'normalize first vector
	FOR j&=1 TO n&
		znorm(1)=znorm(1)+v(n&*j&+1)^2
	NEXT j&
	znorm(1)=SQR(znorm(1))
	FOR j&=1 TO n&
		v(n&*j&+1)=v(n&*j&+1)/znorm(1)
	NEXT j&
	'generate new orthonormal set:
	FOR j&=2 TO n&  'make j-1 gsr coefficients
		FOR k&=1 TO j&-1
			gsc(k&)=0##
			FOR l&=1 TO n&
				gsc(k&)=gsc(k&)+v(n&*l&+j&)*v(n&*l&+k&)
			NEXT l&
		NEXT k&
		FOR k&=1 TO n&  'construct a new vector
			FOR l&=1 TO j&-1
				v(n&*k&+j&)=v(n&*k&+j&)-gsc(l&)*v(n&*k&+l&)
			NEXT l&
		NEXT k&
		znorm(j&)=0##   'calculate the vector's norm
		FOR k&=1 TO n&
			znorm(j&)=znorm(j&)+v(n&*k&+j&)^2
		NEXT k&
		znorm(j&)=SQR(znorm(j&))
		FOR k&=1 TO n&  'normalize the new vector
			v(n&*k&+j&)=v(n&*k&+j&)/znorm(j&)
		NEXT k&
	NEXT j&
	FOR k&=1 TO n&  'update running vector magnitudes
		IF znorm(k&)>0 THEN ltot(k&)=ltot(k&)+LOG(znorm(k&))
	NEXT k&
	INCR m&
	IF (m& MOD io&)=0 THEN  'normalize exponent and print every io& iterations
		PRINT "Time =";CQUD(t);TAB(22);"LEs =";
		lsum=0##: kmax&=0
		FOR k&=1 TO n&
			le=ltot(k&)/t
			lsum=lsum+le
			IF lsum>0 THEN lsum0=lsum: kmax&=k&
			PRINT USING$("##.########  ", le);
		NEXT k&
		IF ltot(1)>0 THEN dky=kmax&-t*lsum0/ltot(kmax&+1) ELSE dky=0
		PRINT USING$("  Dky =##.########", dky) 'Kaplan-Yorke dimension
		IF INKEY$=CHR$(27) THEN EXIT LOOP   'exit when <Esc> key is pressed
	END IF
LOOP
END FUNCTION

SUB DERIVS(x(), xnew(), n&)
	a=1.4##
	b=0.3##
'   Nonlinear Henon map equations:
	xnew(1)=1##-a*x(1)*x(1)+b*x(n&)
	xnew(2)=x(1)
'   Linearized Henon map equations:
	xnew(3)=-2##*a*x(1)*x(3)+b*x(5)
	xnew(4)=-2##*a*x(1)*x(4)+b*x(6)
	xnew(5)=x(3)
	xnew(6)=x(4)
END SUB