'Program WOLFFLOW.BAS calculates the spectrum of LEs for the Lorenz attractor
'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&=3            'number of nonlinear odes
nn&=n&*(n&+1)   'total number of odes (nonlinear + linear)
DIM x(nn&),dxdt(nn&),v(nn&),A(nn&),B(nn&),C(nn&)
DIM D(nn&),ltot(n&),znorm(n&),gsc(n&)
irate&=20       'integration steps per reorthonormalization
h=0.01          'stepsize of integration
io&=5000        'number of iterations between printouts
v(1)=1.0        'initial conditions for nonlinear ODEs
v(2)=1.0        'must be within the basin of attraction
v(3)=1.0
t=0##
FOR i&=n&+1 TO nn&
    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
    'This is a built-in 4th order Runge-Kutta integrator:
    FOR j&=1 TO irate&
        FOR i&=1 TO nn&
            x(i&)=v(i&)
        NEXT i&
        CALL DERIVS(x(), dxdt(), n&)
        FOR i&=1 TO nn&
            A(i&)=dxdt(i&)
        NEXT i&
        FOR i&=1 TO nn&
            x(i&)=v(i&)+(h*A(i&))/2##
        NEXT i&
        CALL DERIVS(x(), dxdt(), n&)
        FOR i&=1 TO nn&
            B(i&)=dxdt(i&)
        NEXT i&
        FOR i&=1 TO nn&
            x(i&)=v(i&)+(h*B(i&))/2##
        NEXT i&
        CALL DERIVS(x(), dxdt(), n&)
        FOR i&=1 TO nn&
            C(i&)=dxdt(i&)
        NEXT i&
        FOR i&=1 TO nn&
            x(i&)=v(i&)+h*C(i&)
        NEXT i&
        CALL DERIVS(x(), dxdt(), n&)
        FOR i&=1 TO nn&
            D(i&)=dxdt(i&)
        NEXT i&
        FOR i&=1 TO nn&
            v(i&)=v(i&)+h*(A(i&)+D(i&)+2##*(B(i&)+C(i&)))/6##
        NEXT i&
        t=t+h
    NEXT j&
    'construct new orthonormal basis by gram-schmidt
    znorm(1)=0##    'normalize the 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 a new orthonormal set
    FOR j&=2 TO n&
        FOR k&=1 TO j&-1    'make j-1 gsr coefficients
            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 the 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(17);"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(), dxdt(), n&)
    b=8##/3##
    sg=10##
    r=28##
'   Nonlinear Lorenz equations:
    dxdt(1)=sg*(x(2)-x(1))
    dxdt(2)=-x(1)*x(3)+r*x(1)-x(2)
    dxdt(3)=x(1)*x(2)-b*x(3)
'   Linearized Lorenz equations:
    FOR i&=0 TO 2
        dxdt(4+i&)=sg*(x(7+i&)-x(4+i&))
        dxdt(7+i&)=(r-x(3))*x(4+i&)-x(7+i&)-x(1)*x(10+i&)
        dxdt(10+i&)=x(2)*x(4+i&)+x(1)*x(7+i&)-b*x(10+i&)
    NEXT i&
END SUB