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