Code:
/*This was a copy and paste from the book, with a few of my alterations. Sorry if it is a little cryptic. I'm still trying to understand it*/
#include <dos.h> /*this came with MinGW and seems not to work*/
#include <stdio.h>
#include <stdlib.h>
#include <graphics.h> /*This failed to replace the 16-bit graph.h header*/
#include <conio[C].h> /*this is the c version of conio.h header*/
#include <wincon.h> //This, along with winbgim.h, was to help replace the old
#include <math.h> // dos.h header. It failed
#include <io.h>
#include <string.h>
#include <winbgim.h>
int PREV, OMAX, D, ODE, SND, PJT, TRD, FTH;
int SAV, T, WID, QM, P, TWOD;
int M, I, I1, I2, O, I3, I4, I5;
int J, WH, FREQ, C4, NC, C, RD, BY, BK;
int COLR[16];
char CODE[515], Q;
char FAV[9] = "XDATA.DAT";
double NMAX, EPS, TWOPI, SEG, NE, X;
double Yy, Z, W, XE, YE, ZE, WE, LSUM, N, NL;
double N1, N2, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX, WMIN, WMAX, XNEW, YNEW;
double ZNEW, WNEW, XP, YP, RAN, XSAVE, YSAVE, ZSAVE, WSAVE, DLX, DLY, DLZ;
double DLW, DL2, DF, RS, L, MX, MY, XL, XH, YL, YH, XA, YA, ZA, TT, PT, TIA;
double XZ, YZ, DUR, D2MAX, DX, DY, DZ, DW, D2, F, TH, PH, XRT, XLT, HSF, AL;
double SINAL, COSAL, DUM, SW, SH, X0, Y0;
double XS[500], YS[500], ZS[500], WS[500], A[505], V[100], XY[5], XN[5];
union REGS regs;
FILE *F1, *F2, *F3;
void fun1300(void); /* Initialize */
void fun1500(void); /* Set parameters */
void fun1700(void); /* Iterate equations */
void fun2100(void); /* Display results */
void fun2400(void); /* Test results */
void fun2600(void); /* Get coefficients */
void fun2800(void); /* Shuffle random numbers */
void fun2900(void); /* Calculate Lyapunov exponent */
void fun3100(void); /* Resize the screen */
void fun3500(void); /* Produce sound */
void fun3600(void); /* Respond to user command */
void fun3900(void); /* Calculate fractal dimension */
void fun4100(void); /* Project onto a sphere */
void fun4200(void); /* Display menu screen */
void fun4700(void); /* Get dimension and order */
void fun4900(void); /* Save attractor to disk file SA.DIC */
void fun5000(void); /* Plot point on screen */
void fun5400(void); /* Plot background grid */
void fun5600(void); /* Set colors */
void fun5800(void); /* Process evaluation command */
void fun5900(void); /* Save favorite attractors to disk file FAVORITE.DIC */
void fun6000(void); /* Update SA.DIC file */
void fun6200(void); /* Special function definitions */
void fun6700(void); /* Project onto a horizontal cylinder */
void fun6800(void); /* Project onto a vertical cylinder */
void fun6900(void); /* Project onto a torus (unity aspect ratio) */
void fun7000(void); /* Save data */
void main()
{
PREV = 5; /* Plot versus fifth previous iterate */
NMAX = 11000; /* Maximum number of iterations */
OMAX = 5; /* Maximum order of polynomial */
D = 2; /* Dimension of system */
EPS = .1; /* Step size for ODE */
ODE = 0; /* System is map */
SND = 0; /* Turn sound off */
PJT = 0; /* Projection is planar */
TRD = 1; /* Display third dimension as shadow */
FTH = 2; /* Display fourth dimension as colors */
SAV = 0; /* Don't save any data */
TWOPI = 6.28318530717959; /* A useful constant (2 pi) */
srand(time(0)); /* Reseed random number generator */
fun4200(); /* Display menu screen */
T = 1;
if (Q == 'X') T = 0; /* Exit immediately on command */
while (T) {
switch (T) {
case 1: fun1300(); /* Initialize */
case 2: fun1500(); /* Set parameters */
case 3: fun1700(); /* Iterate equations */
case 4: fun2100(); /* Display results */
case 5: fun2400(); /* Test results */
}
}
closegraph();
_setcursortype(_NORMALCURSOR);
clrscr();
return(0);
}
void fun1300() /* Initialize */
{
int graphdriver = DETECT, graphmode, error_code;
initgraph(&graphdriver, &graphmode, ""); /* Set graphics mode */
WID = 80; /* Number of text columns */
cleardevice();
gotoxy(WID / 2 - 6, 13);
printf("Searching...");
fun5600(); /* Set colors */
if (QM == 2) {
NE = 0;
fclose(F1);
F1 = fopen("SA.DIC", "a");
fclose(F1);
F1 = fopen("SA.DIC", "r");
}
}
void fun1500() /* Set parameters */
{
X = .05; /* Initial condition */
Yy = .05;
Z = .05;
W = .05;
XE = X + .000001;
YE = Yy;
ZE = Z;
WE = W;
fun2600(); /* Get coefficients */
T = 3;
P = 0;
LSUM = 0;
N = 0;
NL = 0;
N1 = 0;
N2 = 0;
XMIN = 1000000;
XMAX = -XMIN;
YMIN = XMIN;
YMAX = XMAX;
ZMIN = XMIN;
ZMAX = XMAX;
WMIN = XMIN;
WMAX = XMAX;
TWOD = _rotl(1, D);
}
void fun1700() /* Iterate equations */
{
if (ODE > 1)
fun6200(); /* Special function */
else {
M = 1;
XY[1] = X;
XY[2] = Yy;
XY[3] = Z;
XY[4] = W;
for (I = 1; I <= D; I++) {
XN[I] = A[M];
M = M + 1;
for (I1 = 1; I1 <= D; I1++) {
XN[I] = XN[I] + A[M] * XY[I1];
M = M + 1;
for (I2 = I1; I2 <= D; I2++) {
XN[I] = XN[I] + A[M] * XY[I1] * XY[I2];
M = M + 1;
for (I3 = I2; O > 2 && I3 <= D; I3++) {
XN[I] = XN[I] + A[M] * XY[I1] * XY[I2] * XY[I3];
M = M + 1;
for (I4 = I3; O > 3 && I4 <= D; I4++) {
XN[I] = XN[I] + A[M] * XY[I1] * XY[I2] * XY[I3] * XY[I4];
M = M + 1;
for (I5 = I4; O > 4 && I5 <= D; I5++) {
XN[I] = XN[I] + A[M] * XY[I1] * XY[I2] * XY[I3] * XY[I4] * XY[I5];
M = M + 1;
}}}}}
if (ODE == 1) XN[I] = XY[I] + EPS * XN[I];
}
XNEW = XN[1];
YNEW = XN[2];
ZNEW = XN[3];
WNEW = XN[4];
}
N = N + 1;
M = M - 1;
}
void fun2100() /* Display results */
{
if (N >= 100 && N <= 1000) {
if (X < XMIN) XMIN = X;
if (X > XMAX) XMAX = X;
if (Yy < YMIN) YMIN = Yy;
if (Yy > YMAX) YMAX = Yy;
if (Z < ZMIN) ZMIN = Z;
if (Z > ZMAX) ZMAX = Z;
if (W < WMIN) WMIN = W;
if (W > WMAX) WMAX = W;
}
if ((int)N == 1000) fun3100(); /* Resize the screen */
XS[P] = X;
YS[P] = Yy;
ZS[P] = Z;
WS[P] = W;
P = (P + 1) % 500;
I = (P + 500 - PREV) % 500;
if (D == 1) {
XP = XS[I];
YP = XNEW;
}
else {
XP = X;
YP = Yy;
}
if (N >= 1000 && XP > XL && XP < XH && YP > YL && YP < YH) {
if (PJT == 1) fun4100(); /* Project onto a sphere */
if (PJT == 2) fun6700(); /* Project onto a horizontal cylinder */
if (PJT == 3) fun6800(); /* Project onto a vertical cylinder */
if (PJT == 4) fun6900(); /* Project onto a torus */
fun5000(); /* Plot point on screen */
if (SND == 1) fun3500(); /* Produce sound */
}
}
void fun2400() /* Test results */
{
if (fabs(XNEW) + fabs(YNEW) + fabs(ZNEW) + fabs(WNEW) > 1000000) T = 2;
if (QM != 2) { /* Speed up evaluation mode */
fun2900(); /* Calculate Lyapunov exponent */
fun3900(); /* Calculate fractal dimension */
if (QM == 0) { /* Skip tests unless in search mode */
if (N >= NMAX) { /* Strange attractor found */
T = 2;
fun4900(); /* Save attractor to disk file SA.DIC */
}
if (fabs(XNEW - X) + fabs(YNEW - Yy) + fabs(ZNEW - Z) + fabs(WNEW - W) < .000001) T = 2;
if (N > 100 && L < .005) T = 2; /* Limit cycle */
}
}
if (kbhit()) Q = getch(); else Q = 0;
if (Q) fun3600(); /* Respond to user command */
if (SAV > 0) if (N > 1000 && N < 17001) fun7000(); /* Save data */
X = XNEW; /* Update value of X */
Yy = YNEW;
Z = ZNEW;
W = WNEW;
}
void fun2600() /* Get coefficients */
{
if (QM == 2) { /* In evaluate mode */
fgets(CODE, 515, F1);
if (feof(F1)) {
QM = 0;
fun6000(); /* Update SA.DIC file */
}
else {
fun4700(); /* Get dimension and order */
fun5600(); /* Set colors */
}
}
if (QM == 0) { /* In search mode */
O = 2 + (int)floor((OMAX - 1) * (float)rand() / 32768.0);
CODE[0] = 59 + 4 * D + O + 8 * ODE;
if (ODE > 1) CODE[0] = 87 + ODE;
fun4700(); /* Get value of M */
for (I = 1; I <= M; I++) { /* Construct CODE */
fun2800(); /* Shuffle random numbers */
CODE[I] = 65 + (int)floor(25 * RAN);
}
CODE[M + 1] = 0;
}
for (I = 1; I <= M; I++) { /* Convert CODE to coefficient values */
A[I] = (CODE[I] - 77) / 10.0;
}
}
void fun2800() /* Shuffle random numbers */
{
if (V[0] == 0) for (J = 0; J <= 99; J++) {V[J] = (float)rand() / 32768.0;}
J = (int)floor(100 * RAN);
RAN = V[J];
V[J] = (float)rand() / 32768.0;
}
void fun2900() /* Calculate Lyapunov exponent */
{
XSAVE = XNEW;
YSAVE = YNEW;
ZSAVE = ZNEW;
WSAVE = WNEW;
X = XE;
Yy = YE;
Z = ZE;
W = WE;
N = N - 1;
fun1700(); /* Reiterate equations */
DLX = XNEW - XSAVE;
DLY = YNEW - YSAVE;
DLZ = ZNEW - ZSAVE;
DLW = WNEW - WSAVE;
DL2 = DLX * DLX + DLY * DLY + DLZ * DLZ + DLW * DLW;
if (DL2 > 0) { /* Check for division by zero */
DF = 1E12 * DL2;
RS = 1 / sqrt(DF);
XE = XSAVE + RS * (XNEW - XSAVE);
YE = YSAVE + RS * (YNEW - YSAVE);
ZE = ZSAVE + RS * (ZNEW - ZSAVE);
WE = WSAVE + RS * (WNEW - WSAVE);
XNEW = XSAVE;
YNEW = YSAVE;
ZNEW = ZSAVE;
WNEW = WSAVE;
LSUM = LSUM + log(DF);
NL = NL + 1;
L = .721347 * LSUM / NL;
if (ODE == 1 || ODE == 7) L = L / EPS;
if (N > 1000 && (int)N % 10 == 0) {
gotoxy(WID - 4, 1);
printf("%5.2f", L);
}
}
}
void fun3100() /* Resize the screen */
{
int poly[8];
setcolor(WH);
if (D == 1) {
YMIN = XMIN;
YMAX = XMAX;
}
if (XMAX - XMIN < .000001) {
XMIN = XMIN - .0000005;
XMAX = XMAX + .0000005;
}
if (YMAX - YMIN < .000001) {
YMIN = YMIN - .0000005;
YMAX = YMAX + .0000005;
}
if (ZMAX - ZMIN < .000001) {
ZMIN = ZMIN - .0000005;
ZMAX = ZMAX + .0000005;
}
if (WMAX - WMIN < .000001) {
WMIN = WMIN - .0000005;
WMAX = WMAX + .0000005;
}
MX = .1 * (XMAX - XMIN);
MY = .1 * (YMAX - YMIN);
XL = XMIN - MX;
XH = XMAX + MX;
YL = YMIN - MY;
YH = YMAX + 1.5 * MY;
SW = (getmaxx() + 1) / (XH - XL);
SH = (getmaxy() + 1) / (YL - YH);
X0 = SW * XL;
Y0 = SH * YH;
cleardevice();
YH = YH - .5 * MY;
XA = (XL + XH) / 2;
YA = (YL + YH) / 2;
if (D > 2) {
ZA = (ZMAX + ZMIN) / 2;
if (TRD == 1) {
setfillstyle(1, COLR[1]);
poly[0] = 0;
poly[1] = getmaxy() / 24;
poly[2] = getmaxx();
poly[3] = getmaxy() / 24;
poly[4] = getmaxx();
poly[5] = getmaxy();
poly[6] = 0;
poly[7] = getmaxy();
fillpoly(4, poly);
fun5400(); /* Plot background grid */
}
if (TRD == 4) {
setfillstyle(1, WH);
poly[0] = 0;
poly[1] = getmaxy() / 24;
poly[2] = getmaxx();
poly[3] = getmaxy() / 24;
poly[4] = getmaxx();
poly[5] = getmaxy();
poly[6] = 0;
poly[7] = getmaxy();
fillpoly(4, poly);
}
if (TRD == 5) {
moveto(SW * XA - X0, SH * YL - Y0);
lineto(SW * XA - X0, SH * YH - Y0);
}
if (TRD == 6) {
for (I = 1; I <= 3; I++) {
XP = XL + I * (XH - XL) / 4;
moveto(SW * XP - X0, SH * YL - Y0);
lineto(SW * XP - X0, SH * YH - Y0);
YP = YL + I * (YH - YL) / 4;
moveto(SW * XL - X0, SH * YP - Y0);
lineto(SW * XH - X0, SH * YP - Y0);
}
}
}
if (PJT != 1) rectangle(0, getmaxy() / 24, getmaxx(), getmaxy());
if (PJT == 1 && TRD < 5) circle(SW * XA - X0, SH * YA - Y0, 23 * getmaxy() / 48);
TT = 3.1416 / (XMAX - XMIN);
PT = 3.1416 / (YMAX - YMIN);
if (QM == 2) { /* In evaluate mode */
gotoxy(1, 1);
printf("<Space Bar>: Discard <Enter>: Save");
if (WID >= 80) {
gotoxy(49, 1);
printf("<Esc>: Exit");
gotoxy(69, 1);
printf("%d K left", (int)((filelength(fileno(F1)) - ftell(F1)) / 1024.0));
}}
else {
gotoxy(1, 1);
if (strlen(CODE) < WID - 18)
printf(CODE);
else {
printf("%*.*s...", WID - 23, WID - 23, CODE);
}
gotoxy(WID - 17, 1);
printf("F =");
gotoxy(WID - 7, 1);
printf("L = ");
}
TIA = .05; /* Tangent of illumination angle */
XZ = -TIA * (XMAX - XMIN) / (ZMAX - ZMIN);
YZ = TIA * (YMAX - YMIN) / (ZMAX - ZMIN);
}
void fun3500() /* Produce sound */
{
FREQ = 220 * pow(2, (int)(36 * (XNEW - XL) / (XH - XL)) / 12.0);
DUR = 1;
if (D > 1) DUR = pow(2, (int)floor(.5 * (YH - YL) / (YNEW - 9 * YL / 8 + YH / 8)));
sound(FREQ);
delay(55 * DUR);
nosound();
}
void fun3600() /* Respond to user command */
{
if (Q > 96) Q = Q - 32; /* Convert to upper case */
if (QM == 2) fun5800(); /* Process evaluation command */
if (strchr("ACDEHINPRSVX", Q) == 0) fun4200(); /* Display menu screen */
if (Q == 'A') {
T = 1;
QM = 0;
}
if (ODE > 1) D = ODE + 5;
if (ODE == 1) D = D + 2;
if (Q == 'C') if (N > 999) N = 999;
if (Q == 'D') {
D = 1 + D % 12;
T = 1;
}
if (D > 6) {
ODE = D - 5;
D = 4;
}
else {
if (D > 4) {
ODE = 1;
D = D - 2;
}
else ODE = 0;
}
if (Q == 'E') {
T = 1;
QM = 2;
}
if (Q == 'H') {
FTH = (FTH + 1) % 3;
T = 3;
if (N > 999) {
N = 999;
fun5600(); /* Set colors */
}
}
if (Q == 'I') {
if (T != 1) {
closegraph();
textcolor(15);
textbackground(1L);
_setcursortype(_NORMALCURSOR);
clrscr();
printf("Code? ");
I = 0;
CODE[0] = 0;
do {
CODE[I] = getche();
if (CODE[I] == 8 && I >= 0) I = I - 2;
if (CODE[I] == 27) {
I = 0;
CODE[I] = 13;
_setcursortype(_NOCURSOR);
}
}
while (CODE[I++] != 13 && I < 506);
CODE[I - 1] = 0;
if (CODE[0] == 0) {
Q = ' ';
clrscr();}
else {
T = 1;
QM = 1;
fun4700();
}
}
}
if (Q == 'N') {
NMAX = 10 * (NMAX - 1000) + 1000;
if (NMAX > 1E10) NMAX = 2000;
}
if (Q == 'P') {
PJT = (PJT + 1) % 5;
T = 3;
if (N > 999) N = 999;
}
if (Q == 'R') {
TRD = (TRD + 1) % 7;
T = 3;
if (N > 999) {
N = 999;
fun5600(); /* Get dimension and order */
}
}
if (Q == 'S') {
SND = (SND + 1) % 2;
T = 3;
}
if (Q == 'V') {
SAV = (SAV + 1) % 5;
FAV[0] = 87 + SAV % 4;
T = 3;
if (N > 999) N = 999;
}
if (Q == 'X') T = 0;
}
void fun3900() /* Calculate fractal dimension */
{
if (N >= 1000) { /* Wait for transient to settle */
if ((int)N == 1000) {
D2MAX = pow(XMAX - XMIN, 2);
D2MAX = D2MAX + pow(YMAX - YMIN, 2);
D2MAX = D2MAX + pow(ZMAX - ZMIN, 2);
D2MAX = D2MAX + pow(WMAX - WMIN, 2);
}
J = (P + 1 + (int)floor(480 * (float)rand() / 32768.0)) % 500;
DX = XNEW - XS[J];
DY = YNEW - YS[J];
DZ = ZNEW - ZS[J];
DW = WNEW - WS[J];
D2 = DX * DX + DY * DY + DZ * DZ + DW * DW;
if (D2 < .001 * TWOD * D2MAX) N2 = N2 + 1;
if (D2 <= .00001 * TWOD * D2MAX) {
N1 = N1 + 1;
F = .434294 * log(N2 / (N1 - .5));
gotoxy(WID - 14, 1);
printf("%5.2f", F);
}
}
}
void fun4100() /* Project onto a sphere */
{
TH = TT * (XMAX - XP);
PH = PT * (YMAX - YP);
XP = XA + .36 * (XH - XL) * cos(TH) * sin(PH);
YP = YA + .5 * (YH - YL) * cos(PH);
}
void fun4200() /* Display menu screen */
{
closegraph();
textcolor(15);
textbackground(1L);
_setcursortype(_NOCURSOR); /* Turn cursor off */
clrscr();
while (Q == 0 || strchr("AEIX", Q) == 0) {
gotoxy(27, 1);
printf("STRANGE ATTRACTOR PROGRAM\n");
printf("%26cIBM PC TurboC Version 2.0\n", ' ');
printf("%26c(c) 1993 by J. C. Sprott\n", ' ');
printf("\n");
printf("\n");
printf("%26cA: Search for attractors\n", ' ');
printf("%26cC: Clear screen and restart\n", ' ');
if (ODE > 1) {
printf("%26cD: System is 4-D special map %c \n", ' ', 87 + ODE);}
else {
printf("%26cD: System is %d-D polynomial ", ' ', D);
if (ODE == 1) printf("ODE\n"); else printf("map\n");
}
printf("%26cE: Evaluate attractors\n", ' ');
printf("%26cH: Fourth dimension is ", ' ');
if (FTH == 0) printf("projection\n");
if (FTH == 1) printf("bands \n");
if (FTH == 2) printf("colors \n");
printf("%26cI: Input code from keyboard\n", ' ');
printf("%26cN: Number of iterations is 10^%1.0f\n", ' ', log10(NMAX - 1000));
printf("%26cP: Projection is ", ' ');
if (PJT == 0) printf("planar \n");
if (PJT == 1) printf("spherical\n");
if (PJT == 2) printf("horiz cyl\n");
if (PJT == 3) printf("vert cyl \n");
if (PJT == 4) printf("toroidal \n");
printf("%26cR: Third dimension is ", ' ');
if (TRD == 0) printf("projection\n");
if (TRD == 1) printf("shadow \n");
if (TRD == 2) printf("bands \n");
if (TRD == 3) printf("colors \n");
if (TRD == 4) printf("anaglyph \n");
if (TRD == 5) printf("stereogram\n");
if (TRD == 6) printf("slices \n");
printf("%26cS: Sound is ", ' ');
if (SND == 0) printf("off\n");
if (SND == 1) printf("on \n");
printf("%26cV: ", ' ');
if (SAV == 0) printf("No data will be saved \n");
if (SAV > 0) printf("%c will be saved in %cDATA.DAT\n", FAV[0], FAV[0]);
printf("%26cX: Exit program", ' ');
if (kbhit()) Q = getch(); else Q = 0;
if (Q) fun3600(); /* Respond to user command */
}
}
void fun4700() /* Get dimension and order */
{
D = 1 + (int)floor((CODE[0] - 65) / 4);
if (D > 6) {
ODE = CODE[0] - 87;
D = 4;
fun6200(); /* Special function */
}
else {
if (D > 4) {
D = D - 2;
ODE = 1;
}
else ODE = 0;
O = 2 + (CODE[0] - 65) % 4;
M = 1;
for (I = 1; I <= D; I++) {M = M * (O + I);}
if (D > 2) for (I = 3; I <= D; I++) {M = M / (I - 1);}
}
if (strlen(CODE) != M + 1 && QM == 1) {
printf("\a"); /* Illegal code warning */
while (strlen(CODE) < M + 1) strcat(CODE, "M");
if (strlen(CODE) > M + 1) CODE[M + 1] = 0;
}
}
void fun4900() /* Save attractor to disk file SA.DIC */
{
F1 = fopen("SA.DIC", "a");
fprintf(F1, "%s%5.2f%5.2f\n", CODE, F, L);
fclose(F1);
}
void fun5000() /* Plot point on screen */
{
C4 = WH;
if (D > 3) {
if (FTH == 1) if ((int)floor(30 * (W - WMIN) / (WMAX - WMIN)) % 2) return;
if (FTH == 2) C4 = 1 + (int)floor(NC * (W - WMIN) / (WMAX - WMIN) + NC) % NC;
}
if (D < 3) { /* Skip 3-D stuff */
putpixel(SW * XP - X0, SH * YP - Y0, 15);
return;
}
if (TRD == 0) {
putpixel(SW * XP - X0, SH * YP - Y0, C4);
}
if (TRD == 1) {
if (D > 3 && FTH == 2) {
putpixel(SW * XP - X0, SH * YP - Y0, C4);
}
else {
C = getpixel(SW * XP - X0, SH * YP - Y0);
if (C == COLR[2]) {
putpixel(SW * XP - X0, SH * YP - Y0, COLR[3]);}
else {
if (C != COLR[3]) {
putpixel(SW * XP - X0, SH * YP - Y0, COLR[2]);
}
}
}
XP = XP - XZ * (Z - ZMIN);
YP = YP - YZ * (Z - ZMIN);
if (getpixel(SW * XP - X0, SH * YP - Y0) == COLR[1]) {
putpixel(SW * XP - X0, SH * YP - Y0, 0);
}
}
if (TRD == 2) {
if (D > 3 && FTH == 2 && ((int)floor(15 * (Z - ZMIN) / (ZMAX - ZMIN) + 2) % 2) == 1) {
C = C4;}
else {
C = COLR[(int)floor(60 * (Z - ZMIN) / (ZMAX - ZMIN) + 4) % 4];
}
putpixel(SW * XP - X0, SH * YP - Y0, C);
}
if (TRD == 3) {
putpixel(SW * XP - X0, SH * YP - Y0, COLR[(int)floor(NC * (Z - ZMIN) / (ZMAX - ZMIN) + NC) % NC]);
}
if (TRD == 4) {
XRT = XP + XZ * (Z - ZA);
C = getpixel(SW * XRT - X0, SH * YP - Y0);
if (C == WH) {
putpixel(SW * XRT - X0, SH * YP - Y0, RD);
}
if (C == BY) {
putpixel(SW * XRT - X0, SH * YP - Y0, BK);
}
XLT = XP - XZ * (Z - ZA);
C = getpixel(SW * XLT - X0, SH * YP - Y0);
if (C == WH) {
putpixel(SW * XLT - X0, SH * YP - Y0, BY);
}
if (C == RD) {
putpixel(SW * XLT - X0, SH * YP - Y0, BK);
}
}
if (TRD == 5) {
HSF = 2; /* Horizontal scale factor */
XRT = XA + (XP + XZ * (Z - ZA) - XL) / HSF;
putpixel(SW * XRT - X0, SH * YP - Y0, C4);
XLT = XA + (XP - XZ * (Z - ZA) - XH) / HSF;
putpixel(SW * XLT - X0, SH * YP - Y0, C4);
}
if (TRD == 6) {
DZ = (15 * (Z - ZMIN) / (ZMAX - ZMIN) + .5) / 16;
XP = (XP - XL + ((int)floor(16 * DZ) % 4) * (XH - XL)) / 4 + XL;
YP = (YP - YL + (3 - (int)floor(4 * DZ) % 4) * (YH - YL)) / 4 + YL;
putpixel(SW * XP - X0, SH * YP - Y0, C4);
}
}
void fun5400() /* Plot background grid */
{
setcolor(0);
for (I = 0; I <= 15; I++) { /* Draw 15 vertical grid lines */
XP = XMIN + I * (XMAX - XMIN) / 15;
moveto(SW * XP - X0, SH * YMIN - Y0);
lineto(SW * XP - X0, SH * YMAX - Y0);
}
for (I = 0; I <= 10; I++) { /* Draw 10 horizontal grid lines */
YP = YMIN + I * (YMAX - YMIN) / 10;
moveto(SW * XMIN - X0, SH * YP - Y0);
lineto(SW * XMAX - X0, SH * YP - Y0);
}
}
void fun5600() /* Set colors */
{
NC = 15; /* Number of colors */
COLR[0] = 0;
COLR[1] = 8;
COLR[2] = 7;
COLR[3] = 15;
if (TRD == 3 || (D > 3 && FTH == 2 && TRD != 1)) {
for (I = 0; I <= NC; I++) COLR[I] = I + 1;
}
WH = 15;
BK = 8;
RD = 12;
BY = 11;
}
void fun5800() /* Process evaluation command */
{
if (Q == ' ') {
T = 2;
NE = NE + 1;
cleardevice();
}
if (Q == 13) {
T = 2;
NE = NE + 1;
cleardevice();
fun5900(); /* Save favorite attractors to disk */
}
if (Q == 27) {
cleardevice();
fun6000(); /* Update SA.DIC file */
Q = ' ';
QM = 0;
}
else {
if (strchr("CHNPRVS", Q) == 0) Q = 0;
}
}
void fun5900() /* Save favorite attractors to disk file FAVORITE.DIC */
{
F2 = fopen("FAVORITE.DIC", "a");
fprintf(F2, CODE);
fclose(F2);
}
void fun6000() /* Update SA.DIC file */
{
gotoxy(9, 11);
printf("Evaluation complete\n");
gotoxy(8, 12);
printf(" %d cases evaluated", (int)NE);
F2 = fopen("SATEMP.DIC", "w");
if (QM == 2) fprintf(F2, CODE);
while (feof(F1) == 0) {
fgets(CODE, 515, F1);
if (feof(F1) == 0) fprintf(F2, CODE);
}
fcloseall();
remove("SA.DIC");
rename("SATEMP.DIC", "SA.DIC");
}
void fun6200() /* Special function definitions */
{
ZNEW = X * X + Yy * Yy; /* Default 3rd and 4th dimension */
WNEW = (N - 100) / 900;
if (N > 1000) WNEW = (N - 1000) / (NMAX - 1000);
if (ODE == 2) {
M = 10;
XNEW = A[1] + A[2] * X + A[3] * Yy + A[4] * fabs(X) + A[5] * fabs(Yy);
YNEW = A[6] + A[7] * X + A[8] * Yy + A[9] * fabs(X) + A[10] * fabs(Yy);
}
if (ODE == 3) {
M = 14;
XNEW = A[1] + A[2] * X + A[3] * Yy + ((int)(A[4] * X + .5) & (int)(A[5] * Yy + .5)) + ((int)(A[6] * X + .5) | (int)(A[7] * Yy + .5));
YNEW = A[8] + A[9] * X + A[10] * Yy + ((int)(A[11] * X + .5) & (int)(A[12] * Yy + .5)) + ((int)(A[13] * X + .5) | (int)(A[14] * Yy + .5));
}
if (ODE == 4) {
M = 14;
XNEW = A[1] + A[2] * X + A[3] * Yy + A[4] * pow(fabs(X), A[5]) + A[6] * pow(fabs(Yy), A[7]);
YNEW = A[8] + A[9] * X + A[10] * Yy + A[11] * pow(fabs(X), A[12]) + A[13] * pow(fabs(Yy), A[14]);
}
if (ODE == 5) {
M = 18;
XNEW = A[1] + A[2] * X + A[3] * Yy + A[4] * sin(A[5] * X + A[6]) + A[7] * sin(A[8] * Yy + A[9]);
YNEW = A[10] + A[11] * X + A[12] * Yy + A[13] * sin(A[14] * X + A[15]) + A[16] * sin(A[17] * Yy + A[18]);
}
if (ODE == 6) {
M = 6;
if (N < 2) {
AL = TWOPI / (13 + 10 * A[6]);
SINAL = sin(AL);
COSAL = cos(AL);
}
DUM = X + A[2] * sin(A[3] * Yy + A[4]);
XNEW = 10 * A[1] + DUM * COSAL + Yy * SINAL;
YNEW = 10 * A[5] - DUM * SINAL + Yy * COSAL;
}
if (ODE == 7) {
M = 9;
XNEW = X + EPS * A[1] * Yy;
YNEW = Yy + EPS * (A[2] * X + A[3] * X * X * X + A[4] * X * X * Yy + A[5] * X * Yy * Yy + A[6] * Yy + A[7] * Yy * Yy * Yy + A[8] * sin(Z));
ZNEW = Z + EPS * (A[9] + 1.3);
if (ZNEW > TWOPI) ZNEW = ZNEW - TWOPI;
}
}
void fun6700() /* Project onto a horizontal cylinder */
{
PH = PT * (YMAX - YP);
YP = YA + .5 * (YH - YL) * cos(PH);
}
void fun6800() /* Project onto a vertical cylinder */
{
TH = TT * (XMAX - XP);
XP = XA + .5 * (XH - XL) * cos(TH);
}
void fun6900() /* Project onto a torus (unity aspect ratio) */
{
TH = TT * (XMAX - XP);
PH = 2 * PT * (YMAX - YP);
XP = XA + .18 * (XH - XL) * (1 + cos(TH)) * sin(PH);
YP = YA + .25 * (YH - YL) * (1 + cos(TH)) * cos(PH);
}
void fun7000() /* Save data */
{
if ((int)N == 1001) {
fclose(F3);
F3 = fopen(FAV, "w");
}
if (SAV == 1) DUM = XNEW;
if (SAV == 2) DUM = YNEW;
if (SAV == 3) DUM = ZNEW;
if (SAV == 4) DUM = WNEW;
fprintf(F3, "%f\n", DUM);
}
