'Program SYMMETRY.BAS (c) 1997 by J. C. Sprott 'Searches for 3-D strange attractors from ODEs with cyclic symmetry defext a-z 'Use extended (80-bit) precision throughout nd&=5000 'Number of initial iterations to discard h=.02 'Iteration step size dim a(30) as shared ext screen 12 randomize timer 'Reseed the random number generator t0=timer do x=.1*rnd: y=.1*rnd: z=.1*rnd 'Initial conditions (each run is different) call getcode(command\$,code\$) 'Get code for case to analyze call coeff(code\$,a()) 'Get equation coefficients in array a(30) xmin=1e37: ymin=xmin: zmin=xmin: xmax=-xmin: ymax=-ymin: zmax=-zmin xlim=1e4: ylim=1e4: zlim=1e4 do if inkey\$=chr\$(27) then end 'Loop until user presses the key if n&>nd&/10 and n&xlim or abs(y)>ylim or abs(z)>zlim then exit loop if n&>nd& then 'Assume transient has settled dx=16*(z-zav)/zmm xp=213+427*(x-xmin)/xmm-213*(z-zmin)/zmm yp=320-320*(y-ymin)/ymm+160*(z-zmin)/zmm pset(xp-dx,yp),12: pset(xp+dx,yp),11 'Make stereo cyan/red plot call lyapunov(x,y,z,h,l1,xe) end if n&=n&+1 if n&=100*nd& then locate 1,1: print code\$ call savedata(code\$,d0,dd0,d2,dd2,l1,l2,l3,n&) call scn2bmp(left\$(code\$,8)+".BMP") end if loop until (n&>2*nd& and (l1>10 or l1<.005)) or n&>=100*nd& n&=0: xe=0: d0=0: dd0=0: d2=0: dd2=0 loop end sub rk4(x,y,z,h) 'Advance (x,y,z) with fourth order Runge-Kutta k1x=h*fnx(x,y,z) k1y=h*fny(x,y,z) k1z=h*fnz(x,y,z) k2x=h*fnx(x+.5*k1x,y+.5*k1y,z+.5*k1z) k2y=h*fny(x+.5*k1x,y+.5*k1y,z+.5*k1z) k2z=h*fnz(x+.5*k1x,y+.5*k1y,z+.5*k1z) k3x=h*fnx(x+.5*k2x,y+.5*k2y,z+.5*k2z) k3y=h*fny(x+.5*k2x,y+.5*k2y,z+.5*k2z) k3z=h*fnz(x+.5*k2x,y+.5*k2y,z+.5*k2z) k4x=h*fnx(x+k3x,y+k3y,z+k3z) k4y=h*fny(x+k3x,y+k3y,z+k3z) k4z=h*fnz(x+k3x,y+k3y,z+k3z) x=x+(k1x+2*(k2x+k3x)+k4x)/6 y=y+(k1y+2*(k2y+k3y)+k4y)/6 z=z+(k1z+2*(k2z+k3z)+k4z)/6 end sub def fnx(x,y,z)=a(1)+a(2)*x+a(3)*x*x+a(4)*x*y+a(5)*x*z+a(6)*y+a(7)*y*y+a(8)*y*z+a(9)*z+a(10)*z*z def fny(x,y,z)=a(1)+a(2)*y+a(3)*y*y+a(4)*y*z+a(5)*y*x+a(6)*z+a(7)*z*z+a(8)*z*x+a(9)*x+a(10)*x*x def fnz(x,y,z)=a(1)+a(2)*z+a(3)*z*z+a(4)*z*x+a(5)*z*y+a(6)*x+a(7)*x*x+a(8)*x*y+a(9)*y+a(10)*y*y sub getcode(c\$,code\$) 'Get code for case to analyze if len(c\$)=31 then code\$=c\$: exit sub code\$="I" for i%=1 to 10 a%=-12+int(25*rnd) code\$=code\$+chr\$(77+a%) next i% end sub sub coeff(code\$,a(1)) 'Return values of 30 equation coefficients for i%=1 to 10 a(i%)=(asc(mid\$(code\$,i%+1))-77)/10 next i% end sub sub minmax(x,y,z,xmin,xmax,ymin,ymax,zmin,zmax) 'Find min and max of (x,y,z) xmin=min(x,xmin) xmax=max(x,xmax) ymin=min(y,ymin) ymax=max(y,ymax) zmin=min(z,zmin) zmax=max(z,zmax) end sub sub lyapunov(x,y,z,h,l1,xe) 'Get largest Lyapunov exponent static ye,ze,lsum,nl&,j% if xe=0 then xe=x+1e-8##: ye=y: ze=z: lsum=0: nl&=0: j%=0: exit sub j%=(j%+1) mod 499 call rk4(xe,ye,ze,h) dlx=xe-x: dly=ye-y: dlz=ze-z dl2=dlx*dlx+dly*dly+dlz*dlz if cdbl(dl2)>0 then df=1e16##*dl2 rs=1##/sqr(df) xe=x+rs*(xe-x): ye=y+rs*(ye-y): ze=z+rs*(ze-z) lsum=lsum+log(df) nl&=nl&+1 end if if j%=0 then 'Don't print results too often l1=.5*lsum/nl&/abs(h) if abs(l1)<10 then locate 1,45: print" L =";using"##.###";l1; end if end sub sub savedata(code\$,d0,dd0,d2,dd2,l1,l2,l3,n&) 'Save data to NEW.DIC open"NEW.DIC" for append as #1 print#1,using"##.### 0 ##.###";l1;l3 close#1 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 x1% = 0: y1% = 0: x2% = 639: y2% = 479 'Corners of image x2% = x1% + 16 * ceil((x2% - x1%) \ 16 + .01) - 1 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\$("FF FF FF 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\$("00 00 00 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% 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