'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 <Esc> key
        if n&>nd&/10 and n&<nd& then call minmax(x,y,z,xmin,xmax,ymin,ymax,zmin,zmax)
        if n&=nd& then              'Set up the screen
            cls
            t1=timer
            if t1<t0 then t0=t0-24*3600
            incr nt&: locate 1,1: print clng(3600!*nt&/(t1-t0));
            xlim=3*max(abs(xmax),abs(xmin))
            ylim=3*max(abs(ymax),abs(ymin))
            zlim=3*max(abs(zmax),abs(zmin))
            xmm=xmax-xmin
            ymm=ymax-ymin
            zmm=zmax-zmin
            zav=(zmax+zmin)/2
            line(213,320)-(639,320),,,&H4444: locate 22,80: print"X";
            line(213,320)-(213,0),,,&H4444: locate 2,25: print"Y";
            line(213,320)-(0,479),,,&H4444: locate 29,1: print"Z";
        end if
        call rk4(x,y,z,h)           'Advance (x,y,z) by time step h
        if abs(x)>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