« back to all changes in this revision
Viewing changes to miscfrac.c
- browse files at revision 4
- compare with another revision
- download diff
- download tarball
- view history from revision 4
- Committer: Bazaar Package Importer
- Author(s): Riku Voipio
- Date: 2010-11-24 21:24:54 UTC
- mfrom: (1.2.2 upstream)
- Revision ID: james.westby@ubuntu.com-20101124212454-a2vxdh9pk5p5lhlm
Tags: 20.4.10-1
* New upstream version
* enable autobuilding, Closes: #410674, #587959
* update watch file. Closes: #449889
* remove bashism from rules, Closes: #581474
* enable autobuilding, Closes: #410674, #587959
* update watch file. Closes: #449889
* remove bashism from rules, Closes: #581474
added
removed
miscfrac.c
1
/*
2
3
Miscellaneous fractal-specific code (formerly in CALCFRAC.C)
4
5
*/
6
7
#include <string.h>
8
#include <limits.h>
9
/* see Fractint.c for a description of the "include" hierarchy */
10
#include "port.h"
11
#include "prototyp.h"
12
#include "fractype.h"
13
#include "targa_lc.h"
14
15
/* routines in this module */
16
17
static void set_Plasma_palette(void);
18
static U16 _fastcall adjust(int xa,int ya,int x,int y,int xb,int yb);
19
static void _fastcall subDivide(int x1,int y1,int x2,int y2);
20
static int _fastcall new_subD (int x1,int y1,int x2,int y2, int recur);
21
static void verhulst(void);
22
static void Bif_Period_Init(void);
23
static int _fastcall Bif_Periodic(long);
24
static void set_Cellular_palette(void);
25
26
U16 (_fastcall *getpix)(int,int) = (U16(_fastcall *)(int,int))getcolor;
27
28
typedef void (_fastcall *PLOT)(int,int,int);
29
30
/***************** standalone engine for "test" ********************/
31
32
int test(void)
33
{
34
int startrow,startpass,numpasses;
35
startrow = startpass = 0;
36
if (resuming)
37
{
38
start_resume();
39
get_resume(sizeof(startrow),&startrow,sizeof(startpass),&startpass,0);
40
end_resume();
41
}
42
if(teststart()) /* assume it was stand-alone, doesn't want passes logic */
43
return(0);
44
numpasses = (stdcalcmode == '1') ? 0 : 1;
45
for (passes=startpass; passes <= numpasses ; passes++)
46
{
47
for (row = startrow; row <= iystop; row=row+1+numpasses)
48
{
49
for (col = 0; col <= ixstop; col++) /* look at each point on screen */
50
{
51
register int color;
52
init.x = dxpixel();
53
init.y = dypixel();
54
if(keypressed())
55
{
56
testend();
57
alloc_resume(20,1);
58
put_resume(sizeof(row),&row,sizeof(passes),&passes,0);
59
return(-1);
60
}
61
color = testpt(init.x,init.y,parm.x,parm.y,maxit,inside);
62
if (color >= colors) { /* avoid trouble if color is 0 */
63
if (colors < 16)
64
color &= andcolor;
65
else
66
color = ((color-1) % andcolor) + 1; /* skip color zero */
67
}
68
(*plot)(col,row,color);
69
if(numpasses && (passes == 0))
70
(*plot)(col,row+1,color);
71
}
72
}
73
startrow = passes + 1;
74
}
75
testend();
76
return(0);
77
}
78
79
/***************** standalone engine for "plasma" ********************/
80
81
static int iparmx; /* iparmx = parm.x * 8 */
82
static int shiftvalue; /* shift based on #colors */
83
static int recur1=1;
84
static int pcolors;
85
static int recur_level = 0;
86
U16 max_plasma;
87
88
/* returns a random 16 bit value that is never 0 */
89
U16 rand16(void)
90
{
91
U16 value;
92
value = (U16)rand15();
93
value <<= 1;
94
value = (U16)(value + (rand15()&1));
95
if(value < 1)
96
value = 1;
97
return(value);
98
}
99
100
void _fastcall putpot(int x, int y, U16 color)
101
{
102
if(color < 1)
103
color = 1;
104
putcolor(x, y, color >> 8 ? color >> 8 : 1); /* don't write 0 */
105
/* we don't write this if dotmode==11 because the above putcolor
106
was already a "writedisk" in that case */
107
if (dotmode != 11)
108
writedisk(x+sxoffs,y+syoffs,color >> 8); /* upper 8 bits */
109
writedisk(x+sxoffs,y+sydots+syoffs,color&255); /* lower 8 bits */
110
}
111
112
/* fixes border */
113
void _fastcall putpotborder(int x, int y, U16 color)
114
{
115
if((x==0) || (y==0) || (x==xdots-1) || (y==ydots-1))
116
color = (U16)outside;
117
putpot(x,y,color);
118
}
119
120
/* fixes border */
121
void _fastcall putcolorborder(int x, int y, int color)
122
{
123
if((x==0) || (y==0) || (x==xdots-1) || (y==ydots-1))
124
color = outside;
125
if(color < 1)
126
color = 1;
127
putcolor(x,y,color);
128
}
129
130
U16 _fastcall getpot(int x, int y)
131
{
132
U16 color;
133
134
color = (U16)readdisk(x+sxoffs,y+syoffs);
135
color = (U16)((color << 8) + (U16) readdisk(x+sxoffs,y+sydots+syoffs));
136
return(color);
137
}
138
139
static int plasma_check; /* to limit kbd checking */
140
141
static U16 _fastcall adjust(int xa,int ya,int x,int y,int xb,int yb)
142
{
143
S32 pseudorandom;
144
pseudorandom = ((S32)iparmx)*((rand15()-16383));
145
/* pseudorandom = pseudorandom*(abs(xa-xb)+abs(ya-yb));*/
146
pseudorandom = pseudorandom * recur1;
147
pseudorandom = pseudorandom >> shiftvalue;
148
pseudorandom = (((S32)getpix(xa,ya)+(S32)getpix(xb,yb)+1)>>1)+pseudorandom;
149
if(max_plasma == 0)
150
{
151
if (pseudorandom >= pcolors)
152
pseudorandom = pcolors-1;
153
}
154
else if (pseudorandom >= (S32)max_plasma)
155
pseudorandom = max_plasma;
156
if(pseudorandom < 1)
157
pseudorandom = 1;
158
plot(x,y,(U16)pseudorandom);
159
return((U16)pseudorandom);
160
}
161
162
163
static int _fastcall new_subD (int x1,int y1,int x2,int y2, int recur)
164
{
165
int x,y;
166
int nx1;
167
int nx;
168
int ny1, ny;
169
S32 i, v;
170
171
struct sub {
172
BYTE t; /* top of stack */
173
int v[16]; /* subdivided value */
174
BYTE r[16]; /* recursion level */
175
};
176
177
static struct sub subx, suby;
178
179
/*
180
recur1=1;
181
for (i=1;i<=recur;i++)
182
recur1 = recur1 * 2;
183
recur1=320/recur1;
184
*/
185
recur1 = (int)(320L >> recur);
186
suby.t = 2;
187
ny = suby.v[0] = y2;
188
ny1 = suby.v[2] = y1;
189
suby.r[0] = suby.r[2] = 0;
190
suby.r[1] = 1;
191
y = suby.v[1] = (ny1 + ny) >> 1;
192
193
while (suby.t >= 1)
194
{
195
if ((++plasma_check & 0x0f) == 1)
196
if(keypressed())
197
{
198
plasma_check--;
199
return(1);
200
}
201
while (suby.r[suby.t-1] < (BYTE)recur)
202
{
203
/* 1. Create new entry at top of the stack */
204
/* 2. Copy old top value to new top value. */
205
/* This is largest y value. */
206
/* 3. Smallest y is now old mid point */
207
/* 4. Set new mid point recursion level */
208
/* 5. New mid point value is average */
209
/* of largest and smallest */
210
211
suby.t++;
212
ny1 = suby.v[suby.t] = suby.v[suby.t-1];
213
ny = suby.v[suby.t-2];
214
suby.r[suby.t] = suby.r[suby.t-1];
215
y = suby.v[suby.t-1] = (ny1 + ny) >> 1;
216
suby.r[suby.t-1] = (BYTE)(max(suby.r[suby.t], suby.r[suby.t-2])+1);
217
}
218
subx.t = 2;
219
nx = subx.v[0] = x2;
220
nx1 = subx.v[2] = x1;
221
subx.r[0] = subx.r[2] = 0;
222
subx.r[1] = 1;
223
x = subx.v[1] = (nx1 + nx) >> 1;
224
225
while (subx.t >= 1)
226
{
227
while (subx.r[subx.t-1] < (BYTE)recur)
228
{
229
subx.t++; /* move the top ofthe stack up 1 */
230
nx1 = subx.v[subx.t] = subx.v[subx.t-1];
231
nx = subx.v[subx.t-2];
232
subx.r[subx.t] = subx.r[subx.t-1];
233
x = subx.v[subx.t-1] = (nx1 + nx) >> 1;
234
subx.r[subx.t-1] = (BYTE)(max(subx.r[subx.t],
235
subx.r[subx.t-2])+1);
236
}
237
238
if ((i = getpix(nx, y)) == 0)
239
i = adjust(nx,ny1,nx,y ,nx,ny);
240
v = i;
241
if ((i = getpix(x, ny)) == 0)
242
i = adjust(nx1,ny,x ,ny,nx,ny);
243
v += i;
244
if(getpix(x,y) == 0)
245
{
246
if ((i = getpix(x, ny1)) == 0)
247
i = adjust(nx1,ny1,x ,ny1,nx,ny1);
248
v += i;
249
if ((i = getpix(nx1, y)) == 0)
250
i = adjust(nx1,ny1,nx1,y ,nx1,ny);
251
v += i;
252
plot(x,y,(U16)((v + 2) >> 2));
253
}
254
255
if (subx.r[subx.t-1] == (BYTE)recur) subx.t = (BYTE)(subx.t - 2);
256
}
257
258
if (suby.r[suby.t-1] == (BYTE)recur) suby.t = (BYTE)(suby.t - 2);
259
}
260
return(0);
261
}
262
263
static void _fastcall subDivide(int x1,int y1,int x2,int y2)
264
{
265
int x,y;
266
S32 v,i;
267
if ((++plasma_check & 0x7f) == 1)
268
if(keypressed())
269
{
270
plasma_check--;
271
return;
272
}
273
if(x2-x1<2 && y2-y1<2)
274
return;
275
recur_level++;
276
recur1 = (int)(320L >> recur_level);
277
278
x = (x1+x2)>>1;
279
y = (y1+y2)>>1;
280
if((v=getpix(x,y1)) == 0)
281
v=adjust(x1,y1,x ,y1,x2,y1);
282
i=v;
283
if((v=getpix(x2,y)) == 0)
284
v=adjust(x2,y1,x2,y ,x2,y2);
285
i+=v;
286
if((v=getpix(x,y2)) == 0)
287
v=adjust(x1,y2,x ,y2,x2,y2);
288
i+=v;
289
if((v=getpix(x1,y)) == 0)
290
v=adjust(x1,y1,x1,y ,x1,y2);
291
i+=v;
292
293
if(getpix(x,y) == 0)
294
plot(x,y,(U16)((i+2)>>2));
295
296
subDivide(x1,y1,x ,y);
297
subDivide(x ,y1,x2,y);
298
subDivide(x ,y ,x2,y2);
299
subDivide(x1,y ,x ,y2);
300
recur_level--;
301
}
302
303
304
int plasma()
305
{
306
int i,k, n;
307
U16 rnd[4];
308
int OldPotFlag, OldPot16bit;
309
310
OldPotFlag=OldPot16bit=plasma_check = 0;
311
312
if(colors < 4) {
313
static FCODE plasmamsg[]={
314
"\
315
Plasma Clouds can currently only be run in a 4-or-more-color video\n\
316
mode (and color-cycled only on VGA adapters [or EGA adapters in their\n\
317
640x350x16 mode])." };
318
stopmsg(0,plasmamsg);
319
return(-1);
320
}
321
iparmx = (int)(param[0] * 8);
322
if (parm.x <= 0.0) iparmx = 0;
323
if (parm.x >= 100) iparmx = 800;
324
param[0] = (double)iparmx / 8.0; /* let user know what was used */
325
if (param[1] < 0) param[1] = 0; /* limit parameter values */
326
if (param[1] > 1) param[1] = 1;
327
if (param[2] < 0) param[2] = 0; /* limit parameter values */
328
if (param[2] > 1) param[2] = 1;
329
if (param[3] < 0) param[3] = 0; /* limit parameter values */
330
if (param[3] > 1) param[3] = 1;
331
332
if ((!rflag) && param[2] == 1)
333
--rseed;
334
if (param[2] != 0 && param[2] != 1)
335
rseed = (int)param[2];
336
max_plasma = (U16)param[3]; /* max_plasma is used as a flag for potential */
337
338
if(max_plasma != 0)
339
{
340
if (pot_startdisk() >= 0)
341
{
342
/* max_plasma = (U16)(1L << 16) -1; */
343
max_plasma = 0xFFFF;
344
if(outside >= 0)
345
plot = (PLOT)putpotborder;
346
else
347
plot = (PLOT)putpot;
348
getpix = getpot;
349
OldPotFlag = potflag;
350
OldPot16bit = pot16bit;
351
}
352
else
353
{
354
max_plasma = 0; /* can't do potential (startdisk failed) */
355
param[3] = 0;
356
if(outside >= 0)
357
plot = putcolorborder;
358
else
359
plot = putcolor;
360
getpix = (U16(_fastcall *)(int,int))getcolor;
361
}
362
}
363
else
364
{
365
if(outside >= 0)
366
plot = putcolorborder;
367
else
368
plot = putcolor;
369
getpix = (U16(_fastcall *)(int,int))getcolor;
370
}
371
srand(rseed);
372
if (!rflag) ++rseed;
373
374
if (colors == 256) /* set the (256-color) palette */
375
set_Plasma_palette(); /* skip this if < 256 colors */
376
377
if (colors > 16)
378
shiftvalue = 18;
379
else
380
{
381
if (colors > 4)
382
shiftvalue = 22;
383
else
384
{
385
if (colors > 2)
386
shiftvalue = 24;
387
else
388
shiftvalue = 25;
389
}
390
}
391
if(max_plasma != 0)
392
shiftvalue = 10;
393
394
if(max_plasma == 0)
395
{
396
pcolors = min(colors, max_colors);
397
for(n = 0; n < 4; n++)
398
rnd[n] = (U16)(1+(((rand15()/pcolors)*(pcolors-1))>>(shiftvalue-11)));
399
}
400
else
401
for(n = 0; n < 4; n++)
402
rnd[n] = rand16();
403
if(debugflag==3600)
404
for(n = 0; n < 4; n++)
405
rnd[n] = 1;
406
407
plot( 0, 0, rnd[0]);
408
plot(xdots-1, 0, rnd[1]);
409
plot(xdots-1,ydots-1, rnd[2]);
410
plot( 0,ydots-1, rnd[3]);
411
412
recur_level = 0;
413
if (param[1] == 0)
414
subDivide(0,0,xdots-1,ydots-1);
415
else
416
{
417
recur1 = i = k = 1;
418
while(new_subD(0,0,xdots-1,ydots-1,i)==0)
419
{
420
k = k * 2;
421
if (k >(int)max(xdots-1,ydots-1))
422
break;
423
if (keypressed())
424
{
425
n = 1;
426
goto done;
427
}
428
i++;
429
}
430
}
431
if (! keypressed())
432
n = 0;
433
else
434
n = 1;
435
done:
436
if(max_plasma != 0)
437
{
438
potflag = OldPotFlag;
439
pot16bit = OldPot16bit;
440
}
441
plot = putcolor;
442
getpix = (U16(_fastcall *)(int,int))getcolor;
443
return(n);
444
}
445
446
#define dac ((Palettetype *)dacbox)
447
static void set_Plasma_palette()
448
{
449
static Palettetype Red = { 63, 0, 0 };
450
static Palettetype Green = { 0, 63, 0 };
451
static Palettetype Blue = { 0, 0,63 };
452
int i;
453
454
if (mapdacbox || colorpreloaded) return; /* map= specified */
455
456
dac[0].red = 0 ;
457
dac[0].green= 0 ;
458
dac[0].blue = 0 ;
459
for(i=1;i<=85;i++)
460
{
461
#ifdef __SVR4
462
dac[i].red = (BYTE)((i*(int)Green.red + (86-i)*(int)Blue.red)/85);
463
dac[i].green = (BYTE)((i*(int)Green.green + (86-i)*(int)Blue.green)/85);
464
dac[i].blue = (BYTE)((i*(int)Green.blue + (86-i)*(int)Blue.blue)/85);
465
466
dac[i+85].red = (BYTE)((i*(int)Red.red + (86-i)*(int)Green.red)/85);
467
dac[i+85].green = (BYTE)((i*(int)Red.green + (86-i)*(int)Green.green)/85);
468
dac[i+85].blue = (BYTE)((i*(int)Red.blue + (86-i)*(int)Green.blue)/85);
469
470
dac[i+170].red = (BYTE)((i*(int)Blue.red + (86-i)*(int)Red.red)/85);
471
dac[i+170].green = (BYTE)((i*(int)Blue.green + (86-i)*(int)Red.green)/85);
472
dac[i+170].blue = (BYTE)((i*(int)Blue.blue + (86-i)*(int)Red.blue)/85);
473
#else
474
dac[i].red = (BYTE)((i*Green.red + (86-i)*Blue.red)/85);
475
dac[i].green = (BYTE)((i*Green.green + (86-i)*Blue.green)/85);
476
dac[i].blue = (BYTE)((i*Green.blue + (86-i)*Blue.blue)/85);
477
478
dac[i+85].red = (BYTE)((i*Red.red + (86-i)*Green.red)/85);
479
dac[i+85].green = (BYTE)((i*Red.green + (86-i)*Green.green)/85);
480
dac[i+85].blue = (BYTE)((i*Red.blue + (86-i)*Green.blue)/85);
481
dac[i+170].red = (BYTE)((i*Blue.red + (86-i)*Red.red)/85);
482
dac[i+170].green = (BYTE)((i*Blue.green + (86-i)*Red.green)/85);
483
dac[i+170].blue = (BYTE)((i*Blue.blue + (86-i)*Red.blue)/85);
484
#endif
485
}
486
SetTgaColors(); /* TARGA 3 June 89 j mclain */
487
spindac(0,1);
488
}
489
490
/***************** standalone engine for "diffusion" ********************/
491
492
#define RANDOM(x) (rand()%(x))
493
494
int diffusion()
495
{
496
int xmax,ymax,xmin,ymin; /* Current maximum coordinates */
497
int border; /* Distance between release point and fractal */
498
int mode; /* Determines diffusion type: 0 = central (classic) */
499
/* 1 = falling particles */
500
/* 2 = square cavity */
501
int colorshift; /* If zero, select colors at random, otherwise shift */
502
/* the color every colorshift points */
503
504
int colorcount,currentcolor;
505
506
int i;
507
double cosine,sine,angle;
508
int x,y;
509
float r, radius;
510
511
if (diskvideo)
512
notdiskmsg();
513
514
x = y = -1;
515
bitshift = 16;
516
fudge = 1L << 16;
517
518
border = (int)param[0];
519
mode = (int)param[1];
520
colorshift = (int)param[2];
521
522
colorcount = colorshift; /* Counts down from colorshift */
523
currentcolor = 1; /* Start at color 1 (color 0 is probably invisible)*/
524
525
if (mode > 2)
526
mode=0;
527
528
if (border <= 0)
529
border = 10;
530
531
srand(rseed);
532
if (!rflag) ++rseed;
533
534
if (mode == 0) {
535
xmax = xdots / 2 + border; /* Initial box */
536
xmin = xdots / 2 - border;
537
ymax = ydots / 2 + border;
538
ymin = ydots / 2 - border;
539
}
540
if (mode == 1) {
541
xmax = xdots / 2 + border; /* Initial box */
542
xmin = xdots / 2 - border;
543
ymin = ydots - border;
544
}
545
if (mode == 2) {
546
if (xdots > ydots)
547
radius = ydots - border;
548
else
549
radius = xdots - border;
550
}
551
if (resuming) /* restore worklist, if we can't the above will stay in place */
552
{
553
start_resume();
554
if (mode != 2)
555
get_resume(sizeof(xmax),&xmax,sizeof(xmin),&xmin,sizeof(ymax),&ymax,
556
sizeof(ymin),&ymin,0);
557
else
558
get_resume(sizeof(xmax),&xmax,sizeof(xmin),&xmin,sizeof(ymax),&ymax,
559
sizeof(radius),&radius,0);
560
end_resume();
561
}
562
563
switch (mode) {
564
case 0: /* Single seed point in the center */
565
putcolor(xdots / 2, ydots / 2,currentcolor);
566
break;
567
case 1: /* Line along the bottom */
568
for (i=0;i<=xdots;i++)
569
putcolor(i,ydots-1,currentcolor);
570
break;
571
case 2: /* Large square that fills the screen */
572
if (xdots > ydots)
573
for (i=0;i<ydots;i++){
574
putcolor(xdots/2-ydots/2 , i , currentcolor);
575
putcolor(xdots/2+ydots/2 , i , currentcolor);
576
putcolor(xdots/2-ydots/2+i , 0 , currentcolor);
577
putcolor(xdots/2-ydots/2+i , ydots-1 , currentcolor);
578
}
579
else
580
for (i=0;i<xdots;i++){
581
putcolor(0 , ydots/2-xdots/2+i , currentcolor);
582
putcolor(xdots-1 , ydots/2-xdots/2+i , currentcolor);
583
putcolor(i , ydots/2-xdots/2 , currentcolor);
584
putcolor(i , ydots/2+xdots/2 , currentcolor);
585
}
586
break;
587
}
588
589
for(;;)
590
{
591
switch (mode) {
592
case 0: /* Release new point on a circle inside the box */
593
angle=2*(double)rand()/(RAND_MAX/PI);
594
FPUsincos(&angle,&sine,&cosine);
595
x = (int)(cosine*(xmax-xmin) + xdots);
596
y = (int)(sine *(ymax-ymin) + ydots);
597
x = x >> 1; /* divide by 2 */
598
y = y >> 1;
599
break;
600
case 1: /* Release new point on the line ymin somewhere between xmin
601
and xmax */
602
y=ymin;
603
x=RANDOM(xmax-xmin) + (xdots-xmax+xmin)/2;
604
break;
605
case 2: /* Release new point on a circle inside the box with radius
606
given by the radius variable */
607
angle=2*(double)rand()/(RAND_MAX/PI);
608
FPUsincos(&angle,&sine,&cosine);
609
x = (int)(cosine*radius + xdots);
610
y = (int)(sine *radius + ydots);
611
x = x >> 1;
612
y = y >> 1;
613
break;
614
}
615
616
/* Loop as long as the point (x,y) is surrounded by color 0 */
617
/* on all eight sides */
618
619
while((getcolor(x+1,y+1) == 0) && (getcolor(x+1,y) == 0) &&
620
(getcolor(x+1,y-1) == 0) && (getcolor(x ,y+1) == 0) &&
621
(getcolor(x ,y-1) == 0) && (getcolor(x-1,y+1) == 0) &&
622
(getcolor(x-1,y) == 0) && (getcolor(x-1,y-1) == 0))
623
{
624
/* Erase moving point */
625
if (show_orbit)
626
putcolor(x,y,0);
627
628
if (mode==0){ /* Make sure point is inside the box */
629
if (x==xmax)
630
x--;
631
else if (x==xmin)
632
x++;
633
if (y==ymax)
634
y--;
635
else if (y==ymin)
636
y++;
637
}
638
639
if (mode==1) /* Make sure point is on the screen below ymin, but
640
we need a 1 pixel margin because of the next random step.*/
641
{
642
if (x>=xdots-1)
643
x--;
644
else if (x<=1)
645
x++;
646
if (y<ymin)
647
y++;
648
}
649
650
/* Take one random step */
651
x += RANDOM(3) - 1;
652
y += RANDOM(3) - 1;
653
654
/* Check keyboard */
655
if ((++plasma_check & 0x7f) == 1)
656
if(check_key())
657
{
658
alloc_resume(20,1);
659
if (mode!=2)
660
put_resume(sizeof(xmax),&xmax,sizeof(xmin),&xmin,
661
sizeof(ymax),&ymax,sizeof(ymin),&ymin,0);
662
else
663
put_resume(sizeof(xmax),&xmax,sizeof(xmin),&xmin,
664
sizeof(ymax),&ymax,sizeof(radius),&radius,0);
665
666
plasma_check--;
667
return 1;
668
}
669
670
/* Show the moving point */
671
if (show_orbit)
672
putcolor(x,y,RANDOM(colors-1)+1);
673
674
} /* End of loop, now fix the point */
675
676
/* If we're doing colorshifting then use currentcolor, otherwise
677
pick one at random */
678
putcolor(x,y,colorshift?currentcolor:RANDOM(colors-1)+1);
679
680
/* If we're doing colorshifting then check to see if we need to shift*/
681
if (colorshift){
682
if (!--colorcount){ /* If the counter reaches zero then shift*/
683
currentcolor++; /* Increase the current color and wrap */
684
currentcolor%=colors; /* around skipping zero */
685
if (!currentcolor) currentcolor++;
686
colorcount=colorshift; /* and reset the counter */
687
}
688
}
689
690
/* If the new point is close to an edge, we may need to increase
691
some limits so that the limits expand to match the growing
692
fractal. */
693
694
switch (mode) {
695
case 0: if (((x+border)>xmax) || ((x-border)<xmin)
696
|| ((y-border)<ymin) || ((y+border)>ymax))
697
{
698
/* Increase box size, but not past the edge of the screen */
699
ymin--;
700
ymax++;
701
xmin--;
702
xmax++;
703
if ((ymin==0) || (xmin==0))
704
return 0;
705
}
706
break;
707
case 1: /* Decrease ymin, but not past top of screen */
708
if (y-border < ymin)
709
ymin--;
710
if (ymin==0)
711
return 0;
712
break;
713
case 2: /* Decrease the radius where points are released to stay away
714
from the fractal. It might be decreased by 1 or 2 */
715
r = sqr((float)x-xdots/2) + sqr((float)y-ydots/2);
716
if (r<=border*border)
717
return 0;
718
while ((radius-border)*(radius-border) > r)
719
radius--;
720
break;
721
}
722
}
723
}
724
725
726
727
/************ standalone engine for "bifurcation" types ***************/
728
729
/***************************************************************/
730
/* The following code now forms a generalised Fractal Engine */
731
/* for Bifurcation fractal typeS. By rights it now belongs in */
732
/* CALCFRACT.C, but it's easier for me to leave it here ! */
733
734
/* Original code by Phil Wilson, hacked around by Kev Allen. */
735
736
/* Besides generalisation, enhancements include Periodicity */
737
/* Checking during the plotting phase (AND halfway through the */
738
/* filter cycle, if possible, to halve calc times), quicker */
739
/* floating-point calculations for the standard Verhulst type, */
740
/* and new bifurcation types (integer bifurcation, f.p & int */
741
/* biflambda - the real equivalent of complex Lambda sets - */
742
/* and f.p renditions of bifurcations of r*sin(Pi*p), which */
743
/* spurred Mitchel Feigenbaum on to discover his Number). */
744
745
/* To add further types, extend the fractalspecific[] array in */
746
/* usual way, with Bifurcation as the engine, and the name of */
747
/* the routine that calculates the next bifurcation generation */
748
/* as the "orbitcalc" routine in the fractalspecific[] entry. */
749
750
/* Bifurcation "orbitcalc" routines get called once per screen */
751
/* pixel column. They should calculate the next generation */
752
/* from the doubles Rate & Population (or the longs lRate & */
753
/* lPopulation if they use integer math), placing the result */
754
/* back in Population (or lPopulation). They should return 0 */
755
/* if all is ok, or any non-zero value if calculation bailout */
756
/* is desirable (eg in case of errors, or the series tending */
757
/* to infinity). Have fun ! */
758
/***************************************************************/
759
760
#define DEFAULTFILTER 1000 /* "Beauty of Fractals" recommends using 5000
761
(p.25), but that seems unnecessary. Can
762
override this value with a nonzero param1 */
763
764
#define SEED 0.66 /* starting value for population */
765
766
static int far *verhulst_array;
767
unsigned long filter_cycles;
768
static unsigned int half_time_check;
769
static long lPopulation, lRate;
770
double Population, Rate;
771
static int mono, outside_x;
772
static long LPI;
773
774
int Bifurcation(void)
775
{
776
unsigned long array_size;
777
int row, column;
778
column = 0;
779
if (resuming)
780
{
781
start_resume();
782
get_resume(sizeof(column),&column,0);
783
end_resume();
784
}
785
array_size = (iystop + 1) * sizeof(int); /* should be iystop + 1 */
786
if ((verhulst_array = (int far *) farmemalloc(array_size)) == NULL)
787
{
788
static FCODE msg[]={"Insufficient free memory for calculation."};
789
stopmsg(0,msg);
790
return(-1);
791
}
792
793
LPI = (long)(PI * fudge);
794
795
for (row = 0; row <= iystop; row++) /* should be iystop */
796
verhulst_array[row] = 0;
797
798
mono = 0;
799
if (colors == 2)
800
mono = 1;
801
if (mono)
802
{
803
if (inside)
804
{
805
outside_x = 0;
806
inside = 1;
807
}
808
else
809
outside_x = 1;
810
}
811
812
filter_cycles = (parm.x <= 0) ? DEFAULTFILTER : (long)parm.x;
813
half_time_check = FALSE;
814
if (periodicitycheck && (unsigned long)maxit < filter_cycles)
815
{
816
filter_cycles = (filter_cycles - maxit + 1) / 2;
817
half_time_check = TRUE;
818
}
819
820
if (integerfractal)
821
linit.y = ymax - iystop*dely; /* Y-value of */
822
else
823
init.y = (double)(yymax - iystop*delyy); /* bottom pixels */
824
825
while (column <= ixstop)
826
{
827
if(keypressed())
828
{
829
farmemfree((char far *)verhulst_array);
830
alloc_resume(10,1);
831
put_resume(sizeof(column),&column,0);
832
return(-1);
833
}
834
835
if (integerfractal)
836
lRate = xmin + column*delx;
837
else
838
Rate = (double)(xxmin + column*delxx);
839
verhulst(); /* calculate array once per column */
840
841
for (row = iystop; row >= 0; row--) /* should be iystop & >=0 */
842
{
843
int color;
844
color = verhulst_array[row];
845
if(color && mono)
846
color = inside;
847
else if((!color) && mono)
848
color = outside_x;
849
else if (color>=colors)
850
color = colors-1;
851
verhulst_array[row] = 0;
852
(*plot)(column,row,color); /* was row-1, but that's not right? */
853
}
854
column++;
855
}
856
farmemfree((char far *)verhulst_array);
857
return(0);
858
}
859
860
static void verhulst() /* P. F. Verhulst (1845) */
861
{
862
unsigned int pixel_row, errors;
863
unsigned long counter;
864
865
if (integerfractal)
866
lPopulation = (parm.y == 0) ? (long)(SEED*fudge) : (long)(parm.y*fudge);
867
else
868
Population = (parm.y == 0 ) ? SEED : parm.y;
869
870
errors = overflow = FALSE;
871
872
for (counter=0 ; counter < filter_cycles ; counter++)
873
{
874
errors = curfractalspecific->orbitcalc();
875
if (errors)
876
return;
877
}
878
if (half_time_check) /* check for periodicity at half-time */
879
{
880
Bif_Period_Init();
881
for (counter=0 ; counter < (unsigned long)maxit ; counter++)
882
{
883
errors = curfractalspecific->orbitcalc();
884
if (errors) return;
885
if (periodicitycheck && Bif_Periodic(counter)) break;
886
}
887
if (counter >= (unsigned long)maxit) /* if not periodic, go the distance */
888
{
889
for (counter=0 ; counter < filter_cycles ; counter++)
890
{
891
errors = curfractalspecific->orbitcalc();
892
if (errors) return;
893
}
894
}
895
}
896
897
if (periodicitycheck) Bif_Period_Init();
898
for (counter=0 ; counter < (unsigned long)maxit ; counter++)
899
{
900
errors = curfractalspecific->orbitcalc();
901
if (errors) return;
902
903
/* assign population value to Y coordinate in pixels */
904
if (integerfractal)
905
pixel_row = iystop - (int)((lPopulation - linit.y) / dely); /* iystop */
906
else
907
pixel_row = iystop - (int)((Population - init.y) / delyy);
908
909
/* if it's visible on the screen, save it in the column array */
910
if (pixel_row <= (unsigned int)iystop) /* JCO 6/6/92 */
911
verhulst_array[ pixel_row ] ++;
912
if (periodicitycheck && Bif_Periodic(counter))
913
{
914
if (pixel_row <= (unsigned int)iystop) /* JCO 6/6/92 */
915
verhulst_array[ pixel_row ] --;
916
break;
917
}
918
}
919
}
920
static long lBif_closenuf, lBif_savedpop; /* poss future use */
921
static double Bif_closenuf, Bif_savedpop;
922
static int Bif_savedinc;
923
static long Bif_savedand;
924
925
static void Bif_Period_Init()
926
{
927
Bif_savedinc = 1;
928
Bif_savedand = 1;
929
if (integerfractal)
930
{
931
lBif_savedpop = -1;
932
lBif_closenuf = dely / 8;
933
}
934
else
935
{
936
Bif_savedpop = -1.0;
937
Bif_closenuf = (double)delyy / 8.0;
938
}
939
}
940
941
static int _fastcall Bif_Periodic (long time) /* Bifurcation Population Periodicity Check */
942
/* Returns : 1 if periodicity found, else 0 */
943
{
944
if ((time & Bif_savedand) == 0) /* time to save a new value */
945
{
946
if (integerfractal) lBif_savedpop = lPopulation;
947
else Bif_savedpop = Population;
948
if (--Bif_savedinc == 0)
949
{
950
Bif_savedand = (Bif_savedand << 1) + 1;
951
Bif_savedinc = 4;
952
}
953
}
954
else /* check against an old save */
955
{
956
if (integerfractal)
957
{
958
if (labs(lBif_savedpop-lPopulation) <= lBif_closenuf)
959
return(1);
960
}
961
else
962
{
963
if (fabs(Bif_savedpop-Population) <= Bif_closenuf)
964
return(1);
965
}
966
}
967
return(0);
968
}
969
970
/**********************************************************************/
971
/* */
972
/* The following are Bifurcation "orbitcalc" routines... */
973
/* */
974
/**********************************************************************/
975
#ifdef XFRACT
976
int BifurcLambda() /* Used by lyanupov */
977
{
978
Population = Rate * Population * (1 - Population);
979
return (fabs(Population) > BIG);
980
}
981
#endif
982
983
/* Modified formulas below to generalize bifurcations. JCO 7/3/92 */
984
985
#define LCMPLXtrig0(arg,out) Arg1->l = (arg); ltrig0(); (out)=Arg1->l
986
#define CMPLXtrig0(arg,out) Arg1->d = (arg); dtrig0(); (out)=Arg1->d
987
988
int BifurcVerhulstTrig()
989
{
990
/* Population = Pop + Rate * fn(Pop) * (1 - fn(Pop)) */
991
tmp.x = Population;
992
tmp.y = 0;
993
CMPLXtrig0(tmp, tmp);
994
Population += Rate * tmp.x * (1 - tmp.x);
995
return (fabs(Population) > BIG);
996
}
997
998
int LongBifurcVerhulstTrig()
999
{
1000
#ifndef XFRACT
1001
ltmp.x = lPopulation;
1002
ltmp.y = 0;
1003
LCMPLXtrig0(ltmp, ltmp);
1004
ltmp.y = ltmp.x - multiply(ltmp.x,ltmp.x,bitshift);
1005
lPopulation += multiply(lRate,ltmp.y,bitshift);
1006
#endif
1007
return (overflow);
1008
}
1009
1010
int BifurcStewartTrig()
1011
{
1012
/* Population = (Rate * fn(Population) * fn(Population)) - 1.0 */
1013
tmp.x = Population;
1014
tmp.y = 0;
1015
CMPLXtrig0(tmp, tmp);
1016
Population = (Rate * tmp.x * tmp.x) - 1.0;
1017
return (fabs(Population) > BIG);
1018
}
1019
1020
int LongBifurcStewartTrig()
1021
{
1022
#ifndef XFRACT
1023
ltmp.x = lPopulation;
1024
ltmp.y = 0;
1025
LCMPLXtrig0(ltmp, ltmp);
1026
lPopulation = multiply(ltmp.x,ltmp.x,bitshift);
1027
lPopulation = multiply(lPopulation,lRate, bitshift);
1028
lPopulation -= fudge;
1029
#endif
1030
return (overflow);
1031
}
1032
1033
int BifurcSetTrigPi()
1034
{
1035
tmp.x = Population * PI;
1036
tmp.y = 0;
1037
CMPLXtrig0(tmp, tmp);
1038
Population = Rate * tmp.x;
1039
return (fabs(Population) > BIG);
1040
}
1041
1042
int LongBifurcSetTrigPi()
1043
{
1044
#ifndef XFRACT
1045
ltmp.x = multiply(lPopulation,LPI,bitshift);
1046
ltmp.y = 0;
1047
LCMPLXtrig0(ltmp, ltmp);
1048
lPopulation = multiply(lRate,ltmp.x,bitshift);
1049
#endif
1050
return (overflow);
1051
}
1052
1053
int BifurcAddTrigPi()
1054
{
1055
tmp.x = Population * PI;
1056
tmp.y = 0;
1057
CMPLXtrig0(tmp, tmp);
1058
Population += Rate * tmp.x;
1059
return (fabs(Population) > BIG);
1060
}
1061
1062
int LongBifurcAddTrigPi()
1063
{
1064
#ifndef XFRACT
1065
ltmp.x = multiply(lPopulation,LPI,bitshift);
1066
ltmp.y = 0;
1067
LCMPLXtrig0(ltmp, ltmp);
1068
lPopulation += multiply(lRate,ltmp.x,bitshift);
1069
#endif
1070
return (overflow);
1071
}
1072
1073
int BifurcLambdaTrig()
1074
{
1075
/* Population = Rate * fn(Population) * (1 - fn(Population)) */
1076
tmp.x = Population;
1077
tmp.y = 0;
1078
CMPLXtrig0(tmp, tmp);
1079
Population = Rate * tmp.x * (1 - tmp.x);
1080
return (fabs(Population) > BIG);
1081
}
1082
1083
int LongBifurcLambdaTrig()
1084
{
1085
#ifndef XFRACT
1086
ltmp.x = lPopulation;
1087
ltmp.y = 0;
1088
LCMPLXtrig0(ltmp, ltmp);
1089
ltmp.y = ltmp.x - multiply(ltmp.x,ltmp.x,bitshift);
1090
lPopulation = multiply(lRate,ltmp.y,bitshift);
1091
#endif
1092
return (overflow);
1093
}
1094
1095
#define LCMPLXpwr(arg1,arg2,out) Arg2->l = (arg1); Arg1->l = (arg2);\
1096
lStkPwr(); Arg1++; Arg2++; (out) = Arg2->l
1097
1098
long beta;
1099
1100
int BifurcMay()
1101
{ /* X = (lambda * X) / (1 + X)^beta, from R.May as described in Pickover,
1102
Computers, Pattern, Chaos, and Beauty, page 153 */
1103
tmp.x = 1.0 + Population;
1104
tmp.x = pow(tmp.x, -beta); /* pow in math.h included with mpmath.h */
1105
Population = (Rate * Population) * tmp.x;
1106
return (fabs(Population) > BIG);
1107
}
1108
1109
int LongBifurcMay()
1110
{
1111
#ifndef XFRACT
1112
ltmp.x = lPopulation + fudge;
1113
ltmp.y = 0;
1114
lparm2.x = beta * fudge;
1115
LCMPLXpwr(ltmp, lparm2, ltmp);
1116
lPopulation = multiply(lRate,lPopulation,bitshift);
1117
lPopulation = divide(lPopulation,ltmp.x,bitshift);
1118
#endif
1119
return (overflow);
1120
}
1121
1122
int BifurcMaySetup()
1123
{
1124
1125
beta = (long)param[2];
1126
if(beta < 2)
1127
beta = 2;
1128
param[2] = (double)beta;
1129
1130
timer(0,curfractalspecific->calctype);
1131
return(0);
1132
}
1133
1134
/* Here Endeth the Generalised Bifurcation Fractal Engine */
1135
1136
/* END Phil Wilson's Code (modified slightly by Kev Allen et. al. !) */
1137
1138
1139
/******************* standalone engine for "popcorn" ********************/
1140
1141
int popcorn() /* subset of std engine */
1142
{
1143
int start_row;
1144
start_row = 0;
1145
if (resuming)
1146
{
1147
start_resume();
1148
get_resume(sizeof(start_row),&start_row,0);
1149
end_resume();
1150
}
1151
kbdcount=max_kbdcount;
1152
plot = noplot;
1153
tempsqrx = ltempsqrx = 0; /* PB added this to cover weird BAILOUTs */
1154
for (row = start_row; row <= iystop; row++)
1155
{
1156
reset_periodicity = 1;
1157
for (col = 0; col <= ixstop; col++)
1158
{
1159
if (StandardFractal() == -1) /* interrupted */
1160
{
1161
alloc_resume(10,1);
1162
put_resume(sizeof(row),&row,0);
1163
return(-1);
1164
}
1165
reset_periodicity = 0;
1166
}
1167
}
1168
calc_status = 4;
1169
return(0);
1170
}
1171
1172
/******************* standalone engine for "lyapunov" *********************/
1173
/*** Roy Murphy [76376,721] ***/
1174
/*** revision history: ***/
1175
/*** initial version: Winter '91 ***/
1176
/*** Fall '92 integration of Nicholas Wilt's ASM speedups ***/
1177
/*** Jan 93' integration with calcfrac() yielding boundary tracing, ***/
1178
/*** tesseral, and solid guessing, and inversion, inside=nnn ***/
1179
/*** save_release behavior: ***/
1180
/*** 1730 & prior: ignores inside=, calcmode='1', (a,b)->(x,y) ***/
1181
/*** 1731: other calcmodes and inside=nnn ***/
1182
/*** 1732: the infamous axis swap: (b,a)->(x,y), ***/
1183
/*** the order parameter becomes a long int ***/
1184
/**************************************************************************/
1185
int lyaLength, lyaSeedOK;
1186
int lyaRxy[34];
1187
1188
#define WES 1 /* define WES to be 0 to use Nick's lyapunov.obj */
1189
#if WES
1190
int lyapunov_cycles(double, double);
1191
#else
1192
int lyapunov_cycles(int, double, double, double);
1193
#endif
1194
1195
int lyapunov_cycles_in_c(long, double, double);
1196
1197
int lyapunov () {
1198
double a, b;
1199
1200
if (keypressed()) {
1201
return -1;
1202
}
1203
overflow=FALSE;
1204
if (param[1]==1) Population = (1.0+rand())/(2.0+RAND_MAX);
1205
else if (param[1]==0) {
1206
if (fabs(Population)>BIG || Population==0 || Population==1)
1207
Population = (1.0+rand())/(2.0+RAND_MAX);
1208
}
1209
else Population = param[1];
1210
(*plot)(col, row, 1);
1211
if (invert) {
1212
invertz2(&init);
1213
a = init.y;
1214
b = init.x;
1215
}
1216
else {
1217
a = dypixel();
1218
b = dxpixel();
1219
}
1220
#ifndef XFRACT
1221
/* the assembler routines don't work for a & b outside the
1222
ranges 0 < a < 4 and 0 < b < 4. So, fall back on the C
1223
routines if part of the image sticks out.
1224
*/
1225
#if WES
1226
color=lyapunov_cycles(a, b);
1227
#else
1228
if (lyaSeedOK && a>0 && b>0 && a<=4 && b<=4)
1229
color=lyapunov_cycles(filter_cycles, Population, a, b);
1230
else
1231
color=lyapunov_cycles_in_c(filter_cycles, a, b);
1232
#endif
1233
#else
1234
color=lyapunov_cycles_in_c(filter_cycles, a, b);
1235
#endif
1236
if (inside>0 && color==0)
1237
color = inside;
1238
else if (color>=colors)
1239
color = colors-1;
1240
(*plot)(col, row, color);
1241
return color;
1242
}
1243
1244
1245
int lya_setup () {
1246
/* This routine sets up the sequence for forcing the Rate parameter
1247
to vary between the two values. It fills the array lyaRxy[] and
1248
sets lyaLength to the length of the sequence.
1249
1250
The sequence is coded in the bit pattern in an integer.
1251
Briefly, the sequence starts with an A the leading zero bits
1252
are ignored and the remaining bit sequence is decoded. The
1253
sequence ends with a B. Not all possible sequences can be
1254
represented in this manner, but every possible sequence is
1255
either represented as itself, as a rotation of one of the
1256
representable sequences, or as the inverse of a representable
1257
sequence (swapping 0s and 1s in the array.) Sequences that
1258
are the rotation and/or inverses of another sequence will generate
1259
the same lyapunov exponents.
1260
1261
A few examples follow:
1262
number sequence
1263
0 ab
1264
1 aab
1265
2 aabb
1266
3 aaab
1267
4 aabbb
1268
5 aabab
1269
6 aaabb (this is a duplicate of 4, a rotated inverse)
1270
7 aaaab
1271
8 aabbbb etc.
1272
*/
1273
1274
long i;
1275
int t;
1276
1277
if ((filter_cycles=(long)param[2])==0)
1278
filter_cycles=maxit/2;
1279
lyaSeedOK = param[1]>0 && param[1]<=1 && debugflag!=90;
1280
lyaLength = 1;
1281
1282
i = (long)param[0];
1283
#ifndef XFRACT
1284
if (save_release<1732) i &= 0x0FFFFL; /* make it a short to reporduce prior stuff*/
1285
#endif
1286
lyaRxy[0] = 1;
1287
for (t=31; t>=0; t--)
1288
if (i & (1<<t)) break;
1289
for (; t>=0; t--)
1290
lyaRxy[lyaLength++] = (i & (1<<t)) != 0;
1291
lyaRxy[lyaLength++] = 0;
1292
if (save_release<1732) /* swap axes prior to 1732 */
1293
for (t=lyaLength; t>=0; t--)
1294
lyaRxy[t] = !lyaRxy[t];
1295
if (save_release<1731) { /* ignore inside=, stdcalcmode */
1296
stdcalcmode='1';
1297
if (inside==1) inside = 0;
1298
}
1299
if (inside<0) {
1300
static FCODE msg[]=
1301
{"Sorry, inside options other than inside=nnn are not supported by the lyapunov"};
1302
stopmsg(0,(char far *)msg);
1303
inside=1;
1304
}
1305
if (usr_stdcalcmode == 'o') { /* Oops,lyapunov type */
1306
usr_stdcalcmode = '1'; /* doesn't use new & breaks orbits */
1307
stdcalcmode = '1';
1308
}
1309
return 1;
1310
}
1311
1312
int lyapunov_cycles_in_c(long filter_cycles, double a, double b) {
1313
int color, count, lnadjust;
1314
long i;
1315
double lyap, total, temp;
1316
/* e10=22026.4657948 e-10=0.0000453999297625 */
1317
1318
total = 1.0;
1319
lnadjust = 0;
1320
for (i=0; i<filter_cycles; i++) {
1321
for (count=0; count<lyaLength; count++) {
1322
Rate = lyaRxy[count] ? a : b;
1323
if (curfractalspecific->orbitcalc()) {
1324
overflow = TRUE;
1325
goto jumpout;
1326
}
1327
}
1328
}
1329
for (i=0; i < maxit/2; i++) {
1330
for (count = 0; count < lyaLength; count++) {
1331
Rate = lyaRxy[count] ? a : b;
1332
if (curfractalspecific->orbitcalc()) {
1333
overflow = TRUE;
1334
goto jumpout;
1335
}
1336
temp = fabs(Rate-2.0*Rate*Population);
1337
if ((total *= (temp))==0) {
1338
overflow = TRUE;
1339
goto jumpout;
1340
}
1341
}
1342
while (total > 22026.4657948) {
1343
total *= 0.0000453999297625;
1344
lnadjust += 10;
1345
}
1346
while (total < 0.0000453999297625) {
1347
total *= 22026.4657948;
1348
lnadjust -= 10;
1349
}
1350
}
1351
1352
jumpout:
1353
if (overflow || total <= 0 || (temp = log(total) + lnadjust) > 0)
1354
color = 0;
1355
else {
1356
if (LogFlag)
1357
lyap = -temp/((double) lyaLength*i);
1358
else
1359
lyap = 1 - exp(temp/((double) lyaLength*i));
1360
color = 1 + (int)(lyap * (colors-1));
1361
}
1362
return color;
1363
}
1364
1365
1366
/******************* standalone engine for "cellular" ********************/
1367
/* Originally coded by Ken Shirriff.
1368
Modified beyond recognition by Jonathan Osuch.
1369
Original or'd the neighborhood, changed to sum the neighborhood
1370
Changed prompts and error messages
1371
Added CA types
1372
Set the palette to some standard? CA colors
1373
Changed *cell_array to near and used dstack so put_line and get_line
1374
could be used all the time
1375
Made space bar generate next screen
1376
Increased string/rule size to 16 digits and added CA types 9/20/92
1377
*/
1378
1379
#define BAD_T 1
1380
#define BAD_MEM 2
1381
#define STRING1 3
1382
#define STRING2 4
1383
#define TABLEK 5
1384
#define TYPEKR 6
1385
#define RULELENGTH 7
1386
#define INTERUPT 8
1387
1388
#define CELLULAR_DONE 10
1389
1390
#ifndef XFRACT
1391
static BYTE *cell_array[2];
1392
#else
1393
static BYTE far *cell_array[2];
1394
#endif
1395
1396
S16 r, k_1, rule_digits;
1397
int lstscreenflag;
1398
1399
void abort_cellular(int err, int t)
1400
{
1401
int i;
1402
switch (err)
1403
{
1404
case BAD_T:
1405
{
1406
char msg[30];
1407
sprintf(msg,"Bad t=%d, aborting\n", t);
1408
stopmsg(0,(char far *)msg);
1409
}
1410
break;
1411
case BAD_MEM:
1412
{
1413
static FCODE msg[]={"Insufficient free memory for calculation" };
1414
stopmsg(0,msg);
1415
}
1416
break;
1417
case STRING1:
1418
{
1419
static FCODE msg[]={"String can be a maximum of 16 digits" };
1420
stopmsg(0,msg);
1421
}
1422
break;
1423
case STRING2:
1424
{
1425
static FCODE msg[]={"Make string of 0's through 's" };
1426
msg[27] = (char)(k_1 + 48); /* turn into a character value */
1427
stopmsg(0,msg);
1428
}
1429
break;
1430
case TABLEK:
1431
{
1432
static FCODE msg[]={"Make Rule with 0's through 's" };
1433
msg[27] = (char)(k_1 + 48); /* turn into a character value */
1434
stopmsg(0,msg);
1435
}
1436
break;
1437
case TYPEKR:
1438
{
1439
static FCODE msg[]={"Type must be 21, 31, 41, 51, 61, 22, 32, \
1440
42, 23, 33, 24, 25, 26, 27" };
1441
stopmsg(0,msg);
1442
}
1443
break;
1444
case RULELENGTH:
1445
{
1446
static FCODE msg[]={"Rule must be digits long" };
1447
i = rule_digits / 10;
1448
if(i==0)
1449
msg[14] = (char)(rule_digits + 48);
1450
else {
1451
msg[13] = (char)(i+48);
1452
msg[14] = (char)((rule_digits % 10) + 48);
1453
}
1454
stopmsg(0,msg);
1455
}
1456
break;
1457
case INTERUPT:
1458
{
1459
static FCODE msg[]={"Interrupted, can't resume" };
1460
stopmsg(0,msg);
1461
}
1462
break;
1463
case CELLULAR_DONE:
1464
break;
1465
}
1466
if(cell_array[0] != NULL)
1467
#ifndef XFRACT
1468
cell_array[0] = NULL;
1469
#else
1470
farmemfree((char far *)cell_array[0]);
1471
#endif
1472
if(cell_array[1] != NULL)
1473
#ifndef XFRACT
1474
cell_array[1] = NULL;
1475
#else
1476
farmemfree((char far *)cell_array[1]);
1477
#endif
1478
}
1479
1480
int cellular () {
1481
S16 start_row;
1482
S16 filled, notfilled;
1483
U16 cell_table[32];
1484
U16 init_string[16];
1485
U16 kr, k;
1486
U32 lnnmbr;
1487
U16 i, twor;
1488
S16 t, t2;
1489
S32 randparam;
1490
double n;
1491
char buf[30];
1492
1493
set_Cellular_palette();
1494
1495
randparam = (S32)param[0];
1496
lnnmbr = (U32)param[3];
1497
kr = (U16)param[2];
1498
switch (kr) {
1499
case 21:
1500
case 31:
1501
case 41:
1502
case 51:
1503
case 61:
1504
case 22:
1505
case 32:
1506
case 42:
1507
case 23:
1508
case 33:
1509
case 24:
1510
case 25:
1511
case 26:
1512
case 27:
1513
break;
1514
default:
1515
abort_cellular(TYPEKR, 0);
1516
return -1;
1517
/* break; */
1518
}
1519
1520
r = (S16)(kr % 10); /* Number of nearest neighbors to sum */
1521
k = (U16)(kr / 10); /* Number of different states, k=3 has states 0,1,2 */
1522
k_1 = (S16)(k - 1); /* Highest state value, k=3 has highest state value of 2 */
1523
rule_digits = (S16)((r * 2 + 1) * k_1 + 1); /* Number of digits in the rule */
1524
1525
if ((!rflag) && randparam == -1)
1526
--rseed;
1527
if (randparam != 0 && randparam != -1) {
1528
n = param[0];
1529
sprintf(buf,"%.16g",n); /* # of digits in initial string */
1530
t = (S16)strlen(buf);
1531
if (t>16 || t <= 0) {
1532
abort_cellular(STRING1, 0);
1533
return -1;
1534
}
1535
for (i=0;i<16;i++)
1536
init_string[i] = 0; /* zero the array */
1537
t2 = (S16) ((16 - t)/2);
1538
for (i=0;i<(U16)t;i++) { /* center initial string in array */
1539
init_string[i+t2] = (U16)(buf[i] - 48); /* change character to number */
1540
if (init_string[i+t2]>(U16)k_1) {
1541
abort_cellular(STRING2, 0);
1542
return -1;
1543
}
1544
}
1545
}
1546
1547
srand(rseed);
1548
if (!rflag) ++rseed;
1549
1550
/* generate rule table from parameter 1 */
1551
#ifndef XFRACT
1552
n = param[1];
1553
#else
1554
/* gcc can't manage to convert a big double to an unsigned long properly. */
1555
if (param[1]>0x7fffffff) {
1556
n = (param[1]-0x7fffffff);
1557
n += 0x7fffffff;
1558
} else {
1559
n = param[1];
1560
}
1561
#endif
1562
if (n == 0) { /* calculate a random rule */
1563
n = rand()%(int)k;
1564
for (i=1;i<(U16)rule_digits;i++) {
1565
n *= 10;
1566
n += rand()%(int)k;
1567
}
1568
param[1] = n;
1569
}
1570
sprintf(buf,"%.*g",rule_digits ,n);
1571
t = (S16)strlen(buf);
1572
if (rule_digits < t || t < 0) { /* leading 0s could make t smaller */
1573
abort_cellular(RULELENGTH, 0);
1574
return -1;
1575
}
1576
for (i=0;i<(U16)rule_digits;i++) /* zero the table */
1577
cell_table[i] = 0;
1578
for (i=0;i<(U16)t;i++) { /* reverse order */
1579
cell_table[i] = (U16)(buf[t-i-1] - 48); /* change character to number */
1580
if (cell_table[i]>(U16)k_1) {
1581
abort_cellular(TABLEK, 0);
1582
return -1;
1583
}
1584
}
1585
1586
1587
start_row = 0;
1588
#ifndef XFRACT
1589
/* two 4096 byte arrays, at present at most 2024 + 1 bytes should be */
1590
/* needed in each array (max screen width + 1) */
1591
cell_array[0] = (BYTE *)&dstack[0]; /* dstack is in general.asm */
1592
cell_array[1] = (BYTE *)&boxy[0]; /* boxy is in general.asm */
1593
#else
1594
cell_array[0] = (BYTE far *)farmemalloc(ixstop+1);
1595
cell_array[1] = (BYTE far *)farmemalloc(ixstop+1);
1596
#endif
1597
if (cell_array[0]==NULL || cell_array[1]==NULL) {
1598
abort_cellular(BAD_MEM, 0);
1599
return(-1);
1600
}
1601
1602
/* nxtscreenflag toggled by space bar in fractint.c, 1 for continuous */
1603
/* 0 to stop on next screen */
1604
1605
filled = 0;
1606
notfilled = (S16)(1-filled);
1607
if (resuming && !nxtscreenflag && !lstscreenflag) {
1608
start_resume();
1609
get_resume(sizeof(start_row),&start_row,0);
1610
end_resume();
1611
get_line(start_row,0,ixstop,cell_array[filled]);
1612
}
1613
else if (nxtscreenflag && !lstscreenflag) {
1614
start_resume();
1615
end_resume();
1616
get_line(iystop,0,ixstop,cell_array[filled]);
1617
param[3] += iystop + 1;
1618
start_row = -1; /* after 1st iteration its = 0 */
1619
}
1620
else {
1621
if(rflag || randparam==0 || randparam==-1){
1622
for (col=0;col<=ixstop;col++) {
1623
cell_array[filled][col] = (BYTE)(rand()%(int)k);
1624
}
1625
} /* end of if random */
1626
1627
else {
1628
for (col=0;col<=ixstop;col++) { /* Clear from end to end */
1629
cell_array[filled][col] = 0;
1630
}
1631
i = 0;
1632
for (col=(ixstop-16)/2;col<(ixstop+16)/2;col++) { /* insert initial */
1633
cell_array[filled][col] = (BYTE)init_string[i++]; /* string */
1634
}
1635
} /* end of if not random */
1636
if (lnnmbr != 0)
1637
lstscreenflag = 1;
1638
else
1639
lstscreenflag = 0;
1640
put_line(start_row,0,ixstop,cell_array[filled]);
1641
}
1642
start_row++;
1643
1644
/* This section calculates the starting line when it is not zero */
1645
/* This section can't be resumed since no screen output is generated */
1646
/* calculates the (lnnmbr - 1) generation */
1647
if (lstscreenflag) { /* line number != 0 & not resuming & not continuing */
1648
U32 big_row;
1649
for (big_row = (U32)start_row; big_row < lnnmbr; big_row++) {
1650
static FCODE msg[]={"Cellular thinking (higher start row takes longer)"};
1651
1652
thinking(1,msg);
1653
if(rflag || randparam==0 || randparam==-1){
1654
/* Use a random border */
1655
for (i=0;i<=(U16)r;i++) {
1656
cell_array[notfilled][i]=(BYTE)(rand()%(int)k);
1657
cell_array[notfilled][ixstop-i]=(BYTE)(rand()%(int)k);
1658
}
1659
}
1660
else {
1661
/* Use a zero border */
1662
for (i=0;i<=(U16)r;i++) {
1663
cell_array[notfilled][i]=0;
1664
cell_array[notfilled][ixstop-i]=0;
1665
}
1666
}
1667
1668
t = 0; /* do first cell */
1669
for (twor=(U16)(r+r),i=0;i<=twor;i++)
1670
t = (S16)(t + (S16)cell_array[filled][i]);
1671
if (t>rule_digits || t<0) {
1672
thinking(0, NULL);
1673
abort_cellular(BAD_T, t);
1674
return(-1);
1675
}
1676
cell_array[notfilled][r] = (BYTE)cell_table[t];
1677
1678
/* use a rolling sum in t */
1679
for (col=r+1;col<ixstop-r;col++) { /* now do the rest */
1680
t = (S16)(t + cell_array[filled][col+r] - cell_array[filled][col-r-1]);
1681
if (t>rule_digits || t<0) {
1682
thinking(0, NULL);
1683
abort_cellular(BAD_T, t);
1684
return(-1);
1685
}
1686
cell_array[notfilled][col] = (BYTE)cell_table[t];
1687
}
1688
1689
filled = notfilled;
1690
notfilled = (S16)(1-filled);
1691
if (keypressed()) {
1692
thinking(0, NULL);
1693
abort_cellular(INTERUPT, 0);
1694
return -1;
1695
}
1696
}
1697
start_row = 0;
1698
thinking(0, NULL);
1699
lstscreenflag = 0;
1700
}
1701
1702
/* This section does all the work */
1703
contloop:
1704
for (row = start_row; row <= iystop; row++) {
1705
1706
if(rflag || randparam==0 || randparam==-1){
1707
/* Use a random border */
1708
for (i=0;i<=(U16)r;i++) {
1709
cell_array[notfilled][i]=(BYTE)(rand()%(int)k);
1710
cell_array[notfilled][ixstop-i]=(BYTE)(rand()%(int)k);
1711
}
1712
}
1713
else {
1714
/* Use a zero border */
1715
for (i=0;i<=(U16)r;i++) {
1716
cell_array[notfilled][i]=0;
1717
cell_array[notfilled][ixstop-i]=0;
1718
}
1719
}
1720
1721
t = 0; /* do first cell */
1722
for (twor=(U16)(r+r),i=0;i<=twor;i++)
1723
t = (S16)(t + (S16)cell_array[filled][i]);
1724
if (t>rule_digits || t<0) {
1725
thinking(0, NULL);
1726
abort_cellular(BAD_T, t);
1727
return(-1);
1728
}
1729
cell_array[notfilled][r] = (BYTE)cell_table[t];
1730
1731
/* use a rolling sum in t */
1732
for (col=r+1;col<ixstop-r;col++) { /* now do the rest */
1733
t = (S16)(t + cell_array[filled][col+r] - cell_array[filled][col-r-1]);
1734
if (t>rule_digits || t<0) {
1735
thinking(0, NULL);
1736
abort_cellular(BAD_T, t);
1737
return(-1);
1738
}
1739
cell_array[notfilled][col] = (BYTE)cell_table[t];
1740
}
1741
1742
filled = notfilled;
1743
notfilled = (S16)(1-filled);
1744
put_line(row,0,ixstop,cell_array[filled]);
1745
if (keypressed()) {
1746
abort_cellular(CELLULAR_DONE, 0);
1747
alloc_resume(10,1);
1748
put_resume(sizeof(row),&row,0);
1749
return -1;
1750
}
1751
}
1752
if(nxtscreenflag) {
1753
param[3] += iystop + 1;
1754
start_row = -1; /* after 1st iteration its = 0 */
1755
goto contloop;
1756
}
1757
abort_cellular(CELLULAR_DONE, 0);
1758
return 1;
1759
}
1760
1761
int CellularSetup(void)
1762
{
1763
if (!resuming) {
1764
nxtscreenflag = 0; /* initialize flag */
1765
}
1766
timer(0,curfractalspecific->calctype);
1767
return(0);
1768
}
1769
1770
static void set_Cellular_palette()
1771
{
1772
static Palettetype Red = { 42, 0, 0 };
1773
static Palettetype Green = { 10,35,10 };
1774
static Palettetype Blue = { 13,12,29 };
1775
static Palettetype Yellow = { 60,58,18 };
1776
static Palettetype Brown = { 42,21, 0 };
1777
1778
if (mapdacbox) return; /* map= specified */
1779
1780
dac[0].red = 0 ;
1781
dac[0].green= 0 ;
1782
dac[0].blue = 0 ;
1783
1784
dac[1].red = Red.red;
1785
dac[1].green = Red.green;
1786
dac[1].blue = Red.blue;
1787
1788
dac[2].red = Green.red;
1789
dac[2].green = Green.green;
1790
dac[2].blue = Green.blue;
1791
1792
dac[3].red = Blue.red;
1793
dac[3].green = Blue.green;
1794
dac[3].blue = Blue.blue;
1795
1796
dac[4].red = Yellow.red;
1797
dac[4].green = Yellow.green;
1798
dac[4].blue = Yellow.blue;
1799
1800
dac[5].red = Brown.red;
1801
dac[5].green = Brown.green;
1802
dac[5].blue = Brown.blue;
1803
1804
SetTgaColors();
1805
spindac(0,1);
1806
}
1807
1808
/* frothy basin routines */
1809
1810
#define FROTH_BITSHIFT 28
1811
#define FROTH_D_TO_L(x) ((long)((x)*(1L<<FROTH_BITSHIFT)))
1812
#define FROTH_CLOSE 1e-6 /* seems like a good value */
1813
#define FROTH_LCLOSE FROTH_D_TO_L(FROTH_CLOSE)
1814
#define SQRT3 1.732050807568877193
1815
#define FROTH_SLOPE SQRT3
1816
#define FROTH_LSLOPE FROTH_D_TO_L(FROTH_SLOPE)
1817
#define FROTH_CRITICAL_A 1.028713768218725 /* 1.0287137682187249127 */
1818
#define froth_top_x_mapping(x) ((x)*(x)-(x)-3*fsp->fl.f.a*fsp->fl.f.a/4)
1819
1820
1821
static char froth3_256c[] = "froth3.map";
1822
static char froth6_256c[] = "froth6.map";
1823
static char froth3_16c[] = "froth316.map";
1824
static char froth6_16c[] = "froth616.map";
1825
1826
struct froth_double_struct {
1827
double a;
1828
double halfa;
1829
double top_x1;
1830
double top_x2;
1831
double top_x3;
1832
double top_x4;
1833
double left_x1;
1834
double left_x2;
1835
double left_x3;
1836
double left_x4;
1837
double right_x1;
1838
double right_x2;
1839
double right_x3;
1840
double right_x4;
1841
};
1842
1843
struct froth_long_struct {
1844
long a;
1845
long halfa;
1846
long top_x1;
1847
long top_x2;
1848
long top_x3;
1849
long top_x4;
1850
long left_x1;
1851
long left_x2;
1852
long left_x3;
1853
long left_x4;
1854
long right_x1;
1855
long right_x2;
1856
long right_x3;
1857
long right_x4;
1858
};
1859
1860
struct froth_struct {
1861
int repeat_mapping;
1862
int altcolor;
1863
int attractors;
1864
int shades;
1865
union { /* This was made into a union to save 56 malloc()'ed bytes. */
1866
struct froth_double_struct f;
1867
struct froth_long_struct l;
1868
} fl;
1869
};
1870
1871
struct froth_struct *fsp=NULL; /* froth_struct pointer */
1872
1873
/* color maps which attempt to replicate the images of James Alexander. */
1874
static void set_Froth_palette(void)
1875
{
1876
char *mapname;
1877
1878
if (colorstate != 0) /* 0 means dacbox matches default */
1879
return;
1880
if (colors >= 16)
1881
{
1882
if (colors >= 256)
1883
{
1884
if (fsp->attractors == 6)
1885
mapname = froth6_256c;
1886
else
1887
mapname = froth3_256c;
1888
}
1889
else /* colors >= 16 */
1890
{
1891
if (fsp->attractors == 6)
1892
mapname = froth6_16c;
1893
else
1894
mapname = froth3_16c;
1895
}
1896
if (ValidateLuts(mapname) != 0)
1897
return;
1898
colorstate = 0; /* treat map it as default */
1899
spindac(0,1);
1900
}
1901
}
1902
1903
int froth_setup(void)
1904
{
1905
double sin_theta, cos_theta, x0, y0;
1906
1907
sin_theta = SQRT3/2; /* sin(2*PI/3) */
1908
cos_theta = -0.5; /* cos(2*PI/3) */
1909
1910
/* check for NULL as safety net */
1911
if (fsp == NULL)
1912
fsp = (struct froth_struct *)malloc(sizeof (struct froth_struct));
1913
if (fsp == NULL)
1914
{
1915
static FCODE msg[]=
1916
{"Sorry, not enough memory to run the frothybasin fractal type"};
1917
stopmsg(0,(char far *)msg);
1918
return 0;
1919
}
1920
1921
/* for the all important backwards compatibility */
1922
if (save_release <= 1821) /* book version is 18.21 */
1923
{
1924
/* use old release parameters */
1925
1926
fsp->repeat_mapping = ((int)param[0] == 6 || (int)param[0] == 2); /* map 1 or 2 times (3 or 6 basins) */
1927
fsp->altcolor = (int)param[1];
1928
param[2] = 0; /* throw away any value used prior to 18.20 */
1929
1930
fsp->attractors = !fsp->repeat_mapping ? 3 : 6;
1931
1932
/* use old values */ /* old names */
1933
fsp->fl.f.a = 1.02871376822; /* A */
1934
fsp->fl.f.halfa = fsp->fl.f.a/2; /* A/2 */
1935
1936
fsp->fl.f.top_x1 = -1.04368901270; /* X1MIN */
1937
fsp->fl.f.top_x2 = 1.33928675524; /* X1MAX */
1938
fsp->fl.f.top_x3 = -0.339286755220; /* XMIDT */
1939
fsp->fl.f.top_x4 = -0.339286755220; /* XMIDT */
1940
1941
fsp->fl.f.left_x1 = 0.07639837810; /* X3MAX2 */
1942
fsp->fl.f.left_x2 = -1.11508950586; /* X2MIN2 */
1943
fsp->fl.f.left_x3 = -0.27580275066; /* XMIDL */
1944
fsp->fl.f.left_x4 = -0.27580275066; /* XMIDL */
1945
1946
fsp->fl.f.right_x1 = 0.96729063460; /* X2MAX1 */
1947
fsp->fl.f.right_x2 = -0.22419724936; /* X3MIN1 */
1948
fsp->fl.f.right_x3 = 0.61508950585; /* XMIDR */
1949
fsp->fl.f.right_x4 = 0.61508950585; /* XMIDR */
1950
1951
}
1952
else /* use new code */
1953
{
1954
if (param[0] != 2)
1955
param[0] = 1;
1956
fsp->repeat_mapping = (int)param[0] == 2;
1957
if (param[1] != 0)
1958
param[1] = 1;
1959
fsp->altcolor = (int)param[1];
1960
fsp->fl.f.a = param[2];
1961
1962
fsp->attractors = fabs(fsp->fl.f.a) <= FROTH_CRITICAL_A ? (!fsp->repeat_mapping ? 3 : 6)
1963
: (!fsp->repeat_mapping ? 2 : 3);
1964
1965
/* new improved values */
1966
/* 0.5 is the value that causes the mapping to reach a minimum */
1967
x0 = 0.5;
1968
/* a/2 is the value that causes the y value to be invariant over the mappings */
1969
y0 = fsp->fl.f.halfa = fsp->fl.f.a/2;
1970
fsp->fl.f.top_x1 = froth_top_x_mapping(x0);
1971
fsp->fl.f.top_x2 = froth_top_x_mapping(fsp->fl.f.top_x1);
1972
fsp->fl.f.top_x3 = froth_top_x_mapping(fsp->fl.f.top_x2);
1973
fsp->fl.f.top_x4 = froth_top_x_mapping(fsp->fl.f.top_x3);
1974
1975
/* rotate 120 degrees counter-clock-wise */
1976
fsp->fl.f.left_x1 = fsp->fl.f.top_x1 * cos_theta - y0 * sin_theta;
1977
fsp->fl.f.left_x2 = fsp->fl.f.top_x2 * cos_theta - y0 * sin_theta;
1978
fsp->fl.f.left_x3 = fsp->fl.f.top_x3 * cos_theta - y0 * sin_theta;
1979
fsp->fl.f.left_x4 = fsp->fl.f.top_x4 * cos_theta - y0 * sin_theta;
1980
1981
/* rotate 120 degrees clock-wise */
1982
fsp->fl.f.right_x1 = fsp->fl.f.top_x1 * cos_theta + y0 * sin_theta;
1983
fsp->fl.f.right_x2 = fsp->fl.f.top_x2 * cos_theta + y0 * sin_theta;
1984
fsp->fl.f.right_x3 = fsp->fl.f.top_x3 * cos_theta + y0 * sin_theta;
1985
fsp->fl.f.right_x4 = fsp->fl.f.top_x4 * cos_theta + y0 * sin_theta;
1986
1987
}
1988
1989
/* if 2 attractors, use same shades as 3 attractors */
1990
fsp->shades = (colors-1) / max(3,fsp->attractors);
1991
1992
/* rqlim needs to be at least sq(1+sqrt(1+sq(a))), */
1993
/* which is never bigger than 6.93..., so we'll call it 7.0 */
1994
if (rqlim < 7.0)
1995
rqlim=7.0;
1996
set_Froth_palette();
1997
/* make the best of the .map situation */
1998
orbit_color = fsp->attractors != 6 && colors >= 16 ? (fsp->shades<<1)+1 : colors-1;
1999
2000
if (integerfractal)
2001
{
2002
struct froth_long_struct tmp_l;
2003
2004
tmp_l.a = FROTH_D_TO_L(fsp->fl.f.a);
2005
tmp_l.halfa = FROTH_D_TO_L(fsp->fl.f.halfa);
2006
2007
tmp_l.top_x1 = FROTH_D_TO_L(fsp->fl.f.top_x1);
2008
tmp_l.top_x2 = FROTH_D_TO_L(fsp->fl.f.top_x2);
2009
tmp_l.top_x3 = FROTH_D_TO_L(fsp->fl.f.top_x3);
2010
tmp_l.top_x4 = FROTH_D_TO_L(fsp->fl.f.top_x4);
2011
2012
tmp_l.left_x1 = FROTH_D_TO_L(fsp->fl.f.left_x1);
2013
tmp_l.left_x2 = FROTH_D_TO_L(fsp->fl.f.left_x2);
2014
tmp_l.left_x3 = FROTH_D_TO_L(fsp->fl.f.left_x3);
2015
tmp_l.left_x4 = FROTH_D_TO_L(fsp->fl.f.left_x4);
2016
2017
tmp_l.right_x1 = FROTH_D_TO_L(fsp->fl.f.right_x1);
2018
tmp_l.right_x2 = FROTH_D_TO_L(fsp->fl.f.right_x2);
2019
tmp_l.right_x3 = FROTH_D_TO_L(fsp->fl.f.right_x3);
2020
tmp_l.right_x4 = FROTH_D_TO_L(fsp->fl.f.right_x4);
2021
2022
fsp->fl.l = tmp_l;
2023
}
2024
return 1;
2025
}
2026
2027
void froth_cleanup(void)
2028
{
2029
if (fsp != NULL)
2030
free(fsp);
2031
/* set to NULL as a flag that froth_cleanup() has been called */
2032
fsp = NULL;
2033
}
2034
2035
2036
/* Froth Fractal type */
2037
int calcfroth(void) /* per pixel 1/2/g, called with row & col set */
2038
{
2039
int found_attractor=0;
2040
2041
if (check_key()) {
2042
return -1;
2043
}
2044
2045
if (fsp == NULL)
2046
{ /* error occured allocating memory for fsp */
2047
return 0;
2048
}
2049
2050
orbit_ptr = 0;
2051
coloriter = 0;
2052
if(showdot>0)
2053
(*plot) (col, row, showdot%colors);
2054
if (!integerfractal) /* fp mode */
2055
{
2056
if(invert)
2057
{
2058
invertz2(&tmp);
2059
old = tmp;
2060
}
2061
else
2062
{
2063
old.x = dxpixel();
2064
old.y = dypixel();
2065
}
2066
2067
while (!found_attractor
2068
&& ((tempsqrx=sqr(old.x)) + (tempsqry=sqr(old.y)) < rqlim)
2069
&& (coloriter < maxit))
2070
{
2071
/* simple formula: z = z^2 + conj(z*(-1+ai)) */
2072
/* but it's the attractor that makes this so interesting */
2073
new.x = tempsqrx - tempsqry - old.x - fsp->fl.f.a*old.y;
2074
old.y += (old.x+old.x)*old.y - fsp->fl.f.a*old.x;
2075
old.x = new.x;
2076
if (fsp->repeat_mapping)
2077
{
2078
new.x = sqr(old.x) - sqr(old.y) - old.x - fsp->fl.f.a*old.y;
2079
old.y += (old.x+old.x)*old.y - fsp->fl.f.a*old.x;
2080
old.x = new.x;
2081
}
2082
2083
coloriter++;
2084
2085
if (show_orbit) {
2086
if (keypressed())
2087
break;
2088
plot_orbit(old.x, old.y, -1);
2089
}
2090
2091
if (fabs(fsp->fl.f.halfa-old.y) < FROTH_CLOSE
2092
&& old.x >= fsp->fl.f.top_x1 && old.x <= fsp->fl.f.top_x2)
2093
{
2094
if ((!fsp->repeat_mapping && fsp->attractors == 2)
2095
|| (fsp->repeat_mapping && fsp->attractors == 3))
2096
found_attractor = 1;
2097
else if (old.x <= fsp->fl.f.top_x3)
2098
found_attractor = 1;
2099
else if (old.x >= fsp->fl.f.top_x4) {
2100
if (!fsp->repeat_mapping)
2101
found_attractor = 1;
2102
else
2103
found_attractor = 2;
2104
}
2105
}
2106
else if (fabs(FROTH_SLOPE*old.x - fsp->fl.f.a - old.y) < FROTH_CLOSE
2107
&& old.x <= fsp->fl.f.right_x1 && old.x >= fsp->fl.f.right_x2)
2108
{
2109
if (!fsp->repeat_mapping && fsp->attractors == 2)
2110
found_attractor = 2;
2111
else if (fsp->repeat_mapping && fsp->attractors == 3)
2112
found_attractor = 3;
2113
else if (old.x >= fsp->fl.f.right_x3) {
2114
if (!fsp->repeat_mapping)
2115
found_attractor = 2;
2116
else
2117
found_attractor = 4;
2118
}
2119
else if (old.x <= fsp->fl.f.right_x4) {
2120
if (!fsp->repeat_mapping)
2121
found_attractor = 3;
2122
else
2123
found_attractor = 6;
2124
}
2125
}
2126
else if (fabs(-FROTH_SLOPE*old.x - fsp->fl.f.a - old.y) < FROTH_CLOSE
2127
&& old.x <= fsp->fl.f.left_x1 && old.x >= fsp->fl.f.left_x2)
2128
{
2129
if (!fsp->repeat_mapping && fsp->attractors == 2)
2130
found_attractor = 2;
2131
else if (fsp->repeat_mapping && fsp->attractors == 3)
2132
found_attractor = 2;
2133
else if (old.x >= fsp->fl.f.left_x3) {
2134
if (!fsp->repeat_mapping)
2135
found_attractor = 3;
2136
else
2137
found_attractor = 5;
2138
}
2139
else if (old.x <= fsp->fl.f.left_x4) {
2140
if (!fsp->repeat_mapping)
2141
found_attractor = 2;
2142
else
2143
found_attractor = 3;
2144
}
2145
}
2146
}
2147
}
2148
else /* integer mode */
2149
{
2150
if(invert)
2151
{
2152
invertz2(&tmp);
2153
lold.x = (long)(tmp.x * fudge);
2154
lold.y = (long)(tmp.y * fudge);
2155
}
2156
else
2157
{
2158
lold.x = lxpixel();
2159
lold.y = lypixel();
2160
}
2161
2162
while (!found_attractor && ((lmagnitud = (ltempsqrx=lsqr(lold.x)) + (ltempsqry=lsqr(lold.y))) < llimit)
2163
&& (lmagnitud >= 0) && (coloriter < maxit))
2164
{
2165
/* simple formula: z = z^2 + conj(z*(-1+ai)) */
2166
/* but it's the attractor that makes this so interesting */
2167
lnew.x = ltempsqrx - ltempsqry - lold.x - multiply(fsp->fl.l.a,lold.y,bitshift);
2168
lold.y += (multiply(lold.x,lold.y,bitshift)<<1) - multiply(fsp->fl.l.a,lold.x,bitshift);
2169
lold.x = lnew.x;
2170
if (fsp->repeat_mapping)
2171
{
2172
lmagnitud = (ltempsqrx=lsqr(lold.x)) + (ltempsqry=lsqr(lold.y));
2173
if ((lmagnitud > llimit) || (lmagnitud < 0))
2174
break;
2175
lnew.x = ltempsqrx - ltempsqry - lold.x - multiply(fsp->fl.l.a,lold.y,bitshift);
2176
lold.y += (multiply(lold.x,lold.y,bitshift)<<1) - multiply(fsp->fl.l.a,lold.x,bitshift);
2177
lold.x = lnew.x;
2178
}
2179
coloriter++;
2180
2181
if (show_orbit) {
2182
if (keypressed())
2183
break;
2184
iplot_orbit(lold.x, lold.y, -1);
2185
}
2186
2187
if (labs(fsp->fl.l.halfa-lold.y) < FROTH_LCLOSE
2188
&& lold.x > fsp->fl.l.top_x1 && lold.x < fsp->fl.l.top_x2)
2189
{
2190
if ((!fsp->repeat_mapping && fsp->attractors == 2)
2191
|| (fsp->repeat_mapping && fsp->attractors == 3))
2192
found_attractor = 1;
2193
else if (lold.x <= fsp->fl.l.top_x3)
2194
found_attractor = 1;
2195
else if (lold.x >= fsp->fl.l.top_x4) {
2196
if (!fsp->repeat_mapping)
2197
found_attractor = 1;
2198
else
2199
found_attractor = 2;
2200
}
2201
}
2202
else if (labs(multiply(FROTH_LSLOPE,lold.x,bitshift)-fsp->fl.l.a-lold.y) < FROTH_LCLOSE
2203
&& lold.x <= fsp->fl.l.right_x1 && lold.x >= fsp->fl.l.right_x2)
2204
{
2205
if (!fsp->repeat_mapping && fsp->attractors == 2)
2206
found_attractor = 2;
2207
else if (fsp->repeat_mapping && fsp->attractors == 3)
2208
found_attractor = 3;
2209
else if (lold.x >= fsp->fl.l.right_x3) {
2210
if (!fsp->repeat_mapping)
2211
found_attractor = 2;
2212
else
2213
found_attractor = 4;
2214
}
2215
else if (lold.x <= fsp->fl.l.right_x4) {
2216
if (!fsp->repeat_mapping)
2217
found_attractor = 3;
2218
else
2219
found_attractor = 6;
2220
}
2221
}
2222
else if (labs(multiply(-FROTH_LSLOPE,lold.x,bitshift)-fsp->fl.l.a-lold.y) < FROTH_LCLOSE)
2223
{
2224
if (!fsp->repeat_mapping && fsp->attractors == 2)
2225
found_attractor = 2;
2226
else if (fsp->repeat_mapping && fsp->attractors == 3)
2227
found_attractor = 2;
2228
else if (lold.x >= fsp->fl.l.left_x3) {
2229
if (!fsp->repeat_mapping)
2230
found_attractor = 3;
2231
else
2232
found_attractor = 5;
2233
}
2234
else if (lold.x <= fsp->fl.l.left_x4) {
2235
if (!fsp->repeat_mapping)
2236
found_attractor = 2;
2237
else
2238
found_attractor = 3;
2239
}
2240
}
2241
}
2242
}
2243
if (show_orbit)
2244
scrub_orbit();
2245
2246
realcoloriter = coloriter;
2247
if ((kbdcount -= abs((int)realcoloriter)) <= 0)
2248
{
2249
if (check_key())
2250
return (-1);
2251
kbdcount = max_kbdcount;
2252
}
2253
2254
/* inside - Here's where non-palette based images would be nice. Instead, */
2255
/* we'll use blocks of (colors-1)/3 or (colors-1)/6 and use special froth */
2256
/* color maps in attempt to replicate the images of James Alexander. */
2257
if (found_attractor)
2258
{
2259
if (colors >= 256)
2260
{
2261
if (!fsp->altcolor)
2262
{
2263
if (coloriter > fsp->shades)
2264
coloriter = fsp->shades;
2265
}
2266
else
2267
coloriter = fsp->shades * coloriter / maxit;
2268
if (coloriter == 0)
2269
coloriter = 1;
2270
coloriter += fsp->shades * (found_attractor-1);
2271
}
2272
else if (colors >= 16)
2273
{ /* only alternate coloring scheme available for 16 colors */
2274
long lshade;
2275
2276
/* Trying to make a better 16 color distribution. */
2277
/* Since their are only a few possiblities, just handle each case. */
2278
/* This is a mostly guess work here. */
2279
lshade = (coloriter<<16)/maxit;
2280
if (fsp->attractors != 6) /* either 2 or 3 attractors */
2281
{
2282
if (lshade < 2622) /* 0.04 */
2283
coloriter = 1;
2284
else if (lshade < 10486) /* 0.16 */
2285
coloriter = 2;
2286
else if (lshade < 23593) /* 0.36 */
2287
coloriter = 3;
2288
else if (lshade < 41943L) /* 0.64 */
2289
coloriter = 4;
2290
else
2291
coloriter = 5;
2292
coloriter += 5 * (found_attractor-1);
2293
}
2294
else /* 6 attractors */
2295
{
2296
if (lshade < 10486) /* 0.16 */
2297
coloriter = 1;
2298
else
2299
coloriter = 2;
2300
coloriter += 2 * (found_attractor-1);
2301
}
2302
}
2303
else /* use a color corresponding to the attractor */
2304
coloriter = found_attractor;
2305
oldcoloriter = coloriter;
2306
}
2307
else /* outside, or inside but didn't get sucked in by attractor. */
2308
coloriter = 0;
2309
2310
color = abs((int)(coloriter));
2311
2312
(*plot)(col, row, color);
2313
2314
return color;
2315
}
2316
2317
/*
2318
These last two froth functions are for the orbit-in-window feature.
2319
Normally, this feature requires StandardFractal, but since it is the
2320
attractor that makes the frothybasin type so unique, it is worth
2321
putting in as a stand-alone.
2322
*/
2323
2324
int froth_per_pixel(void)
2325
{
2326
if (!integerfractal) /* fp mode */
2327
{
2328
old.x = dxpixel();
2329
old.y = dypixel();
2330
tempsqrx=sqr(old.x);
2331
tempsqry=sqr(old.y);
2332
}
2333
else /* integer mode */
2334
{
2335
lold.x = lxpixel();
2336
lold.y = lypixel();
2337
ltempsqrx = multiply(lold.x, lold.x, bitshift);
2338
ltempsqry = multiply(lold.y, lold.y, bitshift);
2339
}
2340
return 0;
2341
}
2342
2343
int froth_per_orbit(void)
2344
{
2345
if (!integerfractal) /* fp mode */
2346
{
2347
new.x = tempsqrx - tempsqry - old.x - fsp->fl.f.a*old.y;
2348
new.y = 2.0*old.x*old.y - fsp->fl.f.a*old.x + old.y;
2349
if (fsp->repeat_mapping)
2350
{
2351
old = new;
2352
new.x = sqr(old.x) - sqr(old.y) - old.x - fsp->fl.f.a*old.y;
2353
new.y = 2.0*old.x*old.y - fsp->fl.f.a*old.x + old.y;
2354
}
2355
2356
if ((tempsqrx=sqr(new.x)) + (tempsqry=sqr(new.y)) >= rqlim)
2357
return 1;
2358
old = new;
2359
}
2360
else /* integer mode */
2361
{
2362
lnew.x = ltempsqrx - ltempsqry - lold.x - multiply(fsp->fl.l.a,lold.y,bitshift);
2363
lnew.y = lold.y + (multiply(lold.x,lold.y,bitshift)<<1) - multiply(fsp->fl.l.a,lold.x,bitshift);
2364
if (fsp->repeat_mapping)
2365
{
2366
if ((ltempsqrx=lsqr(lnew.x)) + (ltempsqry=lsqr(lnew.y)) >= llimit)
2367
return 1;
2368
lold = lnew;
2369
lnew.x = ltempsqrx - ltempsqry - lold.x - multiply(fsp->fl.l.a,lold.y,bitshift);
2370
lnew.y = lold.y + (multiply(lold.x,lold.y,bitshift)<<1) - multiply(fsp->fl.l.a,lold.x,bitshift);
2371
}
2372
if ((ltempsqrx=lsqr(lnew.x)) + (ltempsqry=lsqr(lnew.y)) >= llimit)
2373
return 1;
2374
lold = lnew;
2375
}
2376
return 0;
2377
}
2378
Loggerhead is a web-based interface for Breezy
Version: 2.0.1