'Program HENONGP.BAS calculates Grassberger-Procaccia dimension for the Henon map
'(c) 1997 by J. C. Sprott (http://sprott.physics.wisc.edu/)
'Compile with PowerBASIC 3.2 or later
cls: print"Calculating...  Please wait."
defext a-z
dim h(6149), b&&(1023), d&&(1023), bzoom&&(128), dzoom&&(128)
sf=2^(1/16)        'scale factor
on error resume next
open"henongp.dat" for input as #1
for n%=0 to 1023: input#1,b&&(n%): d&&(n%)=0: next n%
close
open"zoom.dat" for input as #1
for n%=0 to 128: input#1,bzoom&&(n%): dzoom&&(n%)=0: next n%
close
on error goto 0
randomize timer
x=.01*rnd
y=.01*rnd
i&=0
bmin%=32000
do
    xnew=1-1.4*x*x+.3*y
    y=x
    incr i&
    j%=i& mod 6150
    if i&>7000 then            'orbit is on the attractor and array is full
        for k%=0 to 5999    'this is the inner-most time-consuming loop
        jp%=(j%+k%) mod 6150
            dx=x-h(jp%)
            r2#=dx*dx
            jp%=(jp%+1) mod 6150
            dx=xnew-h(jp%)
            r2#=r2#+dx*dx
            r2#=r2#*r2#: r2#=r2#*r2#: r2#=r2#*r2#
            b%=1048-peeki(varptr(r2#)+6)\16        'fast integer log2
            if b%<1024 then incr d&&(b%)
            if b%>100 then
                if b%>116 then
                    incr dzoom&&(128)
                else
                    r2#=r2#*r2#: r2#=r2#*r2#: r2#=r2#*r2#
                    b%=422-peeki(varptr(r2#)+6)\16
                    incr dzoom&&(b%)
                end if
            end if
        next k%
        if (i& mod 3000)=1000 then        'plot and record results occasionally
            for n%=0 to 1023: b&&(n%)=b&&(n%)+d&&(n%): next n%
            for n%=0 to 128: bzoom&&(n%)=bzoom&&(n%)+dzoom&&(n%): next n%
            call d2vsr(b&&())
            call zoom(bzoom&&())
            open"henongp.dat" for output as #1
            for n%=0 to 1023: print#1,b&&(n%): d&&(n%)=0: next n%
            close
            open"zoom.dat" for output as #1
            for n%=0 to 128: print#1,bzoom&&(n%): dzoom&&(n%)=0: next n%
            close
            dt=timer-t0
            if dt<0 then dt=dt+86400: t0=t0-86400
        end if
        call dimscreen(600)
        q$=lcase$(inkey$)
        if q$=chr$(27) then exit loop
        if q$="p" then call scn2bmp("henongp.bmp")
    end if
    x=xnew
    h(j%)=x
loop
end

sub d2vsr(b&&())
    'Takes the binned data in array b&&() and plots correlation dimen vs. r
    'Bin widths are assumed to be sf
    shared sf
    cls
    yp%=479-240*1.21
    b&&=0
    nfirst%=ubound(b&&)
    while b&&(nfirst%)=0: decr nfirst%: wend    'find first nonempty bin
    call graphbox(2*log2(sf)*nfirst%/6.643856,10)
    line(0,yp%)-(639,yp%),,,&H4444
    for n%=nfirst% to 0 step -1
        c&&=b&&
        b&&=b&&+b&&(n%)
        if c&&=0 then iterate
        d2=log2(b&&/c&&)/log2(sf)
        dd2=1.442695*sqr(1/c&&-1/b&&)/log2(sf)
        xp%=639-639*n%/nfirst%
        yp1%=479-240*(d2-dd2)
        yp2%=479-240*(d2+dd2)
        line(xp%,yp1%)-(xp%,yp2%),13
        yp%=479-240*d2
        if n%>100 and n%<117 then colr%=14 else colr%=15
        if n%<nfirst% then line(xpold%,ypold%)-(xp%,yp%),colr%
        xpold%=xp%: ypold%=yp%
    next n%
    bmed&&=b&&/2
    n%=0: c&&=0             'calculate median bin
    while c&&<b&&/2: c&&=c&&+b&&(n%): incr n%: wend
    xp%=639-639*n%/nfirst%
    line(xp%,0)-(xp%,479),,,&H4444
end sub

sub zoom(bzoom&&())
    'Takes the binned data in array bzoom&&() and plots correlation dimen
    'Bin widths are assumed to be sf
    sf=2^(1/128)            'scale factor
    b&&=bzoom&&(128)
    nfirst%=127
    while b&&(nfirst%)=0: decr nfirst%: wend 'find first nonempty bin
    for n%=nfirst% to 0 step -1
        c&&=b&&
        b&&=b&&+bzoom&&(n%)
        if c&&=0 then iterate
        d2=log2(b&&/c&&)/log2(sf)
        xp%=639-n%-150
        yp%=479-240*8*d2+2280
        if n%<nfirst% then line(xpold%,ypold%)-(xp%,yp%)
        xpold%=xp%: ypold%=yp%
    next n%
end sub

SUB DIMSCREEN (t!)
    'Dims the screen after t! seconds of keyboard inactivity
    'Of course you have to call the function occasionally
    STATIC sdim%, tdim!
    IF sdim% THEN           'The screen is dim
        IF INSTAT = 0 THEN  'Get out quickly
            EXIT SUB
        ELSE                'A key was pressed
            sdim% = 0
            tdim! = 0
            PALETTE
        END IF
    ELSE                    'The screen is bright
        IF tdim! = 0 THEN   'Start the counter
            tdim! = TIMER + t!
            IF tdim! > 86400 + t! THEN tdim! = tdim! - 86400 + t!
        END IF
        IF TIMER > tdim! THEN   'Dim it
            sdim% = -1
            PALETTE 8, 56
            PALETTE 9, 8
            PALETTE 10, 16
            PALETTE 11, 24
            PALETTE 12, 32
            PALETTE 13, 40
            PALETTE 14, 48
            PALETTE 15, 56
        END IF
    END IF
END SUB

SUB graphbox(xticks%, yticks%)
    'Puts box on screen x and y ticks for graphing
    SCREEN 12
    LINE (0, 0) - (639, 479), , b
    FOR i% = 1 TO xticks% - 1
        LINE (640 * i% / xticks%, 1) - (640 * i% / xticks%, 10)
        LINE (640 * i% / xticks%, 469) - (640 * i% / xticks%, 478)
    NEXT i%
    FOR j% = 1 TO yticks% - 1
        LINE (1, 480 * j% / yticks%) - (10, 480 * j% / yticks%)
        LINE (629, 480 * j% / yticks%) - (638, 480 * j% / yticks%)
    NEXT j%
END SUB

FUNCTION mkhex$ (a$)
    'Converts a sequence of space-delimited hexadecimal bytes to a string
    'For example, mkhex$("4A 4B 4C") evaluates to "JKL"
    bytes% = (LEN(a$) + 1) / 3
    c$ = ""
    FOR i% = 1 TO bytes%
        b$ = "&H" + MID$(a$, 3 * i% - 2, 3)
        c$ = c$ + CHR$(VAL(b$))
    NEXT
    mkhex$ = c$
END FUNCTION

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
    x1% = 0: y1% = 0: x2% = 639: y2% = 479        'Corners of image
    x2% = x1% + 16 * ceil((x2% - x1%) \ 16 + .01) - 1
    beep
    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%, mkhex$("42 4D");         'ASCII "BM"
    PRINT #f%, MKL$(sf& + 118);         'Size of file in bytes
    PRINT #f%, mkhex$("00 00 00 00");   'Must be zero
    PRINT #f%, mkhex$("76 00 00 00 28 00 00 00");
    PRINT #f%, MKL$(sw%);               'Screen width
    PRINT #f%, MKL$(sh%);               'Screen height
    PRINT #f%, mkhex$("01 00 04 00 00 00 00 00");
    PRINT #f%, MKL$(sf&);               'Size of image in bytes
    PRINT #f%, mkhex$("00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00");
    PRINT #f%, mkhex$("00 00 00 00");    'Palette RGB values (0 - 15)
    PRINT #f%, mkhex$("AA 00 00 00");
    PRINT #f%, mkhex$("00 AA 00 00");
    PRINT #f%, mkhex$("AA AA 00 00");
    PRINT #f%, mkhex$("00 00 AA 00");
    PRINT #f%, mkhex$("AA 00 AA 00");
    PRINT #f%, mkhex$("00 55 AA 00");
    PRINT #f%, mkhex$("AA AA AA 00");
    PRINT #f%, mkhex$("FF FF FF 00");
    PRINT #f%, mkhex$("FF 55 55 00");
    PRINT #f%, mkhex$("55 FF 55 00");
    PRINT #f%, mkhex$("FF FF 55 00");
    PRINT #f%, mkhex$("55 55 FF 00");
    PRINT #f%, mkhex$("FF 55 FF 00");
    PRINT #f%, mkhex$("55 FF FF 00");
    PRINT #f%, mkhex$("FF FF FF 00");
    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%
    beep
END SUB