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