home *** CD-ROM | disk | FTP | other *** search
/ Media Share 9 / MEDIASHARE_09.ISO / progmisc / frasrc18.zip / MISCFRAC.C < prev    next >
C/C++ Source or Header  |  1993-03-28  |  58KB  |  2,151 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 suffix[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.    {
  536.       notdiskmsg();
  537.       return(-1);
  538.    }
  539.  
  540.    bitshift = 16;
  541.    fudge = 1L << 16;
  542.  
  543.    border = param[0];
  544.    mode = param[1];
  545.    if (mode > 2)
  546.       mode=0;
  547.  
  548.    if (border <= 0)
  549.       border = 10;
  550.  
  551.    srand(rseed);
  552.    if (!rflag) ++rseed;
  553.  
  554.    if (mode == 0) {
  555.       xmax = xdots / 2 + border;  /* Initial box */
  556.       xmin = xdots / 2 - border;
  557.       ymax = ydots / 2 + border;
  558.       ymin = ydots / 2 - border;
  559.    }
  560.    if (mode == 1) {
  561.       xmax = xdots / 2 + border;  /* Initial box */
  562.       xmin = xdots / 2 - border;
  563.       ymin = ydots - 20;
  564.    }
  565.    if (mode == 2) {
  566.       if (xdots>ydots)
  567.          radius = ydots - 5;
  568.       else
  569.          radius = xdots - 5;
  570.    }
  571.    if (resuming) /* restore worklist, if we can't the above will stay in place */
  572.    {
  573.       start_resume();
  574.       if (mode != 2)
  575.          get_resume(sizeof(int),&xmax,sizeof(int),&xmin,sizeof(int),&ymax,
  576.              sizeof(int),&ymin,0);
  577.       else
  578.          get_resume(sizeof(int),&xmax,sizeof(int),&xmin,sizeof(int),&ymax,
  579.              sizeof(int),&radius,0);
  580.       end_resume();
  581.    }
  582.  
  583.    if (mode==0)
  584.       putcolor(xdots / 2, ydots / 2,RANDOM(colors-1)+1);  /* Seed point */
  585.  
  586.    if (mode==1)
  587.       for (i=0;i<=xdots;i++)
  588.          putcolor(i,ydots-1,colors-1);
  589.  
  590.    if (mode==2){
  591.       if (xdots>ydots){
  592.          for (i=0;i<ydots;i++){
  593.             putcolor(xdots/2-ydots/2 , i , colors-1);
  594.             putcolor(xdots/2+ydots/2 , i , colors-1);
  595.             putcolor(xdots/2-ydots/2+i , 0 , colors-1);
  596.             putcolor(xdots/2-ydots/2+i , ydots-1 , colors-1);
  597.          }
  598.       }
  599.       else {
  600.          for (i=0;i<xdots;i++){
  601.             putcolor(0 , ydots/2-xdots/2+i , colors-1);
  602.             putcolor(xdots-1 , ydots/2-xdots/2+i , colors-1);
  603.             putcolor(i , ydots/2-xdots/2 , colors-1);
  604.             putcolor(i , ydots/2+xdots/2 , colors-1);
  605.          }
  606.       }
  607.    }
  608.  
  609.    while (1)
  610.    {
  611.       /* Release new point on circle just inside the box */
  612.  
  613.       if (mode==0)
  614.       {
  615.          if (floatflag)
  616.          {
  617.             angle=2*(double)rand()/(RAND_MAX/PI);
  618.             FPUsincos(&angle,&sine,&cosine);
  619.             x = cosine*(xmax-xmin) + xdots;
  620.             y = sine  *(ymax-ymin) + ydots;
  621.          }
  622.          else
  623.          {
  624.             SinCos086(multiply((long)rand15(),FOURPI,16),&lsine,&lcosine);
  625.             x = (lcosine*(long)(xmax-xmin) >> 16) + xdots;
  626.             y = (lsine  *(long)(ymax-ymin) >> 16) + ydots;
  627.          }
  628.  
  629.          x = x >> 1; /* divide by 2 */
  630.          y = y >> 1;
  631.       }
  632.       if (mode==1)
  633.       {
  634.          y=ymin;
  635.          x=RANDOM(xmax-xmin) + (xdots-xmax+xmin)/2;
  636.       }
  637.       if (mode==2)
  638.       {
  639.          if (floatflag)
  640.          {
  641.             angle=2*(double)rand()/(RAND_MAX/PI);
  642.             FPUsincos(&angle,&sine,&cosine);
  643.             x = cosine*radius + xdots;
  644.             y = sine  *radius + ydots;
  645.          }
  646.          else
  647.          {
  648.             SinCos086(multiply((long)rand15(),FOURPI,16),&lsine,&lcosine);
  649.             x = (lcosine*(long)(radius) >> 16) + xdots;
  650.             y = (lsine  *(long)(radius) >> 16) + ydots;
  651.          }
  652.          x = x >> 1;
  653.          y = y >> 1;
  654.       }
  655.       /* Loop as long as the point (x,y) is surrounded by color 0 */
  656.       /* on all eight sides                                       */
  657.       while((getcolor(x+1,y+1) == 0) && (getcolor(x+1,y) == 0) &&
  658.           (getcolor(x+1,y-1) == 0) && (getcolor(x  ,y+1) == 0) &&
  659.           (getcolor(x  ,y-1) == 0) && (getcolor(x-1,y+1) == 0) &&
  660.           (getcolor(x-1,y) == 0) && (getcolor(x-1,y-1) == 0))
  661.       {
  662.          /* Erase moving point */
  663.          if (show_orbit)
  664.             putcolor(x,y,0);
  665.  
  666.          /* Make sure point is inside the box (if mode==0)*/
  667.          if (mode==0){
  668.             if (x==xmax)
  669.                x--;
  670.             else if (x==xmin)
  671.                x++;
  672.             if (y==ymax)
  673.                y--;
  674.             else if (y==ymin)
  675.                y++;
  676.          }
  677.  
  678.          if (mode==1)
  679.          {
  680.             if (x>xdots-2)
  681.                x--;
  682.             else if (x<1)
  683.                x++;
  684.             if (y<ymin)
  685.                y++;
  686.          }
  687.  
  688.          /* Take one random step */
  689.          x += RANDOM(3) - 1;
  690.          y += RANDOM(3) - 1;
  691.  
  692.          /* Check keyboard */
  693.          if ((++plasma_check & 0x7f) == 1)
  694.             if(check_key())
  695.             {
  696.                alloc_resume(20,1);
  697.                if (mode!=2)
  698.                   put_resume(sizeof(int),&xmax,sizeof(int),&xmin, sizeof(int),&ymax,
  699.                       sizeof(int),&ymin,0);
  700.                else
  701.                   put_resume(sizeof(int),&xmax,sizeof(int),&xmin, sizeof(int),&ymax,
  702.                       sizeof(int),&radius,0);
  703.  
  704.                plasma_check--;
  705.                return 1;
  706.             }
  707.  
  708.          /* Show the moving point */
  709.          if (show_orbit)
  710.             putcolor(x,y,RANDOM(colors-1)+1);
  711.  
  712.       }
  713.       putcolor(x,y,RANDOM(colors-1)+1);
  714.  
  715.       /* Is point too close to the edge? */
  716.  
  717.       if (mode==0)
  718.       {
  719.          if (((x+border)>xmax) || ((x-border)<xmin)
  720.              || ((y-border)<ymin) || ((y+border)>ymax))
  721.          {
  722.             /* Increase box size, but not past the edge of the screen */
  723.             if (ymin != 1)
  724.             {
  725.                ymin--;
  726.                ymax++;
  727.             }
  728.             if (xmin != 1)
  729.             {
  730.                xmin--;
  731.                xmax++;
  732.             }
  733.             if ((ymin==1) || (xmin==1))
  734.                return 0;
  735.          }
  736.       }
  737.       if (mode==1)
  738.       {
  739.          if (y < ymin+5)
  740.             ymin = y - 5;
  741.          if (ymin<2)
  742.             return 0;
  743.       }
  744.       if (mode==2)
  745.       {
  746.          if (abs(x-xdots/2)<5 && abs(y-ydots/2)<5)
  747.             return 0;
  748.  
  749.          r = (x-xdots/2)*(x-xdots/2)+(y-ydots/2)*(y-ydots/2);
  750.          fSqrt14(r,f_tmp);
  751.          r = 2 * f_tmp;
  752.          if (r < radius)
  753.             radius = r;
  754.       }
  755.    }
  756. }
  757.  
  758. /************ standalone engine for "bifurcation" types ***************/
  759.  
  760. /***************************************************************/
  761. /* The following code now forms a generalised Fractal Engine   */
  762. /* for Bifurcation fractal typeS.  By rights it now belongs in */
  763. /* CALCFRACT.C, but it's easier for me to leave it here !      */
  764.  
  765. /* Original code by Phil Wilson, hacked around by Kev Allen.   */
  766.  
  767. /* Besides generalisation, enhancements include Periodicity    */
  768. /* Checking during the plotting phase (AND halfway through the */
  769. /* filter cycle, if possible, to halve calc times), quicker    */
  770. /* floating-point calculations for the standard Verhulst type, */
  771. /* and new bifurcation types (integer bifurcation, f.p & int   */
  772. /* biflambda - the real equivalent of complex Lambda sets -    */
  773. /* and f.p renditions of bifurcations of r*sin(Pi*p), which    */
  774. /* spurred Mitchel Feigenbaum on to discover his Number).      */
  775.  
  776. /* To add further types, extend the fractalspecific[] array in */
  777. /* usual way, with Bifurcation as the engine, and the name of  */
  778. /* the routine that calculates the next bifurcation generation */
  779. /* as the "orbitcalc" routine in the fractalspecific[] entry.  */
  780.  
  781. /* Bifurcation "orbitcalc" routines get called once per screen */
  782. /* pixel column.  They should calculate the next generation    */
  783. /* from the doubles Rate & Population (or the longs lRate &    */
  784. /* lPopulation if they use integer math), placing the result   */
  785. /* back in Population (or lPopulation).  They should return 0  */
  786. /* if all is ok, or any non-zero value if calculation bailout  */
  787. /* is desirable (eg in case of errors, or the series tending   */
  788. /* to infinity).        Have fun !               */
  789. /***************************************************************/
  790.  
  791. #define DEFAULTFILTER 1000     /* "Beauty of Fractals" recommends using 5000
  792.                    (p.25), but that seems unnecessary. Can
  793.                    override this value with a nonzero param1 */
  794.  
  795. #define SEED 0.66        /* starting value for population */
  796.  
  797. static int far *verhulst_array;
  798. unsigned int filter_cycles;
  799. static unsigned int half_time_check;
  800. static long   lPopulation, lRate;
  801. double Population,  Rate;
  802. static int    mono, outside_x;
  803. static long   LPI;
  804.  
  805. int Bifurcation(void)
  806. {
  807.    unsigned long array_size;
  808.    int row, column;
  809.    column = 0;
  810.    if (resuming)
  811.    {
  812.       start_resume();
  813.       get_resume(sizeof(int),&column,0);
  814.       end_resume();
  815.    }
  816.    array_size =  (iystop + 1) * sizeof(int); /* should be iystop + 1 */
  817.    if ((verhulst_array = (int far *) farmemalloc(array_size)) == NULL)
  818.    {
  819.       static char far msg[]={"Insufficient free memory for calculation."};
  820.       stopmsg(0,msg);
  821.       return(-1);
  822.    }
  823.  
  824.    LPI = PI * fudge;
  825.  
  826.    for (row = 0; row <= iystop; row++) /* should be iystop */
  827.       verhulst_array[row] = 0;
  828.  
  829.    mono = 0;
  830.    if (colors == 2)
  831.       mono = 1;
  832.    if (mono)
  833.    {
  834.       if (inside)
  835.       {
  836.      outside_x = 0;
  837.      inside = 1;
  838.       }
  839.       else
  840.      outside_x = 1;
  841.    }
  842.  
  843.    filter_cycles = (parm.x <= 0) ? DEFAULTFILTER : parm.x;
  844.    half_time_check = FALSE;
  845.    if (periodicitycheck && maxit < filter_cycles)
  846.    {
  847.       filter_cycles = (filter_cycles - maxit + 1) / 2;
  848.       half_time_check = TRUE;
  849.    }
  850.  
  851.    if (integerfractal)
  852.       linit.y = ly0[iystop];   /* Y-value of    */
  853.    else
  854.       init.y = dy0[iystop];   /* bottom pixels */
  855.  
  856.    while (column <= ixstop)
  857.    {
  858.       if(check_key())
  859.       {
  860.      farmemfree((char far *)verhulst_array);
  861.      alloc_resume(10,1);
  862.      put_resume(sizeof(int),&column,0);
  863.      return(-1);
  864.       }
  865.       if (integerfractal)
  866.      lRate = lx0[column];
  867.       else
  868.      Rate = dx0[column];
  869.       verhulst();     /* calculate array once per column */
  870.  
  871.       for (row = iystop; row >= 0; row--) /* should be iystop & >=0 */
  872.       {
  873.      int color;
  874.      color = verhulst_array[row];
  875.      if(color && mono)
  876.         color = inside;
  877.      else if((!color) && mono)
  878.         color = outside_x;
  879.      else if (color>=colors)
  880.         color = colors-1;
  881.      verhulst_array[row] = 0;
  882.      (*plot)(column,row,color); /* was row-1, but that's not right? */
  883.       }
  884.       column++;
  885.    }
  886.    farmemfree((char far *)verhulst_array);
  887.    return(0);
  888. }
  889.  
  890. static void verhulst()        /* P. F. Verhulst (1845) */
  891. {
  892.    unsigned int pixel_row, counter, errors;
  893.  
  894.     if (integerfractal)
  895.        lPopulation = (parm.y == 0) ? SEED * fudge : parm.y * fudge;
  896.     else
  897.        Population = (parm.y == 0 ) ? SEED : parm.y;
  898.  
  899.    errors = overflow = FALSE;
  900.  
  901.    for (counter=0 ; counter < filter_cycles ; counter++)
  902.    {
  903.       errors = (*(curfractalspecific->orbitcalc))();
  904.       if (errors)
  905.      return;
  906.    }
  907.    if (half_time_check) /* check for periodicity at half-time */
  908.    {
  909.       Bif_Period_Init();
  910.       for (counter=0 ; counter < maxit ; counter++)
  911.       {
  912.      errors = (*(curfractalspecific->orbitcalc))();
  913.      if (errors) return;
  914.      if (periodicitycheck && Bif_Periodic(counter)) break;
  915.       }
  916.       if (counter >= maxit)   /* if not periodic, go the distance */
  917.       {
  918.      for (counter=0 ; counter < filter_cycles ; counter++)
  919.      {
  920.         errors = (*(curfractalspecific->orbitcalc))();
  921.         if (errors) return;
  922.      }
  923.       }
  924.    }
  925.  
  926.    if (periodicitycheck) Bif_Period_Init();
  927.    for (counter=0 ; counter < maxit ; counter++)
  928.    {
  929.       errors = (*(curfractalspecific->orbitcalc))();
  930.       if (errors) return;
  931.  
  932.       /* assign population value to Y coordinate in pixels */
  933.       if (integerfractal)
  934.      pixel_row = iystop - (lPopulation - linit.y) / dely; /* iystop */
  935.       else
  936.      pixel_row = iystop - (int)((Population - init.y) / deltaY);
  937.  
  938.       /* if it's visible on the screen, save it in the column array */
  939.       if (pixel_row <= iystop) /* JCO 6/6/92 */
  940.      verhulst_array[ pixel_row ] ++;
  941.       if (periodicitycheck && Bif_Periodic(counter))
  942.       {
  943.      if (pixel_row <= iystop) /* JCO 6/6/92 */
  944.         verhulst_array[ pixel_row ] --;
  945.      break;
  946.       }
  947.    }
  948. }
  949. static    long    lBif_closenuf, lBif_savedpop;    /* poss future use  */
  950. static    double     Bif_closenuf,    Bif_savedpop;
  951. static    int     Bif_savedinc,    Bif_savedand;
  952.  
  953. static void Bif_Period_Init()
  954. {
  955.    Bif_savedinc = 1;
  956.    Bif_savedand = 1;
  957.    if (integerfractal)
  958.    {
  959.       lBif_savedpop = -1;
  960.       lBif_closenuf = dely / 8;
  961.    }
  962.    else
  963.    {
  964.       Bif_savedpop = -1.0;
  965.       Bif_closenuf = deltaY / 8.0;
  966.    }
  967. }
  968.  
  969. static int _fastcall Bif_Periodic (time)  /* Bifurcation Population Periodicity Check */
  970. int time;        /* Returns : 1 if periodicity found, else 0 */
  971. {
  972.    if ((time & Bif_savedand) == 0)    /* time to save a new value */
  973.    {
  974.       if (integerfractal) lBif_savedpop = lPopulation;
  975.       else             Bif_savedpop =  Population;
  976.       if (--Bif_savedinc == 0)
  977.       {
  978.      Bif_savedand = (Bif_savedand << 1) + 1;
  979.      Bif_savedinc = 4;
  980.       }
  981.    }
  982.    else             /* check against an old save */
  983.    {
  984.       if (integerfractal)
  985.       {
  986.      if (labs(lBif_savedpop-lPopulation) <= lBif_closenuf)
  987.         return(1);
  988.       }
  989.       else
  990.       {
  991.      if (fabs(Bif_savedpop-Population) <= Bif_closenuf)
  992.         return(1);
  993.       }
  994.    }
  995.    return(0);
  996. }
  997.  
  998. /**********************************************************************/
  999. /*                                                      */
  1000. /* The following are Bifurcation "orbitcalc" routines...              */
  1001. /*                                                      */
  1002. /**********************************************************************/
  1003. #ifdef XFRACT
  1004. int BifurcLambda() /* Used by lyanupov */
  1005.   {
  1006.     Population = Rate * Population * (1 - Population);
  1007.     return (fabs(Population) > BIG);
  1008.   }
  1009. #endif
  1010.  
  1011. /* Modified formulas below to generalize bifurcations. JCO 7/3/92 */
  1012.  
  1013. #define LCMPLXtrig0(arg,out) Arg1->l = (arg); ltrig0(); (out)=Arg1->l
  1014. #define  CMPLXtrig0(arg,out) Arg1->d = (arg); dtrig0(); (out)=Arg1->d
  1015.  
  1016. int BifurcVerhulstTrig()
  1017.   {
  1018. /*  Population = Pop + Rate * fn(Pop) * (1 - fn(Pop)) */
  1019.     tmp.x = Population;
  1020.     tmp.y = 0;
  1021.     CMPLXtrig0(tmp, tmp);
  1022.     Population += Rate * tmp.x * (1 - tmp.x);
  1023.     return (fabs(Population) > BIG);
  1024.   }
  1025.  
  1026. int LongBifurcVerhulstTrig()
  1027.   {
  1028. #ifndef XFRACT
  1029.     ltmp.x = lPopulation;
  1030.     ltmp.y = 0;
  1031.     LCMPLXtrig0(ltmp, ltmp);
  1032.     ltmp.y = ltmp.x - multiply(ltmp.x,ltmp.x,bitshift);
  1033.     lPopulation += multiply(lRate,ltmp.y,bitshift);
  1034.     return (overflow);
  1035. #endif
  1036.   }
  1037.  
  1038. int BifurcStewartTrig()
  1039.   {
  1040. /*  Population = (Rate * fn(Population) * fn(Population)) - 1.0 */
  1041.     tmp.x = Population;
  1042.     tmp.y = 0;
  1043.     CMPLXtrig0(tmp, tmp);
  1044.     Population = (Rate * tmp.x * tmp.x) - 1.0;
  1045.     return (fabs(Population) > BIG);
  1046.   }
  1047.  
  1048. int LongBifurcStewartTrig()
  1049.   {
  1050. #ifndef XFRACT
  1051.     ltmp.x = lPopulation;
  1052.     ltmp.y = 0;
  1053.     LCMPLXtrig0(ltmp, ltmp);
  1054.     lPopulation = multiply(ltmp.x,ltmp.x,bitshift);
  1055.     lPopulation = multiply(lPopulation,lRate,       bitshift);
  1056.     lPopulation -= fudge;
  1057.     return (overflow);
  1058. #endif
  1059.   }
  1060.  
  1061. int BifurcSetTrigPi()
  1062.   {
  1063.     tmp.x = Population * PI;
  1064.     tmp.y = 0;
  1065.     CMPLXtrig0(tmp, tmp);
  1066.     Population = Rate * tmp.x;
  1067.     return (fabs(Population) > BIG);
  1068.   }
  1069.  
  1070. int LongBifurcSetTrigPi()
  1071.   {
  1072. #ifndef XFRACT
  1073.     ltmp.x = multiply(lPopulation,LPI,bitshift);
  1074.     ltmp.y = 0;
  1075.     LCMPLXtrig0(ltmp, ltmp);
  1076.     lPopulation = multiply(lRate,ltmp.x,bitshift);
  1077.     return (overflow);
  1078. #endif
  1079.   }
  1080.  
  1081. int BifurcAddTrigPi()
  1082.   {
  1083.     tmp.x = Population * PI;
  1084.     tmp.y = 0;
  1085.     CMPLXtrig0(tmp, tmp);
  1086.     Population += Rate * tmp.x;
  1087.     return (fabs(Population) > BIG);
  1088.   }
  1089.  
  1090. int LongBifurcAddTrigPi()
  1091.   {
  1092. #ifndef XFRACT
  1093.     ltmp.x = multiply(lPopulation,LPI,bitshift);
  1094.     ltmp.y = 0;
  1095.     LCMPLXtrig0(ltmp, ltmp);
  1096.     lPopulation += multiply(lRate,ltmp.x,bitshift);
  1097.     return (overflow);
  1098. #endif
  1099.   }
  1100.  
  1101. int BifurcLambdaTrig()
  1102.   {
  1103. /*  Population = Rate * fn(Population) * (1 - fn(Population)) */
  1104.     tmp.x = Population;
  1105.     tmp.y = 0;
  1106.     CMPLXtrig0(tmp, tmp);
  1107.     Population = Rate * tmp.x * (1 - tmp.x);
  1108.     return (fabs(Population) > BIG);
  1109.   }
  1110.  
  1111. int LongBifurcLambdaTrig()
  1112.   {
  1113. #ifndef XFRACT
  1114.     ltmp.x = lPopulation;
  1115.     ltmp.y = 0;
  1116.     LCMPLXtrig0(ltmp, ltmp);
  1117.     ltmp.y = ltmp.x - multiply(ltmp.x,ltmp.x,bitshift);
  1118.     lPopulation = multiply(lRate,ltmp.y,bitshift);
  1119.     return (overflow);
  1120. #endif
  1121.   }
  1122.  
  1123. #define LCMPLXpwr(arg1,arg2,out)    Arg2->l = (arg1); Arg1->l = (arg2);\
  1124.      lStkPwr(); Arg1++; Arg2++; (out) = Arg2->l
  1125.  
  1126. long beta;
  1127.  
  1128. int BifurcMay()
  1129.   { /* X = (lambda * X) / (1 + X)^beta, from R.May as described in Pickover,
  1130.             Computers, Pattern, Chaos, and Beauty, page 153 */
  1131.     tmp.x = 1.0 + Population;
  1132.     tmp.x = pow(tmp.x, -beta); /* pow in math.h included with mpmath.h */
  1133.     Population = (Rate * Population) * tmp.x;
  1134.     return (fabs(Population) > BIG);
  1135.   }
  1136.  
  1137. int LongBifurcMay()
  1138.   {
  1139. #ifndef XFRACT
  1140.     ltmp.x = lPopulation + fudge;
  1141.     ltmp.y = 0;
  1142.     lparm2.x = beta * fudge;
  1143.     LCMPLXpwr(ltmp, lparm2, ltmp);
  1144.     lPopulation = multiply(lRate,lPopulation,bitshift);
  1145.     lPopulation = divide(lPopulation,ltmp.x,bitshift);
  1146.     return (overflow);
  1147. #endif
  1148.   }
  1149.  
  1150. int BifurcMaySetup()
  1151.   {
  1152.  
  1153.    beta = (long)param[2];
  1154.    if(beta < 2)
  1155.       beta = 2;
  1156.    param[2] = (double)beta;
  1157.  
  1158.    timer(0,curfractalspecific->calctype);
  1159.    return(0);
  1160.   }
  1161.  
  1162. /* Here Endeth the Generalised Bifurcation Fractal Engine   */
  1163.  
  1164. /* END Phil Wilson's Code (modified slightly by Kev Allen et. al. !) */
  1165.  
  1166.  
  1167. /******************* standalone engine for "popcorn" ********************/
  1168.  
  1169. int popcorn()    /* subset of std engine */
  1170. {
  1171.    int start_row;
  1172.    start_row = 0;
  1173.    if (resuming)
  1174.    {
  1175.       start_resume();
  1176.       get_resume(sizeof(int),&start_row,0);
  1177.       end_resume();
  1178.    }
  1179.    kbdcount=max_kbdcount;
  1180.    plot = noplot;
  1181.    tempsqrx = ltempsqrx = 0; /* PB added this to cover weird BAILOUTs */
  1182.    for (row = start_row; row <= iystop; row++)
  1183.    {
  1184.       reset_periodicity = 1;
  1185.       for (col = 0; col <= ixstop; col++)
  1186.       {
  1187.      if (StandardFractal() == -1) /* interrupted */
  1188.      {
  1189.         alloc_resume(10,1);
  1190.         put_resume(sizeof(int),&row,0);
  1191.         return(-1);
  1192.      }
  1193.      reset_periodicity = 0;
  1194.       }
  1195.    }
  1196.    calc_status = 4;
  1197.    return(0);
  1198. }
  1199.  
  1200. /******************* standalone engine for "lyapunov" *********************/
  1201. /*** Roy Murphy [76376,721]                        ***/
  1202. /*** revision history:                            ***/
  1203. /*** initial version: Winter '91                    ***/
  1204. /***    Fall '92 integration of Nicholas Wilt's ASM speedups        ***/
  1205. /***    Jan 93' integration with calcfrac() yielding boundary tracing,    ***/
  1206. /***    tesseral, and solid guessing, and inversion, inside=nnn        ***/
  1207. /*** save_release behavior:                        ***/
  1208. /***    1730 & prior: ignores inside=, calcmode='1', (a,b)->(x,y)    ***/
  1209. /***    1731: other calcmodes and inside=nnn                ***/
  1210. /***    1732: the infamous axis swap: (b,a)->(x,y),            ***/
  1211. /***        the order parameter becomes a long int            ***/
  1212. /**************************************************************************/
  1213. int lyaLength, lyaSeedOK;
  1214. int lyaRxy[34];
  1215.  
  1216. #define WES 1   /* define WES to be 0 to use Nick's lyapunov.obj */
  1217. #if WES
  1218. extern int lyapunov_cycles(double, double);
  1219. #else
  1220. extern int lyapunov_cycles(int, double, double, double);
  1221. #endif
  1222.  
  1223. int lyapunov_cycles_in_c(int, double, double);
  1224.  
  1225. lyapunov () {
  1226.     double a, b;
  1227.  
  1228.     if (check_key()) {
  1229.     return -1;
  1230.         }
  1231.     overflow=FALSE;
  1232.     if (param[1]==1) Population = (1.0+rand())/(2.0+RAND_MAX);
  1233.     else if (param[1]==0) {
  1234.     if (fabs(Population)>BIG || Population==0 || Population==1)
  1235.         Population = (1.0+rand())/(2.0+RAND_MAX);
  1236.     }
  1237.     else Population = param[1];
  1238.     (*plot)(col, row, 1);
  1239.     if (invert) {
  1240.     invertz2(&init);
  1241.     a = init.y;
  1242.     b = init.x;
  1243.     }
  1244.     else {
  1245.     a = dy0[row]+dy1[col];
  1246.     b = dx0[col]+dx1[row];
  1247.     }
  1248. #ifndef XFRACT
  1249.     /*  the assembler routines don't work for a & b outside the
  1250.     ranges 0 < a < 4 and 0 < b < 4. So, fall back on the C
  1251.     routines if part of the image sticks out.
  1252.     */
  1253. #if WES
  1254.         color=lyapunov_cycles(a, b);
  1255. #else
  1256.     if (lyaSeedOK && a>0 && b>0 && a<=4 && b<=4)
  1257.     color=lyapunov_cycles(filter_cycles, Population, a, b);
  1258.     else
  1259.     color=lyapunov_cycles_in_c(filter_cycles, a, b);
  1260. #endif
  1261. #else
  1262.     color=lyapunov_cycles_in_c(filter_cycles, a, b);
  1263. #endif
  1264.     if (inside>0 && color==0)
  1265.     color = inside;
  1266.     else if (color>=colors)
  1267.     color = colors-1;
  1268.     (*plot)(col, row, color);
  1269.     return color;
  1270. }
  1271.  
  1272.  
  1273. lya_setup () {
  1274.     /* This routine sets up the sequence for forcing the Rate parameter
  1275.     to vary between the two values.  It fills the array lyaRxy[] and
  1276.     sets lyaLength to the length of the sequence.
  1277.  
  1278.     The sequence is coded in the bit pattern in an integer.
  1279.     Briefly, the sequence starts with an A the leading zero bits
  1280.     are ignored and the remaining bit sequence is decoded.  The
  1281.     sequence ends with a B.  Not all possible sequences can be
  1282.     represented in this manner, but every possible sequence is
  1283.     either represented as itself, as a rotation of one of the
  1284.     representable sequences, or as the inverse of a representable
  1285.     sequence (swapping 0s and 1s in the array.)  Sequences that
  1286.     are the rotation and/or inverses of another sequence will generate
  1287.     the same lyapunov exponents.
  1288.  
  1289.     A few examples follow:
  1290.         number    sequence
  1291.           0    ab
  1292.           1    aab
  1293.           2    aabb
  1294.           3    aaab
  1295.           4    aabbb
  1296.           5    aabab
  1297.           6    aaabb (this is a duplicate of 4, a rotated inverse)
  1298.           7    aaaab
  1299.           8    aabbbb    etc.
  1300.      */
  1301.  
  1302.     long i;
  1303.     int t;
  1304.  
  1305.     if ((filter_cycles=param[2])==0)
  1306.     filter_cycles=maxit/2;
  1307.     lyaSeedOK = param[1]>0 && param[1]<=1 && debugflag!=90;
  1308.     lyaLength = 1;
  1309.  
  1310.     i = param[0];
  1311. #ifndef XFRACT
  1312.     if (save_release<1732) i &= 0x0FFFFL; /* make it a short to reporduce prior stuff*/
  1313. #endif
  1314.     lyaRxy[0] = 1;
  1315.     for (t=31; t>=0; t--)
  1316.     if (i & (1<<t)) break;
  1317.     for (; t>=0; t--)
  1318.     lyaRxy[lyaLength++] = (i & (1<<t)) != 0;
  1319.     lyaRxy[lyaLength++] = 0;
  1320.     if (save_release<1732)        /* swap axes prior to 1732 */
  1321.     for (t=lyaLength; t>=0; t--)
  1322.         lyaRxy[t] = !lyaRxy[t];
  1323.     if (save_release<1731) {        /* ignore inside=, stdcalcmode */
  1324.         stdcalcmode='1';
  1325.     if (inside==1) inside = 0;
  1326.     }
  1327.     if (inside<0) {
  1328.         static char far msg[]=
  1329.         {"Sorry, inside options other than inside=nnn are not supported by the lyapunov"};
  1330.         stopmsg(0,(char far *)msg);
  1331.         inside=1;
  1332.     }
  1333.     return 1;
  1334. }
  1335.  
  1336. int lyapunov_cycles_in_c(int filter_cycles, double a, double b) {
  1337.     int color, count, i, lnadjust;
  1338.     double lyap, total, temp;
  1339.     /* e10=22026.4657948  e-10=0.0000453999297625 */
  1340.  
  1341.     for (i=0; i<filter_cycles; i++) {
  1342.     for (count=0; count<lyaLength; count++) {
  1343.         Rate = lyaRxy[count] ? a : b;
  1344.         if (curfractalspecific->orbitcalc()) {
  1345.         overflow = TRUE;
  1346.         goto jumpout;
  1347.         }
  1348.         }
  1349.     }
  1350.     total = 1.0;
  1351.     lnadjust = 0;
  1352.     for (i=0; i < maxit/2; i++) {
  1353.     for (count = 0; count < lyaLength; count++) {
  1354.         Rate = lyaRxy[count] ? a : b;
  1355.         if (curfractalspecific->orbitcalc()) {
  1356.         overflow = TRUE;
  1357.         goto jumpout;
  1358.         }
  1359.         temp = fabs(Rate-2.0*Rate*Population);
  1360.         if ((total *= (temp))==0) {
  1361.         overflow = TRUE;
  1362.         goto jumpout;
  1363.         }
  1364.         }
  1365.     while (total > 22026.4657948) {
  1366.         total *= 0.0000453999297625;
  1367.         lnadjust += 10;
  1368.         }
  1369.     while (total < 0.0000453999297625) {
  1370.         total *= 22026.4657948;
  1371.         lnadjust -= 10;
  1372.         }
  1373.     }
  1374.  
  1375. jumpout:
  1376.     if (overflow || total <= 0 || (temp = log(total) + lnadjust) > 0)
  1377.     color = 0;
  1378.     else {
  1379.     if (LogFlag)
  1380.     lyap = -temp/((double) lyaLength*i);
  1381.     else
  1382.     lyap = 1 - exp(temp/((double) lyaLength*i));
  1383.     color = 1 + (int)(lyap * (colors-1));
  1384.     }
  1385.     return color;
  1386. }
  1387.  
  1388.  
  1389. /******************* standalone engine for "cellular" ********************/
  1390. /* Originally coded by Ken Shirriff.
  1391.    Modified beyond recognition by Jonathan Osuch.
  1392.      Original or'd the neighborhood, changed to sum the neighborhood
  1393.      Changed prompts and error messages
  1394.      Added CA types
  1395.      Set the palette to some standard? CA colors
  1396.      Changed *cell_array to near and used dstack so put_line and get_line
  1397.        could be used all the time
  1398.      Made space bar generate next screen
  1399.      Increased string/rule size to 16 digits and added CA types 9/20/92
  1400. */
  1401.  
  1402. #define BAD_T         1
  1403. #define BAD_MEM       2
  1404. #define STRING1       3
  1405. #define STRING2       4
  1406. #define TABLEK        5
  1407. #define TYPEKR        6
  1408. #define RULELENGTH    7
  1409.  
  1410. #define CELLULAR_DONE 10
  1411.  
  1412. #ifndef XFRACT
  1413. static BYTE *cell_array[2];
  1414. #else
  1415. static BYTE far *cell_array[2];
  1416. #endif
  1417.  
  1418. S16 r, k_1, rule_digits;
  1419. int lstscreenflag;
  1420.  
  1421. void abort_cellular(int err, int t)
  1422. {
  1423.    int i;
  1424.    switch (err)
  1425.    {
  1426.       case BAD_T:
  1427.          {
  1428.          char msg[30];
  1429.          sprintf(msg,"Bad t=%d, aborting\n", t);
  1430.          stopmsg(0,(char far *)msg);
  1431.          }
  1432.          break;
  1433.       case BAD_MEM:
  1434.          {
  1435.          static char far msg[]={"Insufficient free memory for calculation" };
  1436.          stopmsg(0,msg);
  1437.          }
  1438.          break;
  1439.       case STRING1:
  1440.          {
  1441.          static char far msg[]={"String can be a maximum of 16 digits" };
  1442.          stopmsg(0,msg);
  1443.          }
  1444.          break;
  1445.       case STRING2:
  1446.          {
  1447.          static char far msg[]={"Make string of 0's through  's" };
  1448.          msg[27] = k_1 + 48; /* turn into a character value */
  1449.          stopmsg(0,msg);
  1450.          }
  1451.          break;
  1452.       case TABLEK:
  1453.          {
  1454.          static char far msg[]={"Make Rule with 0's through  's" };
  1455.          msg[27] = k_1 + 48; /* turn into a character value */
  1456.          stopmsg(0,msg);
  1457.          }
  1458.          break;
  1459.       case TYPEKR:
  1460.          {
  1461.          static char far msg[]={"Type must be 21, 31, 41, 51, 61, 22, 32, \
  1462. 42, 23, 33, 24, 25, 26, 27" };
  1463.          stopmsg(0,msg);
  1464.          }
  1465.          break;
  1466.       case RULELENGTH:
  1467.          {
  1468.          static char far msg[]={"Rule must be    digits long" };
  1469.          i = rule_digits / 10;
  1470.          if(i==0)
  1471.             msg[14] = rule_digits + 48;
  1472.          else {
  1473.             msg[13] = i;
  1474.             msg[14] = rule_digits % 10 + 48;
  1475.          }
  1476.          stopmsg(0,msg);
  1477.          }
  1478.          break;
  1479.       case CELLULAR_DONE:
  1480.          break;
  1481.    }
  1482.    if(cell_array[0] != NULL)
  1483. #ifndef XFRACT
  1484.       cell_array[0] = NULL;
  1485. #else
  1486.       farmemfree((char far *)cell_array[0]);
  1487. #endif
  1488.    if(cell_array[1] != NULL)
  1489. #ifndef XFRACT
  1490.       cell_array[1] = NULL;
  1491. #else
  1492.       farmemfree((char far *)cell_array[1]);
  1493. #endif
  1494. }
  1495.  
  1496. cellular () {
  1497.    S16 start_row;
  1498.    S16 filled, notfilled;
  1499.    U16 cell_table[32];
  1500.    U16 init_string[16];
  1501.    U16 kr, k;
  1502.    U32 lnnmbr;
  1503.    U16 i,j,l;
  1504.    S16 t, t2;
  1505.    S32 randparam;
  1506.    double n;
  1507.    char buf[30];
  1508.  
  1509.    set_Cellular_palette();
  1510.  
  1511.    randparam = param[0];
  1512.    lnnmbr = param[3];
  1513.    kr = param[2];
  1514.    switch (kr) {
  1515.      case 21:
  1516.      case 31:
  1517.      case 41:
  1518.      case 51:
  1519.      case 61:
  1520.      case 22:
  1521.      case 32:
  1522.      case 42:
  1523.      case 23:
  1524.      case 33:
  1525.      case 24:
  1526.      case 25:
  1527.      case 26:
  1528.      case 27:
  1529.         break;
  1530.      default:
  1531.         abort_cellular(TYPEKR, 0);
  1532.         return -1;
  1533.         break;
  1534.    }
  1535.  
  1536.    r = kr % 10; /* Number of nearest neighbors to sum */
  1537.    k = kr / 10; /* Number of different states, k=3 has states 0,1,2 */
  1538.    k_1 = k - 1; /* Highest state value, k=3 has highest state value of 2 */
  1539.    rule_digits = (r * 2 + 1) * k_1 + 1; /* Number of digits in the rule */
  1540.  
  1541.    if ((!rflag) && randparam == -1)
  1542.        --rseed;
  1543.    if (randparam != 0 && randparam != -1) {
  1544.       n = param[0];
  1545.       sprintf(buf,"%.16g",n); /* # of digits in initial string */
  1546.       t = strlen(buf);
  1547.       if (t>16 || t <= 0) {
  1548.          abort_cellular(STRING1, 0);
  1549.          return -1;
  1550.       }
  1551.       for (i=0;i<16;i++)
  1552.          init_string[i] = 0; /* zero the array */
  1553.       t2 = (S16) ((16 - t)/2);
  1554.       for (i=0;i<t;i++) { /* center initial string in array */
  1555.          init_string[i+t2] = buf[i] - 48; /* change character to number */
  1556.          if (init_string[i+t2]>k_1) {
  1557.             abort_cellular(STRING2, 0);
  1558.             return -1;
  1559.          }
  1560.       }
  1561.    }
  1562.  
  1563.    srand(rseed);
  1564.    if (!rflag) ++rseed;
  1565.  
  1566. /* generate rule table from parameter 1 */
  1567. #ifndef XFRACT
  1568.    n = param[1];
  1569. #else
  1570.    /* gcc can't manage to convert a big double to an unsigned long properly. */
  1571.    if (param[1]>0x7fffffff) {
  1572.        n = (param[1]-0x7fffffff);
  1573.        n += 0x7fffffff;
  1574.    } else {
  1575.        n = param[1];
  1576.    }
  1577. #endif
  1578.    if (n == 0) { /* calculate a random rule */
  1579.       n = rand15()%k;
  1580.       for (i=1;i<rule_digits;i++) {
  1581.          n *= 10;
  1582.          n += rand15()%k;
  1583.       }
  1584.       param[1] = n;
  1585.    }
  1586.    sprintf(buf,"%.*g",rule_digits ,n);
  1587.    t = strlen(buf);
  1588.    if (rule_digits < t || t < 0) { /* leading 0s could make t smaller */
  1589.       abort_cellular(RULELENGTH, 0);
  1590.       return -1;
  1591.    }
  1592.    for (i=0;i<rule_digits;i++) /* zero the table */
  1593.       cell_table[i] = 0;
  1594.    for (i=0;i<t;i++) { /* reverse order */
  1595.       cell_table[i] = buf[t-i-1] - 48; /* change character to number */
  1596.       if (cell_table[i]>k_1) {
  1597.          abort_cellular(TABLEK, 0);
  1598.          return -1;
  1599.       }
  1600.    }
  1601.  
  1602.  
  1603.    start_row = 0;
  1604. #ifndef XFRACT
  1605.   /* two 4096 byte arrays, at present at most 2024 + 1 bytes should be */
  1606.   /* needed in each array (max screen width + 1) */
  1607.    cell_array[0] = (BYTE *)&dstack[0]; /* dstack is in general.asm */
  1608.    cell_array[1] = (BYTE *)&suffix[0]; /* suffix is in general.asm */
  1609. #else
  1610.    cell_array[0] = (BYTE far *)farmemalloc(ixstop+1);
  1611.    cell_array[1] = (BYTE far *)farmemalloc(ixstop+1);
  1612. #endif
  1613.    if (cell_array[0]==NULL || cell_array[1]==NULL) {
  1614.       abort_cellular(BAD_MEM, 0);
  1615.       return(-1);
  1616.    }
  1617.  
  1618. /* nxtscreenflag toggled by space bar in fractint.c, 1 for continuous */
  1619. /* 0 to stop on next screen */
  1620.  
  1621.    filled = 0;
  1622.    notfilled = 1-filled;
  1623.    if (resuming && !nxtscreenflag && !lstscreenflag) {
  1624.       start_resume();
  1625.       get_resume(sizeof(int),&start_row,0);
  1626.       end_resume();
  1627.       get_line(start_row,0,ixstop,cell_array[filled]);
  1628.    }
  1629.    else if (nxtscreenflag && !lstscreenflag) {
  1630.       start_resume();
  1631.       end_resume();
  1632.       get_line(iystop,0,ixstop,cell_array[filled]);
  1633.       param[3] += iystop + 1;
  1634.       start_row = -1; /* after 1st iteration its = 0 */
  1635.    }
  1636.    else {
  1637.     if(rflag || randparam==0 || randparam==-1){
  1638.       for (col=0;col<=ixstop;col++) {
  1639.          cell_array[filled][col] = rand15()%k;
  1640.       }
  1641.     } /* end of if random */
  1642.  
  1643.     else {
  1644.       for (col=0;col<=ixstop;col++) { /* Clear from end to end */
  1645.          cell_array[filled][col] = 0;
  1646.       }
  1647.       i = 0;
  1648.       for (col=(ixstop-16)/2;col<(ixstop+16)/2;col++) { /* insert initial */
  1649.          cell_array[filled][col] = init_string[i++];    /* string */
  1650.       }
  1651.     } /* end of if not random */
  1652.     if (lnnmbr != 0)
  1653.       lstscreenflag = 1;
  1654.     else
  1655.       lstscreenflag = 0;
  1656.     put_line(start_row,0,ixstop,cell_array[filled]);
  1657.    }
  1658.    start_row++;
  1659.  
  1660. /* This section calculates the starting line when it is not zero */
  1661. /* This section can't be resumed since no screen output is generated */
  1662. /* calculates the (lnnmbr - 1) generation */
  1663.    if (lstscreenflag) { /* line number != 0 & not resuming & not continuing */
  1664.      thinking(1,"Cellular thinking (higher start row takes longer)");
  1665.      for (row = start_row; row < lnnmbr; row++) {
  1666.       if(rflag || randparam==0 || randparam==-1){
  1667.        /* Use a random border */
  1668.        for (i=0;i<=r;i++) {
  1669.          cell_array[notfilled][i]=rand15()%k;
  1670.          cell_array[notfilled][ixstop-i]=rand15()%k;
  1671.        }
  1672.       }
  1673.       else {
  1674.        /* Use a zero border */
  1675.        for (i=0;i<=r;i++) {
  1676.          cell_array[notfilled][i]=0;
  1677.          cell_array[notfilled][ixstop-i]=0;
  1678.        }
  1679.       }
  1680.  
  1681.        for (col=r;col<ixstop-r;col++) {
  1682.          t = 0;
  1683.          for (i=col-r;i<=col+r;i++)
  1684.            t += cell_array[filled][i];
  1685.          if (t>rule_digits || t<0) {
  1686.            thinking(0, NULL);
  1687.            abort_cellular(BAD_T, t);
  1688.            return(-1);
  1689.          }
  1690.          cell_array[notfilled][col] = cell_table[t];
  1691.        }
  1692.        filled = notfilled;
  1693.        notfilled = 1-filled;
  1694.        if (check_key()) {
  1695.           thinking(0, NULL);
  1696.           abort_cellular(CELLULAR_DONE, 0);
  1697.           alloc_resume(10,1);
  1698.           put_resume(sizeof(int),&row,0);
  1699.           return -1;
  1700.        }
  1701.      }
  1702.    start_row = 0;
  1703.    thinking(0, NULL);
  1704.    lstscreenflag = 0;
  1705.    }
  1706.  
  1707. /* This section does all the work */
  1708. contloop:
  1709.    for (row = start_row; row <= iystop; row++) {
  1710.  
  1711.       if(rflag || randparam==0 || randparam==-1){
  1712.        /* Use a random border */
  1713.        for (i=0;i<=r;i++) {
  1714.          cell_array[notfilled][i]=rand15()%k;
  1715.          cell_array[notfilled][ixstop-i]=rand15()%k;
  1716.        }
  1717.       }
  1718.       else {
  1719.        /* Use a zero border */
  1720.        for (i=0;i<=r;i++) {
  1721.          cell_array[notfilled][i]=0;
  1722.          cell_array[notfilled][ixstop-i]=0;
  1723.        }
  1724.       }
  1725.  
  1726.        for (col=r;col<ixstop-r;col++) {
  1727.          t = 0;
  1728.          for (i=col-r;i<=col+r;i++)
  1729.            t += cell_array[filled][i];
  1730.          if (t>rule_digits || t<0) {
  1731.            abort_cellular(BAD_T, t);
  1732.            return(-1);
  1733.          }
  1734.          cell_array[notfilled][col] = cell_table[t];
  1735.        }
  1736.        filled = notfilled;
  1737.        notfilled = 1-filled;
  1738.        put_line(row,0,ixstop,cell_array[filled]);
  1739.        if (check_key()) {
  1740.           abort_cellular(CELLULAR_DONE, 0);
  1741.           alloc_resume(10,1);
  1742.           put_resume(sizeof(int),&row,0);
  1743.           return -1;
  1744.        }
  1745.    }
  1746.    if(nxtscreenflag) {
  1747.      param[3] += iystop + 1;
  1748.      start_row = -1; /* after 1st iteration its = 0 */
  1749.      goto contloop;
  1750.    }
  1751.   abort_cellular(CELLULAR_DONE, 0);
  1752.   return 1;
  1753. }
  1754.  
  1755. int CellularSetup(void)
  1756. {
  1757.    if (!resuming) {
  1758.       nxtscreenflag = 0; /* initialize flag */
  1759.    }
  1760.    timer(0,curfractalspecific->calctype);
  1761.    return(0);
  1762. }
  1763.  
  1764. static void set_Cellular_palette()
  1765. {
  1766.    extern char far *mapdacbox;
  1767.    static Palettetype Red    = {
  1768.       42, 0, 0     };
  1769.    static Palettetype Green  = {
  1770.       10,35,10    };
  1771.    static Palettetype Blue   = {
  1772.       13,12,29    };
  1773.    static Palettetype Yellow = {
  1774.       60,58,18    };
  1775.    static Palettetype Brown = {
  1776.       42,21, 0    };
  1777.    int i;
  1778.  
  1779.    if (mapdacbox) return;        /* map= specified */
  1780.  
  1781.    dacbox[0].red  = 0 ;
  1782.    dacbox[0].green= 0 ;
  1783.    dacbox[0].blue = 0 ;
  1784.  
  1785.    dacbox[1].red     = Red.red;
  1786.    dacbox[1].green = Red.green;
  1787.    dacbox[1].blue  = Red.blue;
  1788.  
  1789.    dacbox[2].red   = Green.red;
  1790.    dacbox[2].green = Green.green;
  1791.    dacbox[2].blue  = Green.blue;
  1792.  
  1793.    dacbox[3].red   = Blue.red;
  1794.    dacbox[3].green = Blue.green;
  1795.    dacbox[3].blue  = Blue.blue;
  1796.  
  1797.    dacbox[4].red   = Yellow.red;
  1798.    dacbox[4].green = Yellow.green;
  1799.    dacbox[4].blue  = Yellow.blue;
  1800.  
  1801.    dacbox[5].red   = Brown.red;
  1802.    dacbox[5].green = Brown.green;
  1803.    dacbox[5].blue  = Brown.blue;
  1804.  
  1805.    SetTgaColors();
  1806.    spindac(0,1);
  1807. }
  1808.  
  1809. /* frothy basin routines */
  1810. static char froth3_256c[] = "froth3.map";
  1811. static char froth6_256c[] = "froth6.map";
  1812. static char froth3_16c[] =  "froth316.map";
  1813. static char froth6_16c[] =  "froth616.map";
  1814. int frothsix=0;
  1815. int froth_altcolor;
  1816. int froth_shades;
  1817. extern int colorstate;
  1818.  
  1819. /* color maps which attempt to replicate the images of James Alexander. */
  1820. static void set_Froth_palette(void)
  1821.    {
  1822.    char *mapname;
  1823.  
  1824.    if (colorstate != 0) /* 0 means dacbox matches default */
  1825.       return;
  1826.  
  1827.    if (colors >= 16)
  1828.       {
  1829.       if (colors >= 256)
  1830.          {
  1831.          if (frothsix)
  1832.             mapname = froth6_256c;
  1833.          else
  1834.             mapname = froth3_256c;
  1835.          }
  1836.       else /* colors >= 16 */
  1837.          {
  1838.          if (frothsix)
  1839.             mapname = froth6_16c;
  1840.          else
  1841.             mapname = froth3_16c;
  1842.          }
  1843.       if (ValidateLuts(mapname) != 0)
  1844.          return;
  1845.       colorstate = 0; /* treat map it as default */
  1846.       spindac(0,1);
  1847.       }
  1848.    }
  1849.  
  1850. int froth_setup(void)
  1851.    {
  1852.    if (param[0] != 3 && param[0] != 6) /* if no match then*/
  1853.       param[0] = 3;                    /* make it 3 */
  1854.    frothsix = param[0] == 6;
  1855.    froth_altcolor = param[1] != 0;
  1856.    froth_shades = (colors-1) / (frothsix ? 6 : 3);
  1857.    /* rqlim needs to be at least 6 or so */
  1858.    if (rqlim < 6.0)
  1859.       rqlim=6.0;
  1860.    set_Froth_palette();
  1861.    /* make the best of the .map situation */
  1862.    orbit_color = !frothsix && colors >= 16 ? (froth_shades<<1)+1 : colors-1;
  1863.    return 1;
  1864.    }
  1865.  
  1866. /* Froth Fractal type */
  1867.   int calcfroth(void)   /* per pixel 1/2/g, called with row & col set */
  1868.      {
  1869.      int found_attractor=0;
  1870.      double x, y, nx, ny, x2, y2;
  1871.      long lx, ly, lnx, lny, lx2, ly2;
  1872.  
  1873. /* These points were determined imperically and verified experimentally */
  1874. /* using the program WL-Plot, a plotting program which has a mode for */
  1875. /* orbits of recursive relations. */
  1876. #define CLOSE    1e-6      /* seems like a good value */
  1877. #define SQRT3    1.732050807568877193
  1878. #define A        1.02871376822
  1879. #define B1       (A/2)
  1880. #define M2       SQRT3
  1881. #define B2       (-A)
  1882. #define M3       (-SQRT3)
  1883. #define B3       (-A)
  1884. #define X1MIN   -1.04368901270
  1885. #define X1MAX    1.33928675524
  1886. #define XMIDT   -0.339286755220
  1887. #define X2MAX1   0.96729063460
  1888. #define XMIDR    0.61508950585
  1889. #define X3MIN1  -0.22419724936
  1890. #define X2MIN2  -1.11508950586
  1891. #define XMIDL   -0.27580275066
  1892. #define X3MAX2   0.07639837810
  1893. #define FROTH_BITSHIFT  28
  1894. /* compiler should handle this at compile time */
  1895. #define D_TO_L(x) ((long)((x)*(1L<<FROTH_BITSHIFT)))
  1896.  
  1897.    orbit_ptr = 0;
  1898.    color = 0;
  1899.    if(showdot>0)
  1900.       (*plot) (col, row, showdot&(colors-1));
  1901.    if (!integerfractal) /* fp mode */
  1902.       {
  1903.       double close=  CLOSE;
  1904.       double a=      A;
  1905.       double b1=     B1;
  1906.       double xmidt=  XMIDT;
  1907.       double m2=     M2;
  1908.       double b2=     B2;
  1909.       double m3=     M3;
  1910.       double b3=     B3;
  1911.       double x1min=  X1MIN;
  1912.       double x1max=  X1MAX;
  1913.       double x2max1= X2MAX1;
  1914.       double xmidr=  XMIDR;
  1915.       double x3min1= X3MIN1;
  1916.       double x2min2= X2MIN2;
  1917.       double xmidl=  XMIDL;
  1918.       double x3max2= X3MAX2;
  1919.  
  1920.       if(invert)
  1921.      {
  1922.      invertz2(&tmp);
  1923.      x = tmp.x;
  1924.      y = tmp.y;
  1925.      }
  1926.       else
  1927.      {
  1928.      x = dx0[col]+dx1[row];
  1929.      y = dy0[row]+dy1[col];
  1930.      }
  1931.  
  1932.       magnitude = (x2=sqr(x)) + (y2=sqr(y));
  1933.       while (!found_attractor && (magnitude < rqlim) && (color < maxit))
  1934.      {
  1935.          /* simple formula: z = z^2 + conj(z*(-1+ai)) */
  1936.      /* but it's the attractor that makes this so interesting */
  1937.      nx = x2 - y2 - x - a*y;
  1938.      ny = (x+x)*y - a*x + y;
  1939.      x = nx;
  1940.      y = ny;
  1941.      if (frothsix) /* repeat mapping */
  1942.         {
  1943.         nx = sqr(x) - sqr(y) - x - a*y;
  1944.         ny = (x+x)*y - a*x + y;
  1945.         x = nx;
  1946.         y = ny;
  1947.         }
  1948.      magnitude = (x2=sqr(x)) + (y2=sqr(y));
  1949.      color++;
  1950.  
  1951.      if (show_orbit)
  1952.         plot_orbit(x, y, -1);
  1953.  
  1954.          if (x > x1min && x < x1max && fabs(b1-y) < close)
  1955.         {
  1956.         if (!frothsix || x < xmidt)
  1957.            found_attractor = 1;
  1958.         else
  1959.            found_attractor = 2;
  1960.  
  1961.         }
  1962.      else if (fabs(m2*x+b2-y) < close)
  1963.         {
  1964.         if (x > xmidr && x < x2max1)
  1965.            found_attractor = !frothsix ? 2 : 4;
  1966.         else if (x > x3min1 && x < xmidr)
  1967.            found_attractor = !frothsix ? 3 : 6;
  1968.         }
  1969.      else if (fabs(m3*x+b3-y) < close)
  1970.         {
  1971.         if (x > x2min2 && x < xmidl)
  1972.            found_attractor = !frothsix ? 2 : 3;
  1973.         else if (x > xmidl && x < x3max2)
  1974.            found_attractor = !frothsix ? 3 : 5;
  1975.         }
  1976.      }
  1977.       }
  1978.    else /* integer mode */
  1979.       {
  1980.       long lclose=   D_TO_L(CLOSE);
  1981.       long la=       D_TO_L(A);
  1982.       long lb1=      D_TO_L(B1);
  1983.       long lxmidt=   D_TO_L(XMIDT);
  1984.       long lm2=      D_TO_L(M2);
  1985.       long lb2=      D_TO_L(B2);
  1986.       long lm3=      D_TO_L(M3);
  1987.       long lb3=      D_TO_L(B3);
  1988.       long lx1min=   D_TO_L(X1MIN);
  1989.       long lx1max=   D_TO_L(X1MAX);
  1990.       long lx2max1=  D_TO_L(X2MAX1);
  1991.       long lxmidr=   D_TO_L(XMIDR);
  1992.       long lx3min1=  D_TO_L(X3MIN1);
  1993.       long lx2min2=  D_TO_L(X2MIN2);
  1994.       long lxmidl=   D_TO_L(XMIDL);
  1995.       long lx3max2=  D_TO_L(X3MAX2);
  1996.  
  1997.       if(invert)
  1998.      {
  1999.      invertz2(&tmp);
  2000.      lx = tmp.x * fudge;
  2001.      ly = tmp.y * fudge;
  2002.      }
  2003.       else
  2004.      {
  2005.      lx = lx0[col] + lx1[row];
  2006.      ly = ly0[row] + ly1[col];
  2007.      }
  2008.  
  2009.       lmagnitud = (lx2=lsqr(lx)) + (ly2=lsqr(ly));
  2010.       while (!found_attractor && (lmagnitud < llimit)
  2011.          && (lmagnitud > 0) && (color < maxit))
  2012.      {
  2013.          /* simple formula: z = z^2 + conj(z*(-1+ai)) */
  2014.      /* but it's the attractor that makes this so interesting */
  2015.      lnx = lx2 - ly2 - lx - multiply(la,ly,bitshift);
  2016.      lny = (multiply(lx,ly,bitshift)<<1) - multiply(la,lx,bitshift) + ly;
  2017.      lx = lnx;
  2018.      ly = lny;
  2019.      if (frothsix)
  2020.         {
  2021.         lmagnitud = (lx2=lsqr(lx)) + (ly2=lsqr(ly));
  2022.         if ((lmagnitud > llimit) || (lmagnitud < 0))
  2023.            break;
  2024.         lnx = lx2 - ly2 - lx - multiply(la,ly,bitshift);
  2025.         lny = (multiply(lx,ly,bitshift)<<1) - multiply(la,lx,bitshift) + ly;
  2026.         lx = lnx;
  2027.         ly = lny;
  2028.         }
  2029.      lmagnitud = (lx2=lsqr(lx)) + (ly2=lsqr(ly));
  2030.      color++;
  2031.  
  2032.      if (show_orbit)
  2033.         iplot_orbit(lx, ly, -1);
  2034.  
  2035.          if (lx > lx1min && lx < lx1max && labs(lb1-ly) < lclose)
  2036.         {
  2037.         if (!frothsix || lx < lxmidt)
  2038.            found_attractor = 1;
  2039.         else
  2040.            found_attractor = 2;
  2041.         }
  2042.      else if (labs(multiply(lm2,lx,bitshift)+lb2-ly) < lclose)
  2043.         {
  2044.         if (lx > lxmidr && lx < lx2max1)
  2045.            found_attractor = !frothsix ? 2 : 4;
  2046.         else if (lx > lx3min1 && lx < lxmidr)
  2047.            found_attractor = !frothsix ? 3 : 6;
  2048.         }
  2049.      else if (labs(multiply(lm3,lx,bitshift)+lb3-ly) < lclose)
  2050.         {
  2051.         if (lx > lx2min2 && lx < lxmidl)
  2052.            found_attractor = !frothsix ? 2 : 3;
  2053.         else if (lx > lxmidl && lx < lx3max2)
  2054.            found_attractor = !frothsix ? 3 : 5;
  2055.         }
  2056.      }
  2057.       }
  2058.    if (show_orbit)
  2059.       scrub_orbit();
  2060.  
  2061.    realcolor = color;
  2062.    if ((kbdcount -= realcolor) <= 0)
  2063.       {
  2064.       if (check_key())
  2065.      return (-1);
  2066.       kbdcount = max_kbdcount;
  2067.       }
  2068.  
  2069. /* inside - Here's where non-palette based images would be nice.  Instead, */
  2070. /* we'll use blocks of (colors-1)/3 or (colors-1)/6 and use special froth  */
  2071. /* color maps in attempt to replicate the images of James Alexander.       */
  2072.    if (found_attractor)
  2073.       {
  2074.       if (colors >= 256)
  2075.      {
  2076.      if (!froth_altcolor)
  2077.         {
  2078.         if (color > froth_shades)
  2079.         color = froth_shades;
  2080.         }
  2081.      else
  2082.         color = froth_shades * color / maxit;
  2083.      if (color == 0)
  2084.         color = 1;
  2085.      color += froth_shades * (found_attractor-1);
  2086.      }
  2087.       else if (colors >= 16)
  2088.      { /* only alternate coloring scheme available for 16 colors */
  2089.      long lshade;
  2090.  
  2091. /* Trying to make a better 16 color distribution. */
  2092. /* Since their are only a few possiblities, just handle each case. */
  2093. /* This is a mostly guess work here. */
  2094.      lshade = ((long)color<<16)/maxit;
  2095.      if (!frothsix)
  2096.         {
  2097.         if (lshade < 2622)       /* 0.04 */
  2098.            color = 1;
  2099.         else if (lshade < 10486) /* 0.16 */
  2100.            color = 2;
  2101.         else if (lshade < 23593) /* 0.36 */
  2102.            color = 3;
  2103.         else if (lshade < 41943) /* 0.64 */
  2104.            color = 4;
  2105.         else
  2106.            color = 5;
  2107.         color += 5 * (found_attractor-1);
  2108.         }
  2109.      else
  2110.         {
  2111.         if (lshade < 10486)      /* 0.16 */
  2112.            color = 1;
  2113.         else
  2114.            color = 2;
  2115.         color += 2 * (found_attractor-1);
  2116.         }
  2117.      }
  2118.       else /* use a color corresponding to the attractor */
  2119.      color = found_attractor;
  2120.       oldcolor = color;
  2121.       }
  2122.    else if (color >= maxit)
  2123.       color = oldcolor; /* inside, but didn't get sucked in by attractor. */
  2124.    else /* outside */
  2125.       color = 0; /* all outside points are color 0 */
  2126.  
  2127.    (*plot)(col, row, color);
  2128.  
  2129.    return color;
  2130.  
  2131. #undef CLOSE
  2132. #undef SQRT3
  2133. #undef A
  2134. #undef B1
  2135. #undef XMIDT
  2136. #undef M2
  2137. #undef B2
  2138. #undef M3
  2139. #undef B3
  2140. #undef X1MIN
  2141. #undef X1MAX
  2142. #undef X2MAX1
  2143. #undef XMIDR
  2144. #undef X3MIN1
  2145. #undef X2MIN2
  2146. #undef XMIDL
  2147. #undef X3MAX2
  2148. #undef FROTH_BITSHIFT
  2149. #undef D_TO_L
  2150.    }
  2151.