File: common\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
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
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);
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);
1232 #endif
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];
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;
1471 #endif
1472 if(cell_array[1] != NULL)
1473 #ifndef XFRACT
1474 cell_array[1] = NULL;
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];
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 */
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 && colorstate != 0) 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 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
2379