'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% 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