'Program 3DLYAP.BAS (c) 2000 by J. C. Sprott
'Calculates Lyapunov exponent for 3d map
'Use PowerBasic Console Compiler
#DEBUG ERROR ON
#DEBUG DISPLAY ON
DEFEXT a-z

FUNCTION PBMAIN()
CLS
dr=1e-8                     'Separation of orbits
x=.1: y=.1: z=.1            'start on attractor
xe=x+dr: ye=y: ze=z         'Nearby initial condition
dr2=dr*dr
DO
    INCR n&&
    xnew=1.2##+3.5##*x-0.1##*z-0.4##*x*x+0.1##*z*z
    ynew=-0.5*z+0.9##*y*z
    z=1##*-0.1##*y-1.6*z*z
    x=xnew
    y=ynew
    xenew=1.2##+3.5##*xe-0.1##*ze-0.4##*xe*xe+0.1##*ze*ze
    yenew=-0.5*ze+0.9##*ye*ze
    ze=1##*-0.1##*ye-1.6*ze*ze
    xe=xenew
    ye=yenew
    dx=x-xe: dy=y-ye: dz=z-ze
    rs=SQR((dx*dx+dy*dy+dz*dz)/dr2)
    xe=x+dx/rs
    ye=y+dy/rs
    ze=z+dz/rs
    ltot=ltot+LOG(rs)
    le=ltot/n&&
    l1tot=l1tot+LOG(ABS(3.5##-0.8##*x))
    l1=l1tot/n&&
    IF (n&& MOD 1000000)=0 THEN
        PRINT FORMAT$(n&&,"##.#####^^^^");
        PRINT FORMAT$(le,"    #.#######");
        PRINT " ";CHR$(241);" ";
        PRINT FORMAT$(1/SQR(n&&),"#.#######");
        PRINT FORMAT$(l1,"    ##.#######");
        PRINT " ";CHR$(241);" ";
        PRINT FORMAT$(1/SQR(n&&),"#.#######")
        IF INKEY$=CHR$(27) THEN EXIT LOOP
    END IF
LOOP
END FUNCTION