'Program FIG13.BAS calculates probability distribution function and variance versus time
'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&=8192      'time to spend at each initial condition
s=100##         'scale factor on graph
pi=4#*ATN(1##)
CLS
CONSOLE NAME "Fig 13"
CONSOLE SET LOC 808,0
GRAPHIC WINDOW "Fig 13", 0, 0, 800, 640 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
CONSOLE SET FOCUS
DIM x(n&), dxdt(n&), p&(900), vartot(100)
OPEN"fig15.dat" FOR INPUT AS #2
    FOR i&=0 TO LOG2(tmax&): INPUT#2,ntot&,t&,vartot(i&): NEXT i&
CLOSE
RANDOMIZE TIMER
OPEN"fig13.dat" FOR APPEND AS #1
WHILE INKEY$<>CHR$(27)  ' try different random initial conditions
    x(1)=0.1-0.2*RND: x(2)=0.1-0.2*RND: x(3)=0.1-0.2*RND
    INCR ntot&
    OPEN"fig15.dat" FOR OUTPUT AS #2
    FOR t&=1 TO tmax&+1
        FOR i&=1 TO 1/h
            CALL rk4(x(),dxdt(),n&,h,x())
            xp&=400+400*(x(2)-.707*x(1))/s
            yp&=320-320*(x(3)-.707*x(1))/s
            GRAPHIC SET PIXEL(xp&,yp&)
        NEXT i&
        IF t&=4000 THEN             'This is for figure 13
            FOR i&=1 TO n&
                PRINT#1,TRIM$(STR$(CLNG(x(i&))))
            NEXT i&
        END IF
        logt&=INT(LOG2(t&))
        IF LOG2(t&)=logt& THEN      'This is for figure 15
            FOR i&=1 TO n&
                vartot(logt&)=vartot(logt&)+x(i&)*x(i&)
            NEXT i&
            PRINT#2,ntot&,t&,vartot(logt&)
        END IF
    NEXT t&
    CLOSE #2
    INCR j&
    LOCATE 1,1
    PRINT j&;
    GRAPHIC REDRAW
    GRAPHIC CLEAR
WEND
CLOSE #1
CLS
n&=0: v=0: kur=0
OPEN"fig13.dat" FOR INPUT AS #1
    WHILE NOT EOF(1)
        INPUT#1,x
        i&=x+400
        IF i&<0 THEN i&=0
        IF i&>800 THEN i&=800
        INCR p&(i&)
        IF i&>0 AND i&<799 THEN pmax&=MAX(pmax&,p&(i&))
        v=v+x*x
        INCR n&
    WEND
CLOSE
sd=SQR(v/n&)
OPEN"fig13.dat" FOR INPUT AS #1
    WHILE NOT EOF(1)
        INPUT#1,x
        kur=kur+(x/sd)^4
    WEND
CLOSE
kur=kur/n&-3
PRINT"N =";n&
PRINT"sd =";sd
PRINT"kur =";kur
FOR i&=1 TO 799
    j&=639-639*p&(i&)/pmax&
    IF i&=1 THEN GRAPHIC SET PIXEL(i&,j&) ELSE GRAPHIC LINE-(i&,j&)
NEXT i&
GRAPHIC REDRAW
GRAPHIC BOX(0,0)-(800,640)
GRAPHIC LINE(400,0)-(400,640)
GRAPHIC COLOR RGB(255,0,0)
FOR i&=1 TO 799
    x=i&-400
    pdf=n&*EXP(-x*x/(2*sd*sd))/(sd*SQR(2*pi))
    j&=639-639*pdf/pmax&
    IF i&=1 THEN GRAPHIC SET PIXEL(i&,j&) ELSE GRAPHIC LINE-(i&,j&)
NEXT i&
GRAPHIC REDRAW
GRAPHIC SAVE"fig13.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&)=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