'Program POINCARE.BAS displays the Poincare section for x(1)=0
'$error all on
DEFEXT a-z
n&=3
h=.02
pi=4##*atn(1##)
SCALE=0.15##
SCREEN 12
CALL invertcolors
LINE(0,0)-(639,479),,b
DIM x(n&), dxdt(n&), xnew(n&)
x3=0.808219##*pi
x(1)=0: x(2)=pi/2##: x(3)=x3
WHILE INKEY$<>CHR$(27)
    CALL rk4(x(),dxdt(),n&,h,xnew())
    IF xnew(1)*x(1)<0 THEN
        xs=xnew(2)-(xnew(2)-x(2))*xnew(1)/(xnew(1)-x(1))
        ys=xnew(3)-(xnew(3)-x(3))*xnew(1)/(xnew(1)-x(1))
        xp=320+240*(xs-1.57##)/SCALE
        yp=240-240*(ys-2.54)/SCALE
        IF xp>0 AND xp<639 THEN pset(xp,yp)
    END IF
    FOR i&=1 TO n&
    	x(i&)=xnew(i&)
        while x(i&)>pi: x(i&)=x(i&)-2##*pi: wend
        while x(i&)<-pi: x(i&)=x(i&)+2##*pi: wend
    NEXT i&
    incr t&
    if t&>3e5*2^((x3-0.808219*pi)/0.0032) then
    	x3=x3+0.0032##
    	x(1)=0: x(2)=pi/2##: x(3)=x3
        t&=0
    end if
WEND
CALL scn2bmp("fig11.bmp")
END

SUB invertcolors ()
    'Sets the RGB palette to iverted gray scale in 16-color VGA screen modes
    FOR pal& = 0 TO 15
        PALETTE pal&, pal&
        c& = 63-INT(4.2 * pal&)
        OUT &h3C8, pal&
        OUT &h3C9, c&
        OUT &h3C9, c&
        OUT &h3C9, c&
    NEXT pal&
END SUB

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

SUB scn2bmp (FILENAME$)
    'Captures the screen to a file filename$ in bitmaped (BMP) format
    'Palette correct in EGA and VGA color modes 7 - 12 only
    WINDOW SCREEN (0, 0)-(1, 1)             'Get width and height of screen
    w& = 1 + PMAP(1, 0)
    h& = 1 + PMAP(1, 1)
    WINDOW
    x1& = 0: y1& = 0: x2& = w&-1: y2& = h&-1        'Corners of image
    x2& = x1& + 16 * CEIL((x2& - x1&) \ 16 + .01) - 1
    SOUND 1000, .3
    SOUND 750, .3
    sw& = 1 + x2& - x1&                'Width of image in pixels
    sh& = 1 + y2& - y1&                'Height of image in pixels
    sf& = sw& * sh& \ 2                'Bytes in image
    f& = FREEFILE
    OPEN FILENAME$ FOR OUTPUT AS f&
    PRINT #f&, MKWRD$(&H4D42);          'ASCII "BM"
    PRINT #f&, MKL$(sf& + 118);         'Size of file in bytes
    PRINT #f&, MKDWD$(0);               'Must be zero
    PRINT #f&, MKWRD$(&h76)+MKWRD$(0)+MKWRD$(&H28)+MKWRD$(0);
    PRINT #f&, MKL$(sw&);               'Screen width
    PRINT #f&, MKL$(sh&);               'Screen height
    PRINT #f&, MKWRD$(1)+MKWRD$(4)+MKDWD$(0);
    PRINT #f&, MKL$(sf&);               'Size of image in bytes
    PRINT #f&, MKDWD$(0)+MKDWD$(0)+MKDWD$(0)+MKDWD$(0);
    FOR i& = 0 TO 15                    'Palette RGB values (0 - 15)
        OUT &h3C7,i&
        red?=4.047619*INP(&h3c9)
        grn?=4.047619*INP(&h3c9)
        blu?=4.047619*INP(&h3c9)
        PRINT #f&, MKBYT$(blu?);MKBYT$(grn?);MKBYT$(red?);MKBYT$(0);
    NEXT i&
    a$ = SPACE$(sw& \ 2)
    FOR y& = y2& TO y1& STEP -1
        x& = x1&
        WHILE x& < x2&
            a? = POINT(x&, y&)
            SHIFT LEFT a?, 4
            INCR x&
            a? = a? OR POINT(x&, y&)
            INCR x&
            MID$(a$, (x& - x1&) \ 2, 1) = MKBYT$(a?)
        WEND
        PRINT #f&, a$;
    NEXT y&
    CLOSE f&
    SOUND 800, .3
    SOUND 600, .3
END SUB