'Program ELEGFRAC.BAS (c) 2018 by J. C. Sprott 'Produces and examines images for the book ELEGANT FRACTALS 'Use PowerBasic Console Compiler '#DEBUG ERROR ON '#DEBUG DISPLAY ON DEFEXT a-z GLOBAL pi##,seed&,c&() FUNCTION PBMAIN() AS LONG ne&=6 'Number of variables np&=12 'Number of parameters w&=1200 'Width of plot h&=900 'Height of plot nd&=4 'Number of digits in printout ipp&=10 'Number of iterations per pixel ni&=1e5 'Number of iterations between printouts mp&=1e4 'Maximum period for periodicity checking bail&=1e3 'Bailout value for escape time plot dr=2^-26 'Separation of orbits for LE calculation tc$="" 'Test case iebest&=4*np& 'A very elegant equation pi=4##*ATN(1##) '3.1415926... DIM x(ne&),xe(ne&),a(np&),xmin(ne&),xmax(ne&),xd(9999,ne&),lnd&&(255) DIM c&(7) 'A simple color palette c&(0)=%RGB_WHITE c&(1)=%RGB_RED c&(2)=%RGB_GREEN c&(3)=%RGB_BLUE c&(4)=%RGB_YELLOW c&(5)=%RGB_CYAN c&(6)=%RGB_MAGENTA c&(7)=RGB(220,220,220) 'Shadow color DESKTOP GET CLIENT TO wc&, hc& CONSOLE GET SIZE TO x&, y& x&=MIN(w&-2,wc&-400) 'Keep console on the screen CONSOLE SET LOC x&,0 COLOR 0,15 CLS LOCAL Built AS IPOWERTIME 'Version number is the date compiled LET Built = CLASS "PowerTime" Built.FileTime = %PB_COMPILETIME ver$=Built.DateString+" Version" DO PRINT TAB(28)" ELEGANT FRACTALS" PRINT TAB(28)"(c)2018 J. C. Sprott" PRINT TAB(39-LEN(ver$)/2)ver$ PRINT PRINT" This program allows you to replicate all the fractal images in the book" PRINT" `Elegant Fractals: Automated Generation of Computer Art' by J. C. Sprott" PRINT" (World Scientific Press, 2018)." PRINT PRINT" It also allows you to produce countless additional examples using the" PRINT" automated search techniques described in the book." LOCATE 12,1 PRINT,," H: Help" PRINT,," M: Modify" PRINT,," S: Search" PRINT,," V: Verify" PRINT,," Esc: Exit program" PRINT PRINT,," Make selection: "; q$=UCASE$(WAITKEY$("hHmMsSvV"+$ESC)) IF q$="H" THEN helpme(q$,ver$) IF q$="M" THEN modify(w&,h&,nd&,ipp&,ni&,mp&,bail&) IF q$=$ESC THEN EXIT FUNCTION LOOP UNTIL q$<>"H" AND q$<>"M" DIM xold(mp&-1),xy(2*w&*h&-1) GRAPHIC WINDOW EXE.FULL$, 0, 0, w&, h& TO hWin??? GRAPHIC ATTACH hWin???, 0, REDRAW GRAPHIC SCALE PIXELS GRAPHIC COLOR %BLACK, %WHITE CONSOLE SET FOCUS CLS SELECT CASE q$ CASE "S" PRINT,,"Type of System" PRINT PRINT" A: 6-D quadratic map N: Generalized Lorenz system" PRINT" B: 6-D cubic map O: 2-D Iterated function system" PRINT" C: 3-D cubic map P: 3-D Iterated function system" PRINT" D: 6-D signum map Q: 2-D quadratic IFS" PRINT" E: 4-D sine map R: 3-D cubic IFS" PRINT" F: 4-D tangent map S: 2-D quintic complex map" PRINT" G: 4-D hyperbolic sine map T: Deterministic cellular automaton" PRINT" H: 4-D hyperbolic tangent map U: 1-D coupled map lattice" PRINT" I: 6-D integer map V: Random walk" PRINT" J: 3-D quadratic flow W: Diffusion limited aggregation" PRINT" K: 3-D cubic flow X: Stoschastic cellular automaton" PRINT" L: 4-D sine flow Y: Percolation" PRINT" M: 4-D hypertangent flow Z: Fractal tessellation" PRINT PRINT" =: A range of the above Esc: Exit program" PRINT PRINT,,"Make selection: "; t$=UCASE$(WAITKEY$("aAbBcCdDeEfFgGhHiIjJkKlLmMnNoOpPqQrRsStTuUvVwWxXyYzZ="+$ESC+$CR)) IF t$=$ESC THEN EXIT FUNCTION IF t$=$CR THEN tstart$="A" tend$="Z" ELSE tstart$=t$ tend$=t$ END IF IF t$="=" THEN PRINT PRINT,,"Starting case: "; tstart$=UCASE$(WAITKEY$) IF tstart$=$CR THEN tstart$="A" PRINT tstart$ PRINT,,"Ending case: "; tend$=UCASE$(WAITKEY$) IF tend$=$CR THEN tend$="Z" PRINT tend$ END IF CLS PRINT,," Type of Search" PRINT PRINT,," A: Alternate coloring" PRINT,," B: Basin of attraction" PRINT,," C: Escape-time contours" PRINT,," D: Escape-time (alt color)" PRINT,," E: Elegant equation" PRINT,," I: Search for icon" PRINT,," J: Icon with alternate color" PRINT,," N: Normal search" PRINT,," S: Select for symmetry" PRINT,,"Esc: Exit program" PRINT PRINT,,"Make selection: "; ic$=UCASE$(WAITKEY$("aAbBcCdDeEiIjJnNsS"+$ESC+$CR)) IF ic$=$ESC THEN EXIT FUNCTION IF ic$=$CR THEN ic$="N" CLS PRINT" Searching type ";tstart$; IF tend$>tstart$ THEN PRINT" to ";tend$; PRINT" for "; SELECT CASE ic$ CASE "A": PRINT"alternate color plot..." CASE "B": PRINT"basin of attraction..." CASE "C": PRINT"escape time contours..." CASE "D": PRINT"escape time contours with alternate color scheme..." CASE "E": PRINT"elegant equation..." CASE "I": PRINT"symmetric icon..." CASE "J": PRINT"icon with alternate coloring... CASE "N": PRINT"normal plot..." CASE "S": PRINT"symmetric image..." END SELECT CASE "V" a$=COMMAND$ 'Allow entry of code on command line IF tc$<>"" THEN a$=tc$ 'An optional hard coded fractal code DO 'Enter code to verify and check its length IF a$="" THEN LINE INPUT"Type the Code for the image you want to display: ",a$ IF LEN(a$)<13 THEN PRINT"String too short!": a$="": ITERATE IF LEN(a$)=14 AND INSTR("123456789&",RIGHT$(a$,1))=0 THEN PRINT"String too long!": a$="": ITERATE IF LEN(a$)>14 THEN PRINT"String too long!": a$="": ITERATE EXIT LOOP LOOP a$=UCASE$(a$) 'Code is not case sensitive t$=LEFT$(a$,1) 'Prefix b$=RIGHT$(a$,1) 'Suffix CLS PRINT,," Type of Display" PRINT IF b$="&" THEN PRINT,," B: Basin of attraction" PRINT,," C: Escape-time contours" PRINT,," D: Escape-time (alt color)" PRINT,,"Esc: Exit program" PRINT PRINT,,"Make selection: "; ic$=UCASE$(WAITKEY$("bBcCdD"+$ESC)) ELSE IF INSTR("STUVWXYZ",t$) THEN EXIT IF PRINT,," A: Alternate color plot" PRINT,," J: Icon with alternate coloring" PRINT,," N: Normal plot PRINT,,"Esc: Exit program PRINT PRINT,,"Make selection: "; ic$=UCASE$(WAITKEY$("aAjJnN"+$ESC)) END IF IF ic$=$ESC THEN EXIT FUNCTION CLS PRINT" Plotting ";a$;"..." ns&=VAL(b$) 'Number of icon segments END SELECT DO REDIM zb(w&,h&) 'Reset the z-buffer REDIM xy(2*w&*h&-1) 'and the screen pixels GRAPHIC CLEAR seed&=-24850##*TIMER 'Reseed the random number generator IF VAL(s$)>0 THEN INCR pc& IF s$="C" THEN INCR sa& IF s$="U" THEN INCR ub& s$="C" 'Default solution type SELECT CASE q$ CASE"S" DO 'Control the frequency of testing certain cases during a random search t$=CHR$(INT(ASC(tstart$)+(ASC(tend$)-ASC(tstart$)+1)*RAN)) SELECT CASE t$ CASE "B","H","K": IF ran<0.5## THEN EXIT LOOP CASE "C","P": IF ran<0.2## THEN EXIT LOOP CASE "D","E","O","Q": IF ran<0.1## THEN EXIT LOOP CASE "R": IF ran<0.05## THEN EXIT LOOP CASE "U": IF ran<0.01## THEN EXIT LOOP CASE "T","W","Y","Z": IF ran<0.002## THEN EXIT LOOP CASE "V","X": IF ran<0.001## THEN EXIT LOOP CASE ELSE: EXIT LOOP END SELECT LOOP a$=t$ FOR i&=1 TO np& r&=INT(26##*ran) 'Choose a random letter of the alphabet (0 to 25) IF t$="T" THEN IF i&=1 THEN r&=r& MOD 8 IF i&>11 THEN r&=12 END IF IF t$="U" THEN IF i&=1 THEN r1&=r& IF i&>11 AND (r1& MOD 13)>3 AND r1&<>25 THEN r&=12 END IF IF t$="V" OR t$="W" THEN IF i&=1 THEN IF ic$="N" THEN r&=0 IF ic$="I" THEN r&=1 IF ic$="A" THEN r&=2 IF ic$="J" THEN r&=3 END IF IF t$="V" AND i&>7 THEN r&=12 IF t$="W" AND i&=12 THEN r&=10*CINT(r&/12##)+2 END IF IF t$="X" OR t$="Y" THEN IF i&=1 THEN r&=12 END IF IF t$="Z" THEN IF INSTR("ABCDEFG",MID$(a$,2,1)) THEN IF i&>9 THEN r&=12 ELSEIF MID$(a$,2,1)="H" THEN IF i&>10 THEN r&=12 ELSE IF i&=12 THEN r&=12 END IF END IF a$=a$+CHR$(r&+65) 'Add it to the string identifier a(i&)=0.1##*(r&-12) 'And set the parameter in the range (-1.2 to 1.3) NEXT i& IF ic$="I" OR ic$="J" THEN 'It's an icon ns&=CEIL(9*ran) 'Number of icon segments a$=a$+CHR$(48+ns&) ELSE ns&=0 END IF IF ic$="B" THEN a$=a$+"&" 'It's a basin of attraction plot IF ic$="C" OR ic$="D" THEN a$=a$+"&" 'It's an escape time plot PRINT " Periodic:"+STR$(pc&)+" Chaotic:"+STR$(sa&)+" Unbounded:"+STR$(ub&); LOCATE,1 CASE"V" b$=RIGHT$(a$,1) FOR i&=1 TO np& r&=ASC(MID$(a$,i&+1,1))-65 a(i&)=0.1##*(r&-12) IF i&=1 AND (t$="V" OR t$="W") THEN IF r&=0 THEN ic$="N" IF r&=1 THEN ic$="I" IF r&=2 THEN ic$="A" IF r&=3 THEN ic$="J" END IF NEXT i& END SELECT dum???=rani???(a()) 'Choose reproducible seed for random number generator seed&=-(dum??? MOD 2^31) ie&=inelegance&(a$,np&) IF ic$="E" THEN IF ie&>iebest& THEN ITERATE 'Disregard inelegant equations FOR i&=1 TO ne& 'Set initial conditions x(i&)=0.1## xmin(i&)=100## xmax(i&)=-xmin(i&) NEXT i& sigsq=0##: skew=0## nmax&=10*w&*h& 'Every pixel can be visited ten times n1&=nmax&/1e3 'Iterates to discard to get on attractor n2&=nmax&/10 'Iterates to discard to get bounds FOR n&=0 TO nmax& gb&=0 IF ic$="C" OR ic$="D" THEN 'Make escape time contours IF q$="S" THEN 'See if it's worth plotting REDIM x(ne&) 'Set all initial conditions to zero x(1)=-1##+2##*ran 'Pick a random point in the plot x(2)=-1##+2##*ran FOR i&=1 TO bail& map(n&,x(),a(),ne&,t$,h,xy(),w&,h&) IF x(1)*x(1)+x(2)*x(2)>1e4 OR VAL(s$)>0 THEN IF i&n2& AND q$="S" THEN FOR i&=1 TO 2 IF x(i&)>xmax(i&) THEN s$="U" 'Solution is off scale IF x(i&)1e4 THEN s$="U" 'Solution is unbounded IF n&>2*w& AND (n& MOD mp&)=0 AND ic$<>"C" AND INSTR("TWXYZ",t$)=0 THEN FOR j&=mp&-1 TO 0 STEP -1 'Test for periodic cycle jj&=(ii&-j&+mp&) MOD mp& IF ROUND(x,16)=ROUND(xold(jj&),16) THEN s$=STR$(j&+1) 'j&+1 cycle found NEXT j& END IF ii&=(ii&+1) MOD mp& xold(ii&)=x IF (n& MOD w&)=w&-1 THEN GRAPHIC REDRAW IF INKEY$=$ESC THEN EXIT FUNCTION IF GRAPHIC$(INKEY$)=$ESC THEN EXIT FUNCTION END IF IF q$="V" AND ic$="B" THEN ic$="C" 'Plot contours even if there is no attractor IF s$<>"C" THEN IF q$="S" THEN EXIT, ITERATE IF q$="V" THEN PRINT" Period ";s$: EXIT, EXIT END IF x=x(1) IF n&=w& OR yp<=0 OR yp>=h& THEN ITERATE IF zp>zb(xp,yp) THEN IF ns& THEN makeicon(xp,yp,w&,h&,ns&) IF ic$="B" AND h<1## AND x(3)*zold>0 THEN EXIT IF 'Make a cross section at z=0 IF ic$="A" OR ic$="J" THEN GRAPHIC SET PIXEL(xp,yp) ELSE setpixel(xp,yp,zp,upp,vp,wp) shadow (xp,yp,zp,w&,h&,h) END IF zb(xp,yp)=zp 'Save z-value of pixel END IF NEXT n& IF ic$<>"C" AND ic$<>"D" THEN IF t$="T" THEN 'Alternate PF and CP calculation for CA pf=0## FOR n&=nmax&-w&*h& TO nmax& ij&=n& MOD (2*w&) IF xy(ij&)<>0## THEN INCR pf xp&=n& MOD w& yp&=INT(n&/w&) MOD h& zb(xp&,yp&)=xy(ij&) NEXT i& pf=pf/(w&*h&) ELSEIF t$="U" THEN 'Alternate PF and CP calculation for CML pf=0## FOR n&=nmax&-w&*h& TO nmax& ij&=n& MOD (2*w&) IF xy(ij&)>0## THEN INCR pf xp&=n& MOD w& yp&=INT(n&/w&) MOD h& zb(xp&,yp&)=MAX(0,xy(ij&)) NEXT i& pf=pf/(w&*h&) ELSE pf=pixfrac(w&,h&) 'Pixel fraction END IF IF bs=1## THEN pf=9##*pf 'Correct pf for nonzero boundary IF t$="X" OR t$="Y" OR t$="Z" THEN 'Alternate CP calculation cp=0 FOR i&=1 TO w&-2 FOR j&=1 TO h&-2 GRAPHIC GET PIXEL(i&,j&) TO c& GRAPHIC GET PIXEL(i&+1,j&) TO c1& GRAPHIC GET PIXEL(i&,j&-1) TO c2& GRAPHIC GET PIXEL(i&-1,j&) TO c3& GRAPHIC GET PIXEL(i&,j&+1) TO c4& IF c1&=c& AND c2&=c& AND c3&=c& AND c4&=c& THEN INCR cp NEXT j& NEXT i& cp=cp/(w&*h&-2*w&-2*h&) ELSE cp=cluster(zb(),w&,h&) 'Cluster probability END IF sk=skew/(nold&-n2&) 'Skewness IF t$="Y" OR t$="Z" THEN sk=0## 'Ignore skewness for these cases END IF IF ic$="A" OR ic$="J" OR t$="Z" THEN IF q$="V" THEN PRINT" Coloring..." FOR i&=0 TO w&-1 FOR j&=0 TO h&-1 GRAPHIC GET PIXEL(i&,j&) TO c& IF c&<>%WHITE THEN ITERATE c&=c&(CEIL(6*ran)) GRAPHIC PAINT BORDER (i&,j&),c&,%RGB_BLACK NEXT j& GRAPHIC REDRAW NEXT i& IF ran>0.5## THEN 'Remove background color GRAPHIC PAINT BORDER (0,0),%WHITE,%BLACK IF t$="Z" THEN GRAPHIC PAINT BORDER (0,0),c&(7),%BLACK 'Gray gasket for tessellation GRAPHIC REDRAW END IF END IF IF q$="S" AND gb&=0 THEN IF pf<0.03 AND t$<>"Z" THEN ITERATE 'Not likely to be good (too thin) IF t$="T" THEN IF pf>0.97 THEN ITERATE 'Cellular automaton is mostly ones IF pf>0.49 AND pf<0.51 THEN ITERATE 'CA is probably periodic END IF IF ic$="S" THEN IF ABS(sk)>0.01 THEN ITERATE 'Keep only symmetric cases 'IF (cp-1##/3##)^2+(pf-1##/3##)^2>(1##+8##*ran)/9## THEN ITERATE 'Prefer cp and pf near 1/3 END IF GRAPHIC REDRAW p$=" "+a$+" IE ="+STR$(ie&) IF ic$="C" OR ic$="D" THEN p$=p$+SPACE$(22) ELSE p$=p$+" PF = "+USING$("#.###",pf)+" CP = "+USING$("#.###",cp)+" SK ="+USING$("##.###",sk) END IF IF q$="V" THEN CLS PRINT p$ IF ic$="B" THEN 'Plot basin of attraction escapetime(x(),a(),ne&,t$,h,w&,h&,bail&,xmin(1),xmax(1),xmin(2),xmax(2),ic$,mp&) FOR i&=1 TO ne&: x(i&)=0.1##: NEXT i& END IF OPEN EXE.NAME$+".txt" FOR APPEND AS #1 PRINT #1, DATE$;" ";TIME$;" ";p$ 'Save results to a file CLOSE #1 IF q$="V" THEN GRAPHIC SAVE EXE.NAME$+".bmp" EXIT LOOP END IF iebest&=MIN(iebest&,ie&) BEEP GRAPHIC SAVE a$+".bmp" SLEEP 5000 'Allow time (5 seconds) to look at result LOOP PRINT" Press any key to show detailed results: "; q$=WAITKEY$ PRINT PRINT" Time x";SPACE$(nd&+4);"LE";SPACE$(nd&+3);"D2" REDIM zb(w&,h&) 'Reset the z-buffer GRAPHIC CLEAR n&&=0 ltot=0## FOR j&=1 TO ne& xe(j&)=x(j&)+dr/SQR(ne&) 'Perturb nearby initial condition NEXT j& totime=0## DO INCR n&& lle=locle(x(),xe(),a(),ne&,dr,t$,h,xy(),w&,h&,n&&) 'Local largest Lyapunov exponent IF ABS(x(1))>1e6 THEN PRINT" Unbounded! ";: WAITKEY$: EXIT LOOP ltot=ltot+lle*h 'Weight LE by time spent there totime=totime+h 'Total elapsed time le=ltot/totime 'Global largest Lyapunov exponent d2=corrdim(x(),xd(),lnd&&(),ne&,ntot&&) 'Correlation dimension xp=(w&-2)*(x(1) - xmin(1))/dx yp=(h&-2)*(xmax(2) - x(2))/dy zp=(x(3) - xmin(3))/(xmax(3)-xmin(3)) upp=(x(4) - xmin(4))/(xmax(4)-xmin(4)) '"up" is a reserved word in BASIC vp=(x(5) - xmin(5))/(xmax(5)-xmin(5)) wp=(x(6) - xmin(6))/(xmax(6)-xmin(6)) IF xp<=0 OR xp>=w& OR yp<=0 OR yp>=h& THEN ITERATE IF zp>zb(xp,yp) AND ic$<>"C" AND ic$<>"D" THEN IF ns& THEN makeicon(xp,yp,w&,h&,ns&) setpixel(xp,yp,zp,upp,vp,wp) zb(xp,yp)=zp 'Save z-value of pixel shadow (xp,yp,zp,w&,h&,h) END IF IF n&& MOD ni& THEN ITERATE q$=INKEY$ IF q$=$ESC THEN EXIT LOOP IF q$=" " THEN WAITKEY$ GRAPHIC REDRAW PRINT USING$("##."+STRING$(MAX(1,INT(LOG10(n&&/ni&))),"#")+"^^^^",totime); PRINT USING$("####."+STRING$(nd&,"#"),x(1),le,d2) LOOP END FUNCTION FUNCTION cluster(zb(),w&,h&) 'Calculates the probability that a point is part of a cluster FOR i&=1 TO w&-1 FOR j&=1 TO h&-1 z=zb(i&,j&) IF z=0## THEN ITERATE INCR n& IF zb(i&,j&-1) AND zb(i&-1,j&) AND zb(i&,j&+1) AND zb(i&+1,j&) THEN INCR c& NEXT j& NEXT i& FUNCTION=c&/n& END FUNCTION FUNCTION corrdim(x(),xd(),lnd&&(),ne&,ntot&&) 'Estimates the correlation dimension r%=INT(10000*ran) 'Pick a random point from the bank ep2=0## 'and calculate its Euclidean separation FOR j&=1 TO ne& 'from the current point dx=x(j&)-xd(r%,j&) ep2=ep2+dx*dx NEXT j& IF ep2<2^-127 THEN ep2=2^-127 'Prevent underflow if separation is too small ep2=ep2*2^-128 k%=255-PEEK(VARPTR(ep2)+8) 'Fast calculation of INT(-LOG2(ep2) INCR lnd&&(k%) 'Add separation to bin k% n&&=0 'Estimate correlation dimension k1%=0 FOR k%=255 TO 0 STEP -1 n&&=n&&+lnd&&(k%) 'This is the correlation sum IF n&&>10^3 AND k%<255 AND k1%=0 THEN n1&&=n&&: k1%=k% IF n&&>MAX(10^6,2*n1&&) AND k%1e4 THEN IF ic$="D" THEN dn=LOG(r2/1e4)/LOG(r2/r2old) 'interpolated correction to n& hue=3##*pi*LOG(n&-dn)/LOG(bail&) 'alternate coloring scheme hsi(hue,1,128,r&,g&,b&) ELSE pickcolor(n&,128,r&,g&,b&) END IF GRAPHIC SET PIXEL(i&,j&),RGB(r&,g&,b&) EXIT FOR END IF r2old=r2 NEXT n& NEXT j& GRAPHIC REDRAW NEXT i& FOR i&=1 TO ne&: x(i&)=0.1##: NEXT i& 'Restore initial conditions END SUB SUB helpme(q$,ver$) 'Opens a help window with additional instructions TXT.WINDOW("ELEGFRAC Help", 0, 0) TO hWin& TXT.PRINT TAB(28)" ELEGANT FRACTALS" TXT.PRINT TAB(28)"(c)2018 J. C. Sprott" TXT.PRINT TAB(39-LEN(ver$)/2)ver$ TXT.PRINT t$="This program has two modes: SEARCH and VERIFY. In the SEARCH mode, it will " t$=t$+"search for new examples of fractals, either of a specific type or for a " t$=t$+"range of types. The cases found are saved as .bmp graphic files, and are " t$=t$+"listed in the elegfrac.txt file. Only those cases that the program judges as " t$=t$+"elegant are saved. In the VERIFY mode, the program will allow you to " t$=t$+"type a specific 13 or 14 character code to reproduce a case that was " t$=t$+"previously found such as those shown in the ELEGANT FRACTALS book (World " t$=t$+"Scientific, 2018). It also includes a command MODIFY that allows you to " t$=t$+"change many of the parameters used in the program such as the size of the " t$=t$+"graphic image. This is especially useful if you want to make very high " t$=t$+"resolution plots for printing or other artistic purposes. Of course you " t$=t$+"will sacrifice calculation speed when you increase the resolution, and so " t$=t$+"you may want to use that sparingly. When searching for new fractals, " t$=t$+"you will probably want to run the program ovenight and then examine the few " t$=t$+"dozen cases found and judged to be elegant. Some systems calculate much " t$=t$+"faster than others. This program is provided as-is without any warranty or " t$=t$+"support. It is copyrighted, but you are free to use it for noncommercial " t$=t$+"purposes. If you want more detail about the program or the fractal images " t$=t$+"that it produces, you will need to purchase the ELEGANT FRACTALS book from " t$=t$+"the publisher. The current Windows executable version of the program along " t$=t$+"with the PowerBASIC source code is available from " t$=t$+"http://sprott.physics.wisc.edu/fractals/elegant/." DO 'Wrap the text to 80 columns with line breaks t$=TRIM$(t$) IF LEN(t$)<80 THEN TXT.PRINT t$ EXIT LOOP END IF c&=INSTR(80-LEN(t$),t$," ") b$=LEFT$(t$,c&-1) TXT.PRINT b$ t$=MID$(t$,LEN(b$)+1) LOOP q$="" WHILE q$="" IF TXT.INSTAT THEN q$=UCASE$(TXT.INKEY$) IF INSTAT THEN q$=UCASE$(INKEY$) SLEEP 100 WEND TXT.END CLS END SUB SUB hsi(h,s,i&,r&,g&,b&) 'Converts hue, saturation, and intensity (h,s,i&) into (r,g,b) colors hp=h WHILE hp<0##: hp=hp+2##*pi: WEND SELECT CASE hp MOD (2##*pi) CASE 0 r&=i&+2*i&*s: g&=i&-i&*s: b&=i&-i&*s CASE 0 TO 2##*pi/3## r&=i&+i&*s*COS(h)/COS(pi/3##-h) g&=i&+i&*s*(1##-COS(h)/COS(pi/3##-h)) b&=i&-i&*s CASE 2##*pi/3## TO 4##*pi/3## r&=i&-i&*s g&=i&+i&*s*COS(h-2##*pi/3##)/COS(pi-h) b&=i&+i&*s*(1##-COS(h-2##*pi/3##)/COS(pi-h)) CASE 4##*pi/3## TO 2##*pi r&=i&+i&*s*(1##-COS(h-4##*pi/3##)/COS(5##*pi/3##-h)) g&=i&-i&*s b&=i&+i&*s*COS(h-4##*pi/3##)/COS(5##*pi/3##-h##) END SELECT END SUB FUNCTION inelegance&(a$,np&) 'Calculates the inelegance of the parameters ie&=0 FOR i&=1 TO np& b$=MID$(a$,i&+1,1) SELECT CASE b$ CASE "A","B","X","Y","Z": ie&=ie&+4 CASE "C","W": ie&=ie&+1 CASE "M": ie&=ie& CASE ELSE: ie&=ie&+3 END SELECT NEXT i& FUNCTION=ie& END FUNCTION SUB koch(xp,yp,rp,c&,w&,h&,v&,th,d&) 'Makes a Koch curve centered at (xp,yp) with radius rp and v& vertices and orientation th filled with color c& 'd& is the number levels of recursion (0 for regular polygon) axiom$="F"+REPEAT$(v&-1,"+F") rule$=axiom$ fo=2##*rp*SIN(pi/v&) 'Length of each leg of polygon dth=2##*pi/v& 'Turning angle (45 degrees) IF d&=0 THEN fc=fo ELSE FOR k&=1 TO d& a$="F+F-F-FF+F+F-F" REPLACE "F" WITH a$ IN rule$ fc=fo/(4##+4##*COS(dth))^k& NEXT k& END IF x=xp+rp*COS(th) 'Starting position y=yp-rp*SIN(th) th=th-pi/2-pi/v& 'moving clockwise FOR k&=1 TO LEN(rule$) c$=MID$(rule$,k&,1) CALL turtle(c$,fc,dth,x,y,th,w&,h&) NEXT k& GRAPHIC PAINT BORDER (xp,yp),c&,%BLACK GRAPHIC REDRAW END SUB FUNCTION locle (x(),xe(),a(),ne&,dr,t$,h,xy(),w&,h&,n&&) 'Calculates the local largest Lyapunov exponent REDIM dx(ne&) n&=n&& MOD 2e9 map (n&,x(),a(),ne&,t$,h,xy(),w&,h&) map (n&,xe(),a(),ne&,t$,h,xy(),w&,h&) rs=0## FOR i&=1 TO ne& dx(i&)=xe(i&)-x(i&) rs=rs+dx(i&)*dx(i&) NEXT i& rs=SQR(rs)/dr IF rs=0## THEN EXIT FUNCTION 'Orbits collapsed together FOR i&=1 TO ne& xe(i&)=x(i&)-dx(i&)/rs NEXT i& locle=LOG(rs)/h END FUNCTION SUB makeicon(xp,yp,w&,h&,ns&) 'Make an icon with ns& segments xpo=xp ypo=yp i&=INT(ns&*ran) r=xpo/w& th=6.2832##*(ypo/(h&-9)+i&)/ns& IF (ns& MOD 2)=0 AND (i& MOD 2)=0 THEN th=6.2832##/ns&-th xp=(1##-r*SIN(th))*w&/2## yp=(1##+r*COS(th))*h&/2## END SUB SUB map(n&,x(),a(),ne&,t$,h,xy(),w&,h&) 'Advances the variables for the map h=1e-3## 'Step size for flow x=x(1): y=x(2): z=x(3): u=x(4): v=x(5): w=x(6) SELECT CASE t$ CASE"A": x(1)=a(1)*x+a(2)*x*x+a(3)*y+a(4)*y*y+a(5)*z+a(6)*z*z+a(7)*u+a(8)*u*u+a(9)*v+a(10)*v*v+a(11)*w+a(12)*w*w CASE"B": x(1)=a(1)*x+a(2)*x*x*x+a(3)*y+a(4)*y*y*y+a(5)*z+a(6)*z*z*z+a(7)*u+a(8)*u*u*u+a(9)*v+a(10)*v*v*v+a(11)*w+a(12)*w*w*w CASE"C": x(1)=x*(a(1)+a(2)*x*x+a(3)*y*y+a(4)*z*z)+y*(a(5)+a(6)*x*x+a(7)*y*y+a(8)*z*z)+z*(a(9)+a(10)*x*x+a(11)*y*y+a(12)*z*z) CASE"D": x(1)=a(1)*x+a(2)*SGN(x)+a(3)*y+a(4)*SGN(y)+a(5)*z+a(6)*SGN(z)+a(7)*u+a(8)*SGN(u)+a(9)*v+a(10)*SGN(v)+a(11)*w+a(12)*SGN(w) CASE"E": x(1)=a(1)*x+a(2)*SIN(a(3)*x)+a(4)*y+a(5)*SIN(a(6)*y)+a(7)*z+a(8)*SIN(a(9)*z)+a(10)*u+a(11)*SIN(a(12)*u) CASE"F": x(1)=a(1)*x+a(2)*TAN(a(3)*x)+a(4)*y+a(5)*TAN(a(6)*y)+a(7)*z+a(8)*TAN(a(9)*z)+a(10)*u+a(11)*TAN(a(12)*u) CASE"G": x(1)=a(1)*x+a(2)*sinh(a(3)*x)+a(4)*y+a(5)*sinh(a(6)*y)+a(7)*z+a(8)*sinh(a(9)*z)+a(10)*u+a(11)*sinh(a(12)*u) CASE"H": x(1)=a(1)*x+a(2)*tanh(a(3)*x)+a(4)*y+a(5)*tanh(a(6)*y)+a(7)*z+a(8)*tanh(a(9)*z)+a(10)*u+a(11)*tanh(a(12)*u) CASE"I": x(1)=a(1)*x+a(2)*CINT(x)+a(3)*y+a(4)*CINT(y)+a(5)*z+a(6)*CINT(z)+a(7)*u+a(8)*CINT(u)+a(9)*v+a(10)*CINT(v)+a(11)*w+a(12)*CINT(w) CASE"J" dx=(a(1)+a(2)*x+a(3)*y+a(4)*z)*y dy=(a(5)+a(6)*x+a(7)*y+a(8)*z)*z dz=(a(9)+a(10)*x+a(11)*y+a(12)*z)*x x(1)=x+h*dx x(2)=y+h*dy x(3)=z+h*dz x(4)=dx x(5)=dy x(6)=dz EXIT SUB CASE"K" dx=y dy=z dz=x*(a(1)+a(2)*x*x+a(3)*y*y+a(4)*z*z)+y*(a(5)+a(6)*x*x+a(7)*y*y+a(8)*z*z)+z*(a(9)+a(10)*x*x+a(11)*y*y+a(12)*z*z) x(1)=x+h*dx x(2)=y+h*dy x(3)=z+h*dz x(4)=dx x(5)=dy x(6)=dz EXIT SUB CASE"L" dx=y dy=z dz=u du=a(1)*x+a(2)*SIN(a(3)*x)+a(4)*y+a(5)*SIN(a(6)*y)+a(7)*z+a(8)*SIN(a(9)*z)+a(10)*u+a(11)*SIN(a(12)*u) x(1)=x+h*dx x(2)=y+h*dy x(3)=z+h*dz x(4)=u+h*du x(5)=dx x(6)=dy EXIT SUB CASE"M" dx=y dy=z dz=u du=a(1)*x+a(2)*tanh(a(3)*x)+a(4)*y+a(5)*tanh(a(6)*y)+a(7)*z+a(8)*tanh(a(9)*z)+a(10)*u+a(11)*tanh(a(12)*u) x(1)=x+h*dx x(2)=y+h*dy x(3)=z+h*dz x(4)=u+h*du x(5)=dx x(6)=dy EXIT SUB CASE"N" dx=a(1)*x+a(2)*y+a(3)*z dy=a(4)*x+a(5)*y+a(6)*z+a(7)*x*z dz=a(8)*x+a(9)*y+a(10)*z+a(11)*x*y+a(12)*y*y x(1)=x+h*dx x(2)=y+h*dy x(3)=z+h*dz x(4)=dx x(5)=dy x(6)=dz EXIT SUB CASE"O" SELECT CASE CEIL(4*ran) CASE 1: x(1)=a(1)+a(2)*x+a(3)*y CASE 2: x(1)=a(4)+a(5)*x+a(6)*y CASE 3: x(1)=a(7)+a(8)*x+a(9)*y CASE 4: x(1)=a(10)+a(11)*x+a(12)*y END SELECT CASE"P" SELECT CASE CEIL(3*ran) CASE 1: x(1)=a(1)+a(2)*x+a(3)*y+a(4)*z CASE 2: x(1)=a(5)+a(6)*x+a(7)*y+a(8)*z CASE 3: x(1)=a(9)+a(10)*x+a(11)*y+a(12)*z END SELECT CASE"Q" SELECT CASE CEIL(2*ran) CASE 1: x(1)=a(1)+a(2)*x+a(3)*y+a(4)*x*x+a(5)*x*y+a(6)*y*y CASE 2: x(1)=a(7)+a(8)*x+a(9)*y+a(10)*x*x+a(11)*x*y+a(12)*y*y END SELECT CASE"R" SELECT CASE CEIL(3*ran) CASE 1: x(1)=a(1)+a(2)*x*x*x+a(3)*y*y*y+a(4)*z*z*z CASE 2: x(1)=a(5)+a(6)*x*x*x+a(7)*y*y*y+a(8)*z*z*z CASE 3: x(1)=a(9)+a(10)*x*x*x+a(11)*y*y*y+a(12)*z*z*z END SELECT CASE"S" x(1)=a(1)+a(3)*x-a(4)*y x(2)=a(2)+a(4)*x+a(3)*y a=a(5): b=a(6) x(1)=x(1) + a*x*x - 2*b*x*y - a*y*y x(2)=x(2) + b*x*x + 2*a*x*y - b*y*y a=a(7): b=a(8) x(1)=x(1) + a*x*x*x - 3*b*x*x*y - 3*a*x*y*y + b*y*y*y x(2)=x(2) + b*x*x*x + 3*a*x*x*y - 3*b*x*y*y - a*y*y*y a=a(9): b=a(10) x(1)=x(1) + a*x^4 - 4*b*x*x*x*y - 6*a*x*x*y*y + 4*b*x*y*y*y + a*y^4 x(2)=x(2) + b*x^4 + 4*a*x*x*x*y - 6*b*x*x*y*y - 4*a*x*y*y*y + b*y^4 a=a(11): b=a(12) x(1)=x(1) + a*x^5 - 5*b*x*x*x*x*y - 10*a*x*x*x*y*y + 10*b*x*x*y*y*y + 5*a*x*y^4 - b*y^5 x(2)=x(2) + b*x^5 + 5*a*x*x*x*x*y - 10*b*x*x*x*y*y - 10*a*x*x*y*y*y + 5*b*x*y^4 + a*y^5 h=1## EXIT SUB CASE"T" rule???=rani???(a()) xi&=0 IF ROUND(a(1),1)<-0.8## THEN 'It's a 1-D CA IF n& 0: IF ran 0: IF ran 0: IF ran 0: IF ran1 THEN ITERATE x1=x1+a(j&+4)*xy(ij&) x2=x2+a(j&+7)*xy(ij&) NEXT j& SELECT CASE (12+CINT(10*a(1))) MOD 13 'Select map rule CASE 0: xi=a(2)*tanh(x1)+a(12)*tanh(x2) 'UA CASE 1: xi=a(2)*SIN(x1)+a(12)*SIN(x2) 'UB CASE 2: xi=a(2)*COS(x1)+a(12)*COS(x2) 'UC CASE 3: xi=a(2)*SIN(x1)+a(12)*COS(x2) 'UD CASE 5: xi=SIN(xi) 'UE CASE 6: xi=COS(xi) 'UF CASE 4: xi=1##-2##*(xi)^2 'UG CASE 7: xi=1##-2##*(xi)^4 'UH Higher powers and odd powers rarely give chaos CASE 8: xi=1##-2##*ABS(xi) 'UI CASE 9: xi=1##-((2##*xi+1##) MOD 2##) 'UJ CASE 10: xi=1##-(4##*(xi+1##) MOD 2##) 'UK CASE 11: xi=1##-(6##*(xi+1##) MOD 2##) 'UL CASE 12: xi=1##-(8##*(xi+1##) MOD 2##) 'UM END SELECT END IF IF (12+CINT(10*a(1)))>12 THEN 'Modification for flow lattice ij&=(w&+n&) MOD (2*w&) xi=xy(ij&)+h*xi END IF x(1)=xi 'Update current point ij&=n& MOD (2*w&) xy(ij&)=xi xp&=n& MOD w& 'and plot it yp&=INT(n&/w&) MOD h& hue=0.5##*pi*(2.2##+a(10)+xi) END IF hsi(hue,1,128,r&,g&,b&) c&=RGB(r&,g&,b&) 'in this color IF a(11)<=0 AND xi<0## THEN c&=%WHITE GRAPHIC SET PIXEL(xp&,yp&),c& h=1## EXIT SUB CASE"V" FOR i&=1 TO ne& dx=2##*ran-1## x(i&)=x(i&)+h*dx IF x(i&)<-1## THEN x(i&)=x(i&)+2## IF x(i&)>1## THEN x(i&)=x(i&)-2## NEXT i& h=1## EXIT SUB CASE"W" IF n&0 IF n&=w&*h& THEN r&=MIN(h&,w&)/2##-1 'GRAPHIC ELLIPSE(w&/2+r&-2,h&/2+r&-2)-(w&/2-r&+1,h&/2-r&+1),c& FOR th=0 TO 2##*pi STEP 1##/r& x&=w&/2##+r&*COS(th) y&=h&/2##+r&*SIN(th) xy(x&+w&*y&)=1## NEXT th END IF x&=w&/2: y&=h&/2 END SELECT GRAPHIC REDRAW DO SELECT CASE INT(4*ran) CASE 0: INCR x&: IF x&>w& THEN x&=x&-w& CASE 1: DECR y&: IF y&<2 THEN y&=y&+h& CASE 2: DECR x&: IF x&<2 THEN x&=x&+w& CASE 3: INCR y&: IF y&>h& THEN y&=y&-h& END SELECT ind&=x&+w&*y& IF ind&1.3##+a(11) THEN EXIT IF xy(ind&)=1## IF (n& MOD 10)=0 THEN GRAPHIC REDRAW x(1)=2##*x&/w&-1## x(2)=2##*y&/h&-1## x(3)=1## x(4)=r&/128-1 x(5)=g&/128-1 x(6)=b&/128-1 SELECT CASE CINT(a(12)) 'See if we are done CASE 0: IF ABS(x(1))>0.99## OR ABS(x(2))>=0.99## THEN n&=10*w&*h& CASE ELSE: IF ABS(x(1))<0.01## AND ABS(x(2))<=0.01## THEN n&=10*w&*h& END SELECT h=1## EXIT SUB END IF LOOP CASE"X" IF n&=0 THEN 'Initialize FOR j&=0 TO h&-1 FOR i&=0 TO w&-1 IF ran>(a(12)+1.3) THEN ITERATE nc&=2+CINT(12##+10##*a(11)) nc&=INT(nc&*ran) pickcolor(nc&,128,r&,g&,b&) GRAPHIC SET PIXEL(i&,j&),RGB(r&,g&,b&) NEXT j& GRAPHIC REDRAW NEXT i& END IF FOR k&=1 TO 100*(a(9)+1.3##) r=2.3##+a(10) 'The neighborhood size io&=INT(w&*ran) 'Pick a random point jo&=INT(h&*ran) DO di&=2*r*ran-r 'and a random neighbor dj&=2*r*ran-r LOOP UNTIL di&*di&+dj&*djo&<=r*r 'within a circle of radius r& i&=(w&+io&+di&) MOD w& j&=(h&+jo&+dj&) MOD h& GRAPHIC GET PIXEL(i&,j&) TO c& GRAPHIC SET PIXEL(io&,jo&),c& 'and make the replacement NEXT k& h=1## EXIT SUB CASE"Y" IF n&=0 THEN 'Initialize c&=RGB(120##+100##*a(9),120##+100##*a(10),120##+100##*a(11)) FOR j&=0 TO h&-1 FOR i&=0 TO w&-1 IF ran>0.6##+0.02##*a(12) THEN GRAPHIC SET PIXEL(i&,j&),c& NEXT i& NEXT j& GRAPHIC REDRAW END IF IF n&>=w& THEN n&=2*n& ELSE GRAPHIC GET PIXEL(n&,0) TO b& IF b&<>%WHITE THEN EXIT IF r&=1+(n& MOD 6) GRAPHIC PAINT REPLACE (n&,0),c&(r&),b& GRAPHIC REDRAW END IF CASE"Z" z&=13##+10##*a(1) SELECT CASE z& CASE 1 IF n&=0 THEN GRAPHIC BOX(0,0)-(w&-1,h&-1) f=(13##+10##*a(9))/54## x=f FOR i&=1 TO f*w&*h& IF ran<0.5## THEN x=f*x ELSE x=1##-f*x GRAPHIC LINE(0,x*h&)-(w&,x*h&) GRAPHIC LINE(x*w&,0)-(x*w&,h&) NEXT i& END IF CASE 2 IF n&=0 THEN GRAPHIC BOX(0,0)-(w&-1,h&-1) x=0## FOR i&=1 TO 130+100*a(9) f=ran/2## IF ran<0.5## THEN x=f*x ELSE x=1##-f*x GRAPHIC LINE(0,x*h&)-(w&,x*h&) GRAPHIC LINE(x*w&,0)-(x*w&,h&) NEXT i& END IF CASE 3 IF n&=0 THEN GRAPHIC BOX(0,0)-(w&-1,h&-1) x=0## FOR i&=1 TO 130+100*a(9) f=ran/2## IF ran<0.5## THEN x=f*x ELSE x=1##-f*x GRAPHIC LINE(0,x*h&)-(w&,x*h&) f=ran/2## IF ran<0.5## THEN x=f*x ELSE x=1##-f*x GRAPHIC LINE(x*w&,0)-(x*w&,h&) NEXT i& END IF CASE 4 IF n&=0 THEN GRAPHIC BOX(0,0)-(w&-1,h&-1) FOR i&=1 TO 130+100*a(9) x=ran GRAPHIC LINE(0,x*h&)-(w&,x*h&) x=ran GRAPHIC LINE(x*w&,0)-(x*w&,h&) NEXT i& END IF CASE 5 IF n&=0 THEN GRAPHIC BOX(0,0)-(w&-1,h&-1) FOR i&=1 TO 130+100*a(9) x=ran y=ran th=pi*ran xp=x*w& yp=y*h& DO xp=xp+COS(th) yp=yp+SIN(th) GRAPHIC SET PIXEL(xp,yp) LOOP UNTIL xp<=0 OR yp<=0 OR xp>=w& OR yp>=h& xp=x*w& yp=y*h& DO xp=xp-COS(th) yp=yp-SIN(th) GRAPHIC SET PIXEL(xp,yp) LOOP UNTIL xp<=0 OR yp<=0 OR xp>=w& OR yp>=h& NEXT i& END IF CASE 6 IF n&=0 THEN GRAPHIC BOX(0,0)-(w&-1,h&-1) FOR i&=1 TO 130+100*a(9) x=ran y=ran th=pi*ran xp=x*w& yp=y*h& DO xp=xp+COS(th) yp=yp+SIN(th) GRAPHIC GET PIXEL(xp,yp) TO c& LOOP UNTIL c&<>%WHITE xo=xp yo=yp xp=x*w& yp=y*h& DO xp=xp-COS(th) yp=yp-SIN(th) GRAPHIC GET PIXEL(xp,yp) TO c& LOOP UNTIL c&<>%WHITE GRAPHIC LINE(xo,yo)-(xp,yp) NEXT i& END IF CASE 7 IF n&=0 THEN GRAPHIC BOX(0,0)-(w&-1,h&-1) FOR i&=1 TO 130+100*a(9) x=ran y=ran IF ran<0.5## THEN th=0## ELSE th=pi/2## xp=x*w& yp=y*h& DO xp=xp+COS(th) yp=yp+SIN(th) GRAPHIC GET PIXEL(xp,yp) TO c& LOOP UNTIL c&<>%WHITE xo=xp yo=yp xp=x*w& yp=y*h& DO xp=xp-COS(th) yp=yp-SIN(th) GRAPHIC GET PIXEL(xp,yp) TO c& LOOP UNTIL c&<>%WHITE GRAPHIC LINE(xo,yo)-(xp,yp) NEXT i& END IF CASE 8 IF n&=0 THEN b=(1.3##+a(9))/5.4## t=0## told=0## xm2=1## t1=ATN(1##/b) ym2=EXP(-b*t1)*SIN(t1) t2=pi+ATN(-b) xm1=EXP(-b*t2)*COS(t2) t3=pi+t1 ym1=EXP(-b*t3)*SIN(t3) DO r=EXP(-b*t) x=r*COS(t) y=r*SIN(t) xp&=(x-xm1)*(w&-3)/(xm2-xm1)+1 yp&=(ym2-y)*(h&-3)/(ym2-ym1)+1 GRAPHIC SET PIXEL(xp&,yp&) IF r<2##/MAX(w&,h&) THEN EXIT LOOP IF told<=INT(t) AND t>INT(t) THEN tt=t+3.9##+3##*a(10) r0=EXP(-b*tt) x=r0*COS(tt) y=r0*SIN(tt) xpp&=(x-xm1)*(w&-3)/(xm2-xm1)+1 ypp&=(ym2-y)*(h&-3)/(ym2-ym1)+1 GRAPHIC LINE(xpp&,ypp&)-(xp&,yp&) END IF told=t t=t+h/(r*MAX(w&,h&)) LOOP END IF CASE 9 IF n&=0 THEN sc&=INT(6##*ran) IF ran>0.5## THEN bc&=%WHITE ELSE bc&=c&(7) GRAPHIC BOX(0,0)-(w&-1,h&-1),,bc&,bc& NN=3.2+a(10) c=1.26##+0.2##*a(11) imax&=700+500*a(9) FOR i&=0 TO imax& ai=1##/(zeta(c,NN)*(NN+i&)^c) r=SQR(ai/pi) kmax&=1e4 'Allow ten thousand trials FOR k&=1 TO kmax& x=ran y=ran IF x-r>0## AND y-r>0## AND x+r<1## AND y+r<1## THEN EXIT FOR NEXT k& xp=x*w& yp=y*h& rx=r*w& ry=r*h& IF rx<1## OR ry<1## THEN EXIT FOR fc&=c&(1+((i&+sc&) MOD 6)) GRAPHIC ELLIPSE(xp-rx,yp-ry)-(xp+rx,yp+ry),%BLACK,fc& GRAPHIC REDRAW NEXT i& END IF CASE 10 IF n&=0 THEN sc&=INT(6##*ran) IF ran>0.5## THEN bc&=%WHITE ELSE bc&=c&(7) GRAPHIC BOX(0,0)-(w&-1,h&-1),,bc&,bc& NN=3.2+a(10) c=1.26##+0.2##*a(11) imax&=700+500*a(9) FOR i&=0 TO imax& ai=1##/(zeta(c,NN)*(NN+i&)^c) r=SQR(ai/4) kmax&=1e4 'Allow ten thousand trials FOR k&=1 TO kmax& x=ran y=ran IF x-r>0## AND y-r>0## AND x+r<1## AND y+r<1## THEN EXIT FOR NEXT k& xp=x*w& yp=y*h& rx=r*w& ry=r*h& IF rx<1## OR ry<1## THEN EXIT FOR fc&=c&(1+((i&+sc&) MOD 6)) GRAPHIC BOX(xp-rx,yp-ry)-(xp+rx,yp+ry),,%BLACK,fc& GRAPHIC REDRAW NEXT i& END IF CASE 11 IF n&=0 THEN sc&=INT(6##*ran) IF ran>0.5## THEN bc&=%WHITE ELSE bc&=c&(7) GRAPHIC BOX(0,0)-(w&-1,h&-1),,bc&,bc& NN=3.2+a(10) c=1.26##+0.2##*a(11) imax&=700+500*a(9) DIM xj(imax&),yj(imax&),rj(imax&) FOR i&=0 TO imax& ai=1##/(zeta(c,NN)*(NN+i&)^c) r=SQR(ai/pi) kmax&=1e7 'Allow ten million trials FOR k&=1 TO kmax& x=ran y=ran ol&=0 FOR j&=0 TO i&-1 'test for overlap IF i&=0 THEN EXIT FOR dr=r+rj(j&) dx=x-xj(j&) IF dx>dr THEN ITERATE dy=y-yj(j&) IF dx*dx+dy*dy0## AND y-r>0## AND x+r<1## AND y+r<1## AND ol&=0 THEN EXIT FOR NEXT k& IF k&>=kmax& THEN EXIT FOR 'Assume the search has halted xp=x*w& yp=y*h& rx=r*w& ry=r*h& IF rx<1## OR ry<1## THEN EXIT FOR fc&=c&(1+((i&+sc&) MOD 6)) GRAPHIC ELLIPSE(xp-rx,yp-ry)-(xp+rx,yp+ry),%BLACK,fc& GRAPHIC REDRAW xj(i&)=x yj(i&)=y rj(i&)=r NEXT i& END IF CASE 12 IF n&=0 THEN sc&=INT(6##*ran) IF ran>0.5## THEN bc&=%WHITE ELSE bc&=c&(7) GRAPHIC BOX(0,0)-(w&-1,h&-1),,bc&,bc& NN=3.2+a(10) c=1.26##+0.2##*a(11) imax&=700+500*a(9) DIM xj(imax&),yj(imax&),rj(imax&) FOR i&=0 TO imax& ai=1##/(zeta(c,NN)*(NN+i&)^c) r=SQR(ai/4) kmax&=1e7 'Allow ten million trials FOR k&=1 TO kmax& x=ran y=ran ol&=0 FOR j&=0 TO i&-1 'test for overlap IF i&=0 THEN EXIT FOR dr=r+rj(j&) dx=ABS(x-xj(j&)) dy=ABS(y-yj(j&)) IF dx0## AND y-r>0## AND x+r<1## AND y+r<1## AND ol&=0 THEN EXIT FOR NEXT k& IF k&>=kmax& THEN EXIT FOR 'Assume the search has halted xp=x*w& yp=y*h& rx=r*w& ry=r*h& IF rx<1## OR ry<1## THEN EXIT FOR fc&=c&(1+((i&+sc&) MOD 6)) GRAPHIC BOX(xp-rx,yp-ry)-(xp+rx,yp+ry),,%BLACK,fc& GRAPHIC REDRAW xj(i&)=x yj(i&)=y rj(i&)=r NEXT i& END IF CASE ELSE IF n&=0 THEN sc&=INT(6##*ran) IF ran>0.5## THEN bc&=%WHITE ELSE bc&=c&(7) GRAPHIC BOX(0,0)-(w&-1,h&-1),,bc&,bc& NN=3.2+a(10) c=1.26##+0.2##*a(11) imax&=700+500*a(9) DIM xj(imax&),yj(imax&),rj(imax&) FOR i&=0 TO imax& ai=w&*h&/(zeta(c,NN)*(NN+i&)^c) rp=SQR(ai/pi) kmax&=1e7 'Allow ten million trials FOR k&=1 TO kmax& xp=ran*w& yp=ran*h& ol&=0 FOR j&=0 TO i&-1 'test for overlap IF i&=0 THEN EXIT FOR dr=rp+rj(j&) dx=xp-xj(j&) IF dx>dr THEN ITERATE dy=yp-yj(j&) IF dx*dx+dy*dy0## AND yp-rp>0## AND xp+rp=kmax& THEN EXIT FOR 'Assume the search has halted IF rp<1## THEN EXIT FOR fc&=c&(1+((i&+sc&) MOD 6)) SELECT CASE z& 'Upright pentagon CASE 13: koch(xp,yp,rp,fc&,w&,h&,5,pi/2##,0) CASE 14: koch(xp,yp,rp,fc&,w&,h&,6,0##,0) CASE 15: koch(xp,yp,rp,fc&,w&,h&,7,pi/2##,0) CASE 16: koch(xp,yp,rp,fc&,w&,h&,8,pi/8##,0) CASE 17: v&=5: koch(xp,yp,rp,fc&,w&,h&,v&,2##*pi*ran/v&,0) CASE 18: v&=6: koch(xp,yp,rp,fc&,w&,h&,v&,2##*pi*ran/v&,0) CASE 19: v&=7: koch(xp,yp,rp,fc&,w&,h&,v&,2##*pi*ran/v&,0) CASE 20: v&=8: koch(xp,yp,rp,fc&,w&,h&,v&,2##*pi*ran/v&,0) CASE 21: v&=5+INT(4##*ran): koch(xp,yp,rp,fc&,w&,h&,v&,2##*pi*ran/v&,0) CASE 22: v&=8: koch(xp,yp,rp,fc&,w&,h&,v&,2##*pi*ran/v&,5) CASE ELSE koch(xp,yp,rp,fc&,w&,h&,6,pi/6##,0) xq=xp+rp*SIN(pi/3##) yq=yp-rp*COS(pi/3##) GRAPHIC LINE(xp,yp)-(xq,yq) xq=xp-rp*SIN(pi/3##) yq=yp-rp*COS(pi/3##) GRAPHIC LINE(xp,yp)-(xq,yq) GRAPHIC PAINT BORDER (xp+1,yp+1),%LTGRAY,%BLACK xq=xp yq=yp+rp GRAPHIC LINE(xp,yp)-(xq,yq) GRAPHIC PAINT BORDER (xp-1,yp+1),%GRAY,%BLACK END SELECT GRAPHIC REDRAW xj(i&)=xp yj(i&)=yp rj(i&)=rp NEXT i& END IF END SELECT n&=2*n& GRAPHIC REDRAW END SELECT x(2)=x x(3)=y x(4)=z x(5)=u x(6)=v h=1## END SUB SUB modify(w&,h&,nd&,ipp&,ni&,p&,bail&) 'Allow one to modify the default parameters CLS LINE INPUT"Width of the plot in pixels (default 1200): ",q$ IF q$="" THEN w&=1200 ELSE w&=VAL(q$) LINE INPUT"Height of the plot in pixels (default 900): ",q$ IF q$="" THEN h&=900 ELSE h&=VAL(q$) LINE INPUT"Number of digits in the printout (default 4): ",q$ IF q$="" THEN nd&=4 ELSE nd&=VAL(q$) LINE INPUT"Number of iterations per pixel (default 10): ",q$ IF q$="" THEN ipp&=10 ELSE ipp&=VAL(q$) LINE INPUT"Number of iterations between printouts (default 1e5): ",q$ IF q$="" THEN ni&=1e5 ELSE ni&=VAL(q$) LINE INPUT"Maximum period for periodicity checking (default 1e4): ",q$ IF q$="" THEN mp&=1e4 ELSE mp&=VAL(q$) LINE INPUT"Bailout value for escape time plots (default 1e3): ",q$ IF q$="" THEN bail&=1e3 ELSE bail&=VAL(q$) CLS END SUB SUB pickcolor(nc&,i&,r&,g&,b&) 'Finds the rgb value for a color with index nc& and intensity i& (0 to 255) 'Choose i&=128 for normal intensity h=1.9416##*(nc&-1) 'hue (pi divided by the Golden mean) s=1## 'saturation CALL hsi(h,s,i&,r&,g&,b&) END SUB FUNCTION pixfrac(w&,h&) 'Counts the fraction of pixels in the imaage ignoring the shadow g&=c&(7) 'Shadow intensity FOR i&=0 TO w&-1 FOR j&=0 TO h&-1 GRAPHIC GET PIXEL(i&,j&) TO c& IF c&=%WHITE OR c&=g& THEN ITERATE INCR pf& NEXT j& NEXT i& FUNCTION = pf&/(w&*h&) END FUNCTION FUNCTION ran STATIC '$1000 if ran isn't random over [0,1) 'Ref: W. A. Press and S. A. Teukolsky, Computers in Physics 6, 522 (1992) 'This routine is about 5X slower than the BASIC RND function 'It reseeds itself with the global variable seed& whenever seed& is negative IF seed& <= 0 THEN 'Initialize on first call im1& = 2147483563 im2& = 2147483399 am = 1## / im1& imm1& = im1& - 1 ia1& = 40014 ia2& = 40692 iq1& = 53668 iq2& = 52774 ir1& = 12211 ir2& = 3791 ntab& = 32 ndiv& = 1 + imm1& \ ntab& DIM iv&(ntab&) idum2& = 123456789 FOR j& = 1 TO ntab&: iv&(j&) = 0: NEXT j& iy& = 0 seed& = MAX (- seed&, 1) idum2& = seed& FOR j& = ntab& + 8 TO 1 STEP -1 k& = seed& \ iq1& seed& = ia1& * (seed& - k& * iq1&) - k& * ir1& IF seed& < 0 THEN seed& = seed& + im1& IF j& <= ntab& THEN iv&(j&) = seed& NEXT j& y& = iv&(1) END IF k& = seed& \ iq1& seed& = ia1& * (seed& - k& * iq1&) - k& * ir1& IF seed& < 0 THEN seed& = seed& + im1& k& = idum2& \ iq2& idum2& = ia2& * (idum2& - k& * iq2&) - k& * ir2& IF idum2& < 0 THEN idum2& = idum2& + im2& j& = 1 + iy& \ ndiv& iy& = iv&(j&) - idum2& iv&(j&) = seed& IF iy& < 1 THEN iy& = iy& + imm1& ran = am * iy& END FUNCTION FUNCTION rani???(a()) 'Generates an integer in the range (0,4294967296) based on the parameters 2 through 8 '(useful to seed the portable random number generator) r&&=0 FOR i&=0 TO 6 r&&=r&&+(12##+10##*a(i&+2))*26^i& NEXT i& WHILE r&&>=4294967296 'This is the modulo function r&&=r&&-4294967296 WEND FUNCTION=r&& END FUNCTION SUB setpixel(xp,yp,zp,upp,vp,wp) 'Plots a point at (xp,yp) at depth zp with anti-aliasing x&=INT(xp) y&=INT(yp) dx=xp-x& dy=yp-y& FOR i&=0 TO 1 FOR j&=0 TO 1 g&=255-255*zp*(i&-dx)*(j&-dy) red=g&*(1-wp) grn=g&*(1-vp) blu=g&*(1-upp) fg&=RGB(red,grn,blu) GRAPHIC SET PIXEL(x&+i&,y&+j&),fg& NEXT j& NEXT i& END SUB SUB shadow(xp,yp,zp,w&,h&,h) 'Adds a shadow IF h<1## THEN IF RAN*LOG10(h)<-1## THEN EXIT SUB 'Don't over-shadow a flow g&=c&(7) 'Shadow intensity gs=0: FOR i&=1 TO 12: gs=gs+RAN: NEXT i& 'Approximation to a Gaussian xs=xp+(w&/5)*zp/gs gs=0: FOR i&=1 TO 12: gs=gs+RAN: NEXT i& 'Approximation to a Gaussian ys=yp+(h&/5)*zp/gs GRAPHIC GET PIXEL(xs,ys) TO s& IF s&=%RGB_WHITE THEN 'Make shadow GRAPHIC SET PIXEL(xs,ys),g& END IF END SUB FUNCTION sinh(x) 'Returns hyperbolic sine of x sinh = (EXP(x) - EXP(-x)) / 2## END FUNCTION FUNCTION tanh(x) 'Returns hyperbolic tangent of x IF ABS(x)>22 THEN tanh = SGN(x) ELSE tanh = 1## - 2## / (EXP(2## * x) + 1##) END IF END FUNCTION SUB turtle(c$,fc,dth,x,y,th,w&,h&) 'Interprets turtle commands SELECT CASE c$ CASE "F" 'Move forward a distance fc while drawing a line GRAPHIC SET POS(x,y) x=x+fc*COS(th) y=y-fc*SIN(th) GRAPHIC LINE-(x,y) CASE "f" 'Move forward a distance fc while not drawing a line x=x+fc*COS(th) y=y-fc*SIN(th) CASE "+" 'Turn right through angle dth th=th-dth CASE "-" 'Turn left through angle dth th=th+dth END SELECT END SUB FUNCTION printwrap$(a$,l&) 'Prints a string with words wrapped to a maximum line length of l& DO a$=TRIM$(a$) IF LEN(a$)