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