home *** CD-ROM | disk | FTP | other *** search
/ Shareware Overload / ShartewareOverload.cdr / graf / fract4.zip / CALCFRAC.C next >
C/C++ Source or Header  |  1989-12-23  |  44KB  |  1,417 lines
  1. /*
  2.  
  3. FRACTALS.C and CALCFRAC.C actually calculate the fractal images (well,
  4. SOMEBODY had to do it!).  The modules are set up so that all logic that
  5. is independent of any fractal-specific code is in CALCFRAC.C, and the
  6. code that IS fractal-specific is in FRACTALS.C. Original author Tim Wegner,
  7. but just about ALL the authors have contributed SOME code to this routine
  8. at one time or another, or contributed to one of the many massive
  9. restructurings.
  10.  
  11.    -------------------------------------------------------------------- */
  12.  
  13. #include <stdio.h>
  14. #include <stdlib.h>
  15. #include <float.h>
  16. #include <math.h>
  17. #include <dos.h>
  18. #include <limits.h>
  19. #include <time.h>
  20. #include <stdarg.h>
  21. #include "fractint.h"
  22. #include "fmath.h"
  23. #include "targa_lc.h"
  24.  
  25. extern int bailout;
  26. extern char far plasmamessage[];
  27. long lmagnitud, llimit, llimit2, lclosenuff,l16triglim;
  28. struct complex init,tmp,old,new,saved;
  29. extern double tempsqrx,tempsqry;
  30. extern long ltempsqrx,ltempsqry;
  31. extern int biomorph;
  32. extern struct lcomplex linit;
  33. extern int basin;
  34.  
  35. int color, oldcolor, oldmax, oldmax10, row, col, passes;
  36. int iterations, invert;
  37. double f_radius,f_xcenter, f_ycenter; /* for inversion */
  38. double far *dx0, far *dy0;
  39. long XXOne, FgOne, FgTwo, LowerLimit;
  40.  
  41. extern int LogFlag;
  42. extern char far *LogTable;
  43.  
  44. void (*plot)() = putcolor;
  45.  
  46. extern double inversion[];        /* inversion radius, f_xcenter, f_ycenter */
  47. extern int    xdots, ydots;        /* coordinates of dots on the screen  */
  48. extern int    colors;                /* maximum colors available */
  49. extern int    inside;                /* "inside" color to use    */
  50. extern int    maxit;                /* try this many iterations */
  51. extern int    fractype;            /* fractal type */
  52. extern int    numpasses;            /* 0 = 1 pass, 1 = double pass */
  53. extern int    solidguessing;            /* 1 if solid-guessing */
  54. extern int    debugflag;            /* for debugging purposes */
  55. extern int    timerflag;            /* ditto */
  56. extern    int    diskvideo;            /* for disk-video klooges */
  57.  
  58. extern int rflag, rseed;
  59. extern int decomp[];
  60.  
  61. extern double    param[];        /* parameters */
  62. extern double    potparam[];        /* potential parameters */
  63. extern long    lx0[], ly0[];        /* X, Y points */
  64. extern long    delx,dely;            /* X, Y increments */
  65. extern long    fudge;                /* fudge factor (2**n) */
  66. extern int    bitshift;            /* bit shift for fudge */
  67. extern char potfile[];          /* potential filename */
  68.  
  69.  
  70. #ifndef sqr
  71. #define sqr(x) ((x)*(x))
  72. #endif
  73.  
  74. extern    int extraseg;
  75.  
  76. /* These are local but I don't want to pass them as parameters */
  77. static    unsigned char top;   /* flag to indicate top of calcfract */ 
  78. static    unsigned char c;
  79. static      int i;
  80.  
  81. double xxmin,xxmax,yymin,yymax; /* corners */
  82. struct complex lambda;
  83. double deltaX, deltaY;
  84. double magnitude, rqlim, rqlim2;
  85. int XXdots, YYdots; /* local dots */
  86. extern struct complex parm;
  87. int (*calctype)();
  88. double closenuff;
  89. int pixelpi; /* value of pi in pixels */
  90. FILE *fp_pot;
  91. int potflag; /* tells if continuous potential on  */
  92. unsigned long lm;        /* magnitude limit (CALCMAND) */
  93.  
  94. /* ORBIT variables */
  95. int    show_orbit;            /* flag to turn on and off */
  96. int    orbit_ptr;            /* pointer into save_orbit array */
  97. static int  far *save_orbit;        /* array to save orbit values */
  98. static int orbit_color=15;        /* XOR color */
  99.  
  100. int    ixstart, ixstop, iystart, iystop;    /* (for CALCMAND) start, stop here */
  101. int    symmetry;            /* symmetry flag for calcmand()    */
  102. int    guessing;            /* solid-guessing flag for calcmand() */
  103.  
  104. static    int    integerfractal;        /* TRUE if fractal uses integer math */
  105.  
  106. /* -------------------------------------------------------------------- */
  107. /*        These variables are external for speed's sake only    */
  108. /* -------------------------------------------------------------------- */
  109.  
  110. extern struct lcomplex lnew,llambda;
  111. int periodicitycheck;
  112.  
  113. extern double floatmin, floatmax;
  114.  
  115. int StandardFractal();
  116.  
  117.  
  118. calcfract()
  119. {
  120.    orbit_color = 15;
  121.    top     = 1;
  122.    if(fp_pot)
  123.    {
  124.       fclose(fp_pot);
  125.       fp_pot = NULL;
  126.    }
  127.  
  128.    if (invert)
  129.    {
  130.       floatmin = FLT_MIN;
  131.       floatmax = FLT_MAX;
  132.    }
  133.    ixstart = 0;                /* default values to disable ..    */
  134.    iystart = 0;                /*   solid guessing        */
  135.    ixstop = xdots-1;
  136.    iystop = ydots-1;
  137.    guessing = 0;
  138.  
  139.    basin = 0;
  140.    
  141.    xxmin  = (double)lx0[      0] / fudge;  
  142.    xxmax  = (double)lx0[xdots-1] / fudge;
  143.    yymax  = (double)ly0[      0] / fudge;
  144.    yymin  = (double)ly0[ydots-1] / fudge;
  145.    deltaX = (double)lx0[      1] / fudge - xxmin;
  146.    deltaY = yymax - (double)ly0[      1] / fudge;
  147.  
  148.    parm.x  = param[0];
  149.    parm.y  = param[1];
  150.  
  151.    if(extraseg)
  152.    {
  153. #ifdef __TURBOC__
  154.       dx0 = MK_FP(extraseg,0);
  155. #else
  156.       FP_SEG(dx0)=extraseg;
  157.       FP_OFF(dx0)=0;
  158. #endif
  159.    }else
  160.       return(-1);
  161.    dy0 = dx0 + MAXPIXELS;
  162.  
  163.    if(fabs(potparam[0]) > 0.0)
  164.       potflag = 1;
  165.    else
  166.       potflag = 0;   
  167.  
  168.    if (LogFlag)
  169.       ChkLogFlag();
  170.  
  171.    lm = 4;                /* CALCMAND magnitude limit */
  172.    lm = lm << bitshift;
  173. /*        Continuous potential override (unused at the moment)
  174.    if (potparam[2] > 4 && potparam[2] < 64)
  175.       lm = potparam[2]*fudge;
  176. */
  177.    /* set modulus bailout value if 3rd potparam set */
  178.  
  179.    if (fabs(potparam[2]) > 0.0) 
  180.       rqlim = potparam[2];
  181.    else if (decomp[0] > 0 && decomp[1] > 0) 
  182.       rqlim = (double)decomp[1];
  183.    else if(bailout)    /* user input bailout */
  184.       rqlim = bailout;
  185.    else if(biomorph != -1)   /* biomorph benefits from larger bailout */
  186.       rqlim = 100;   
  187.    else    
  188.       rqlim = fractalspecific[fractype].orbit_bailout;
  189.  
  190.    /* ORBIT stuff */
  191.    save_orbit = (int far *)(dx0 + 2*MAXPIXELS);
  192.    show_orbit = 0;
  193.    orbit_ptr = 0;
  194.    if(colors < 16)
  195.        orbit_color = 1;
  196.  
  197.    integerfractal = fractalspecific[fractype].isinteger;
  198.    if ((!integerfractal) || invert)
  199.    {
  200.       /* set up dx0 and dy0 analogs of lx0 and ly0 in Bert's "extra" segment */
  201.       /* put fractal parameters in doubles */
  202.       for(i=0;i<xdots;i++)
  203.          dx0[i] = (double)lx0[i] / fudge;
  204.       for(i=0;i<ydots;i++)
  205.         dy0[i] = (double)ly0[i] / fudge;
  206.    }
  207.  
  208.    if (integerfractal)         /* the bailout limit can't be too high here */
  209.       if (rqlim > 127.0) rqlim = 127.0;
  210.    if(inversion[0] != 0.0)
  211.    {
  212.       f_radius    = inversion[0];
  213.       f_xcenter   = inversion[1];
  214.       f_ycenter   = inversion[2];
  215.       
  216.       if (inversion[0] < 0.0)  /*  auto calc radius 1/6 screen */
  217.          inversion[0] = f_radius = (min(xxmax - xxmin,yymax - yymin)) / 6.0;
  218.       
  219.       if (invert < 2) { /* xcenter not already set */
  220.          inversion[1] = f_xcenter = (xxmin + xxmax) / 2.0;
  221.          if (fabs(f_xcenter) < (xxmax-xxmin) / 100)
  222.              inversion[1] = f_xcenter = 0.0;
  223.          }
  224.  
  225.       if (invert < 3) { /* ycenter not already set */
  226.          inversion[2] = f_ycenter = (yymin + yymax) / 2.0;
  227.          if (fabs(f_ycenter) < (yymax-yymin) / 100)
  228.              inversion[2] = f_ycenter = 0.0;
  229.          }
  230.  
  231.       invert = 3; /* so values will not be changed if we come back */
  232.    }
  233.  
  234.    setsymmetry(fractalspecific[fractype].symmetry);
  235.  
  236.    if (potfile[0] != 0)             /* potential file? */
  237.    {
  238.       numpasses = 0;                /* disable dual-pass */
  239.       solidguessing = 0;            /* disable solid-guessing */
  240.    }
  241.    
  242.    closenuff = ((delx < dely ? delx : dely) >> 1); /* for periodicity checking */
  243.    closenuff /= fudge;
  244.    rqlim2 = sqrt(rqlim);
  245.    if (integerfractal)         /* for integer routines (lambda) */
  246.    {
  247.       llambda.x = parm.x * fudge;        /* real portion of Lambda */
  248.       llambda.y = parm.y * fudge;        /* imaginary portion of Lambda */
  249.       llimit = rqlim * fudge;        /* stop if magnitude exceeds this */
  250.       if (llimit <= 0) llimit = 0x7fffffff; /* klooge for integer math */
  251.       llimit2 = rqlim2 * fudge;    /* stop if magnitude exceeds this */
  252.       lclosenuff = closenuff * fudge;   /* "close enough" value */
  253.       l16triglim = 8L<<16;   /* domain limit of fast trig functions */
  254.    }
  255.  
  256.    calctype = fractalspecific[fractype].calctype;        /* assume a standard fractal type */
  257.  
  258.    periodicitycheck = 1;        /* turn on periodicity checking */
  259.    
  260.    /* per_image routine is run here */   
  261.    if (fractalspecific[fractype].per_image())  /* a stand-alone type? */
  262.    {
  263.       /* fractal routine (usually StandardFractal) is run here */
  264.       if(solidguessing)     /* use std solid-guessing? */
  265.          timer(solidguess,0);
  266.       else
  267.          timer(calctype,0);
  268.    }
  269.  
  270.    potfile[0] = NULL;
  271.    if (! check_key())
  272.       return(0);
  273.    else
  274.       return(1);
  275. }
  276.  
  277. test()
  278. {
  279.    teststart();
  280.    for (passes=0; passes <= numpasses ; passes++)
  281.    for (row = passes; row < YYdots; row=row+1+numpasses)
  282.    {
  283.       register int col;
  284.       init.y = dy0[row];
  285.       for (col = 0; col < XXdots; col++)       /* look at each point on screen */
  286.       {
  287.          register color;
  288.          init.x = dx0[col];
  289.          if(check_key()) 
  290.          {
  291.             testend();
  292.             return(-1);
  293.          }
  294.          color = testpt(init.x,init.y,parm.x,parm.y,maxit,inside);
  295.          (*plot)(col,row,color);
  296.          if(numpasses && (passes == 0))
  297.              (*plot)(col,row+1,color);
  298.       }
  299.    }
  300.    testend();
  301.    return(0);
  302. }
  303.  
  304. StandardFractal()
  305. {
  306.    int caught_a_cycle;
  307.    int savedand, savedincr;        /* for periodicity checking */
  308.    int quad;
  309.    struct lcomplex lsaved;
  310.  
  311.    /* note: numpasses must == 1 in solidguessing mode */
  312.    /* This is set in solidguess() */   
  313.    for (passes=0; passes <= numpasses ; passes++)
  314.    {
  315.       /* guessing == 1 means first pass only */
  316.       if(guessing == 1 && passes == 1)
  317.          break;
  318.          
  319.       /* guessing == 2 means second pass only */   
  320.       if(guessing == 2 && passes == 0)
  321.          continue;
  322.          
  323.       /* if numpasses == 1, on first pass skip rows */
  324.       for (row = iystart; row <= iystop; row = row + 1 + numpasses - passes)
  325.       {
  326.          oldcolor = 1;
  327.          oldmax = min(maxit, 250);
  328.          oldmax10 = oldmax - 10;
  329.          
  330.          /* really fractal specific, but we'll leave it here */
  331.          if (! integerfractal) 
  332.             init.y = dy0[row];
  333.          else 
  334.             linit.y = ly0[row];
  335.  
  336.          /* if numpasses == 1, on first pass skip cols */
  337.          for (col = ixstart; col <= ixstop; col = col + 1 + numpasses - passes)
  338.          {
  339.             /* on second pass even row,col pts already done */
  340.             if(passes) /* second pass */
  341.                if(row&1 ==0 && col&1 == 0)
  342.                   continue;
  343.             if(check_key())
  344.                   return(-1);
  345.             if (! integerfractal) 
  346.             {
  347.                saved.x = 0;
  348.                saved.y = 0;
  349.             }
  350.             else 
  351.             {
  352.                lsaved.x = 0;
  353.                lsaved.y = 0;
  354.             }
  355.             orbit_ptr = 0;        
  356.             color = 0;
  357.             caught_a_cycle = 1;
  358.             savedand = 1;        /* begin checking every other cycle */
  359.             savedincr = 1;        /* start checking the very first time */
  360.             
  361.             fractalspecific[fractype].per_pixel(); /* initialize the calculations */
  362.             while (++color < maxit)
  363.             {
  364.                /* calculation of one orbit goes here */
  365.                /* input in "old" -- output in "new" */
  366.  
  367.                if(fractalspecific[fractype].orbitcalc()) break;
  368.                if(show_orbit) 
  369.                {
  370.                   if (! integerfractal)
  371.                      plot_orbit(new.x,new.y,-1);
  372.                   else
  373.                      iplot_orbit(lnew.x,lnew.y,-1);
  374.                }
  375.  
  376.                if (oldcolor >= oldmax10)
  377.                   if (periodicitycheck)        /* only if this is OK to do */ 
  378.                   {
  379.                      if ((color & savedand) == 0) /* time to save a new value */
  380.                      {
  381.                         if (! integerfractal)
  382.                            saved = new;        /* floating pt fractals */
  383.                         else 
  384.                            lsaved = lnew;        /* integer fractals */
  385.                         if (--savedincr == 0) /* time to lengthen the periodicity? */
  386.                         { 
  387.                            savedand = (savedand << 1) + 1; /* longer periodicity */
  388.                            savedincr = 4;                  /* restart counter */
  389.                         }
  390.                      }
  391.                      else /* check against an old save */
  392.                      {
  393.                         if (! integerfractal)    /* floating-pt periodicity chk */
  394.                         { 
  395.                            if((fabs(saved.x-new.x)+fabs(saved.y-new.y)) 
  396.                                  < closenuff)
  397.                            {
  398.                               caught_a_cycle = 7;
  399.                               color = maxit;
  400.                            }  
  401.                         }
  402.                         else  /* integer periodicity check */
  403.                         {
  404.                            if(labs(lsaved.x-lnew.x)+labs(lsaved.y-lnew.y) 
  405.                                  < lclosenuff) 
  406.                            {
  407.                               caught_a_cycle = 7;
  408.                               color = maxit;
  409.                            }
  410.                         }
  411.                      }
  412.                   }
  413.             }
  414.             if (oldcolor < maxit && color >= maxit) 
  415.             {
  416.                oldmax = oldcolor;
  417.                oldmax10 = oldmax - 10;
  418.             }
  419.             if(show_orbit)
  420.                scrub_orbit();
  421.             oldcolor = color;
  422.             if(color < maxit)
  423.             {
  424.                caught_a_cycle = color;   
  425.                if (color == 0)color = 1; /* needed to make same as calcmand */
  426.                if(debugflag==98) color = caught_a_cycle;
  427.             }
  428.  
  429.             if(potflag) 
  430.             {
  431.                if(integerfractal)         /* adjust integer fractals */
  432.                { 
  433.                   new.x = lnew.x; new.x /= fudge;
  434.                   new.y = lnew.y; new.y /= fudge;
  435.                }
  436.                magnitude = sqr(new.x)+sqr(new.y);
  437.  
  438.                /* printf("rqlim %e magnit %e new.x %e new.y %e\n",rqlim,magnitude,new.x,new.y); */
  439.  
  440.                /* MUST call potential every pixel!! */
  441.                color = potential(magnitude,color); 
  442.                if(color < 0 || color >= colors) continue; 
  443.             }
  444.             else if (decomp[0] > 0)
  445.             {
  446.                int temp = 0;
  447.                struct lcomplex lalt;
  448.                struct complex alt;
  449.                color = 0;
  450.    
  451.                if (integerfractal) /* the only case */
  452.                {
  453.                   if (lnew.y < 0)
  454.                   {
  455.                      temp = 2;
  456.                      lnew.y = -lnew.y;
  457.                   }
  458.    
  459.                   if (lnew.x < 0)
  460.                   {
  461.                      ++temp;
  462.                      lnew.x = -lnew.x;
  463.                   }
  464.    
  465.                   if (decomp[0] >= 8)
  466.                   {
  467.                      temp <<= 1;
  468.                      if (lnew.x < lnew.y)
  469.                      {
  470.                         ++temp;
  471.                         lalt.x = lnew.x; /* just */
  472.                         lnew.x = lnew.y; /* swap */
  473.                         lnew.y = lalt.x; /* them */
  474.                      }
  475.    
  476.                      if (decomp[0] >= 16)
  477.                      {
  478.                         temp <<= 1;
  479.                         if (lnew.x * 0.414213562373 < lnew.y) /* tan(22.5) */
  480.                         {
  481.                            ++temp;
  482.                            lalt = lnew;
  483.                            lnew.x = lalt.x * 0.707106781187   /* cos(45) */
  484.                                   + lalt.y * 0.707106781187; /* sin(45) */
  485.                            lnew.y = lalt.x * 0.707106781187   /* sin(45) */
  486.                                   - lalt.y * 0.707106781187; /* cos(45) */
  487.                         }
  488.    
  489.                         if (decomp[0] >= 32)
  490.                         {
  491.                            temp <<= 1;
  492.                            if (lnew.x * 0.19891236738 < lnew.y) /* tan(11.25) */
  493.                            {
  494.                               ++temp;
  495.                               lalt = lnew;
  496.                               lnew.x = lalt.x * 0.923879532511   /* cos(22.5) */
  497.                                      + lalt.y * 0.382683432365;  /* sin(22.5) */
  498.                               lnew.y = lalt.x * 0.382683432365   /* sin(22.5) */
  499.                                      - lalt.y * 0.923879532511;  /* cos(22.5) */
  500.                            }
  501.    
  502.                            if (decomp[0] >= 64)
  503.                            {
  504.                               temp <<= 1;
  505.                               if (lnew.x * 0.0984914033572 < lnew.y)
  506.                               {
  507.                                  ++temp;
  508.                                  lalt = lnew;
  509.                                  lnew.x = lalt.x * 0.980785280403
  510.                                         + lalt.y * 0.195090322016;
  511.                                  lnew.y = lalt.x * 0.195090322016
  512.                                         - lalt.y * 0.980785280403;
  513.                               }
  514.    
  515.                               if (decomp[0] >= 128)
  516.                               {
  517.                                  temp <<= 1;
  518.                                  if (lnew.x * 0.0491268497695 < lnew.y)
  519.                                  {
  520.                                     ++temp;
  521.                                     lalt = lnew;
  522.                                     lnew.x = lalt.x * 0.998795456205
  523.                                            + lalt.y * 0.0490676743274;
  524.                                     lnew.y = lalt.x * 0.0490676743274
  525.                                            - lalt.y * 0.998795456205;
  526.                                  }
  527.    
  528.                                  if (decomp[0] == 256)
  529.                                  {
  530.                                     temp <<= 1;
  531.                                     if ((lnew.x * 0.0245486221089 < lnew.y))
  532.                                         ++temp;
  533.                                  }
  534.                               }
  535.                            }
  536.                         }
  537.                      }
  538.                   }
  539.                }
  540.                else /* double case */
  541.                {
  542.                   if (new.y < 0)
  543.                   {
  544.                      temp = 2;
  545.                      new.y = -new.y;
  546.                   }
  547.    
  548.                   if (new.x < 0)
  549.                   {
  550.                      ++temp;
  551.                      new.x = -new.x;
  552.                   }
  553.    
  554.                   if (decomp[0] >= 8)
  555.                   {
  556.                      temp <<= 1;
  557.                      if (new.x < new.y)
  558.                      {
  559.                         ++temp;
  560.                         alt.x = new.x; /* just */
  561.                         new.x = new.y; /* swap */
  562.                         new.y = alt.x; /* them */
  563.                      }
  564.    
  565.                      if (decomp[0] >= 16)
  566.                      {
  567.                         temp <<= 1;
  568.                         if (new.x * 0.414213562373 < new.y) /* tan(22.5) */
  569.                         {
  570.                            ++temp;
  571.                            alt = new;
  572.                            new.x = alt.x * 0.707106781187   /* cos(45) */
  573.                                  + alt.y * 0.707106781187;  /* sin(45) */
  574.                            new.y = alt.x * 0.707106781187   /* sin(45) */
  575.                                  - alt.y * 0.707106781187;  /* cos(45) */
  576.                         }
  577.    
  578.                         if (decomp[0] >= 32)
  579.                         {
  580.                            temp <<= 1;
  581.                            if (new.x * 0.19891236738 < new.y) /* tan(11.25) */
  582.                            {
  583.                               ++temp;
  584.                               alt = new;
  585.                               new.x = alt.x * 0.923879532511   /* cos(22.5) */
  586.                                     + alt.y * 0.382683432365;  /* sin(22.5) */
  587.                               new.y = alt.x * 0.382683432365   /* sin(22.5) */
  588.                                     - alt.y * 0.923879532511;  /* cos(22.5) */
  589.                            }
  590.    
  591.                            if (decomp[0] >= 64)
  592.                            {
  593.                               temp <<= 1;
  594.                               if (new.x * 0.0984914033572 < new.y)
  595.                               {
  596.                                  ++temp;
  597.                                  alt = new;
  598.                                  new.x = alt.x * 0.980785280403
  599.                                        + alt.y * 0.195090322016;
  600.                                  new.y = alt.x * 0.195090322016
  601.                                        - alt.y * 0.980785280403;
  602.                               }
  603.    
  604.                               if (decomp[0] >= 128)
  605.                               {
  606.                                  temp <<= 1;
  607.                                  if (new.x * 0.0491268497695 < new.y)
  608.                                  {
  609.                                     ++temp;
  610.                                     alt = new;
  611.                                     new.x = alt.x * 0.998795456205
  612.                                           + alt.y * 0.0490676743274;
  613.                                     new.y = alt.x * 0.0490676743274
  614.                                           - alt.y * 0.998795456205;
  615.                                  }
  616.    
  617.                                  if (decomp[0] == 256)
  618.                                  {
  619.                                     temp <<= 1;
  620.                                     if ((new.x * 0.0245486221089 < new.y))
  621.                                         ++temp;
  622.                                  }
  623.                               }
  624.                            }
  625.                         }
  626.                      }
  627.                   }
  628.                }
  629.                for (i = 1; temp > 0; ++i)
  630.                {
  631.                   if (temp & 1)
  632.                      color = (1 << i) - 1 - color;
  633.                   temp >>= 1;
  634.                }
  635.                if (decomp[0] == 2)        color &= 1;
  636.                if (colors > decomp[0])    color++;
  637.             }
  638.             else if(biomorph != -1)
  639.             {
  640.                if(integerfractal)
  641.                {
  642.                   if(labs(lnew.x) < llimit2 || labs(lnew.y) < llimit2)
  643.                      color = biomorph;
  644.                }
  645.                else if(fabs(new.x) < rqlim2 || fabs(new.y) < rqlim2)
  646.                   color = biomorph;
  647.             }
  648.             
  649.             if(oldcolor >= maxit && inside >= 0) /* really color, not oldcolor */
  650.                color = inside;
  651.             if(LogFlag)
  652.                color = LogTable[color];
  653.  
  654.             (*plot)(col,row,color);
  655.             if(numpasses && (passes == 0))
  656.             {
  657.                 (*plot)(col  ,row+1,color);
  658.                 (*plot)(col+1,row  ,color);
  659.                 (*plot)(col+1,row+1,color);
  660.             }
  661.          }
  662.       }
  663.    }
  664.    return(0);
  665. }
  666.  
  667. #if 0
  668. MainNewton()
  669. {
  670.    Newton(); /* Lee's Newton */
  671.    return(0);
  672. }
  673. #endif
  674.  
  675. /* Symmetry plot for period PI */
  676. void symPIplot( x, y, color)
  677. int x, y, color ;
  678. {
  679.   while(x < xdots)
  680.   {
  681.      putcolor( x, y, color) ;
  682.      x += pixelpi;
  683.   }   
  684. }
  685. /* Symmetry plot for period PI plus Origin Symmetry */
  686. void symPIplot2J( x, y, color)
  687. int x, y, color ;
  688. {
  689.   while(x < xdots)
  690.   {
  691.      putcolor( x, y, color) ;
  692.      putcolor( xdots - x - 1, (ydots - y - 1), color) ;
  693.      x += pixelpi;
  694.   }   
  695. }
  696. /* Symmetry plot for period PI plus Both Axis Symmetry */
  697. void symPIplot4J( x, y, color)
  698. int x, y, color ;
  699. {
  700.   while(x < xdots>>1)
  701.   {
  702.      putcolor( x            , y            , color) ;
  703.      putcolor( x            , ydots - y - 1, color) ;
  704.      putcolor( xdots - x - 1, y            , color) ;
  705.      putcolor( xdots - x - 1, ydots - y - 1, color) ;
  706.      x += pixelpi;
  707.   }   
  708. }
  709.  
  710. /* Symmetry plot for X Axis Symmetry */
  711. void symplot2( x, y, color)
  712. int x, y, color ;
  713. {
  714.   putcolor( x, y, color) ;
  715.   putcolor( x, (ydots - y - 1), color) ;
  716. }
  717.  
  718. /* Symmetry plot for Y Axis Symmetry */
  719. void symplot2Y( x, y, color)
  720. int x, y, color ;
  721. {
  722.   putcolor( x, y, color) ;
  723.   putcolor( (xdots - x - 1), y, color) ;
  724. }
  725.  
  726. /* Symmetry plot for Origin Symmetry */
  727. void symplot2J( x, y, color)
  728. int x, y, color ;
  729. {
  730.   putcolor( x, y, color) ;
  731.   putcolor( xdots - x - 1, (ydots - y - 1), color) ;
  732. }
  733.  
  734. /* Symmetry plot for Both Axis Symmetry */
  735. void symplot4( x, y, color)
  736. int x, y, color ;
  737. {
  738.   putcolor( x            , y            , color) ;
  739.   putcolor( x            , ydots - y - 1, color) ;
  740.   putcolor( xdots - x - 1, y            , color) ;
  741.   putcolor( xdots - x - 1, ydots - y - 1, color) ;
  742. }
  743.  
  744. /* Symmetry plot for X Axis Symmetry - Newtbasin version */
  745. void symplot2basin( x, y, color)
  746. int x, y, color ;
  747. {
  748.   extern int degree;
  749.   putcolor( x, y, color) ;
  750.   putcolor( x, (ydots - y - 1),( degree+1-color)%degree+1)  ;
  751. }
  752.  
  753. /* Symmetry plot for Both Axis Symmetry  - Newtbasin version */
  754. void symplot4basin( x, y, color)
  755. int x, y, color ;
  756. {
  757.   extern int degree;
  758.   int color1;
  759.   if(color == 0) /* assumed to be "inside" color */
  760.   {
  761.      symplot4( x, y, color);
  762.      return;
  763.   }
  764.   color1 = degree/2+degree+2 - color;
  765.   if(color < degree/2+2)
  766.      color1 = degree/2+2 - color;
  767.   else
  768.      color1 = degree/2+degree+2 - color;
  769.   putcolor( x            , y            , color) ;
  770.   putcolor( x            , ydots - y - 1,  (degree+1 - color)%degree+1) ;
  771.   putcolor( xdots - x - 1, y            ,  color1) ;
  772.   putcolor( xdots - x - 1, ydots - y - 1, (degree+1 - color1)%degree+1) ;
  773. }
  774.  
  775. /* Do nothing plot!!! */
  776. void noplot(int x,int y,int color)
  777. {
  778. }
  779.  
  780. static int iparmx;                    /* iparmx = parm.x * 16 */
  781. static int shiftvalue;                /* shift based on #colors */
  782.  
  783. typedef struct palett
  784. {
  785.    unsigned char red;
  786.    unsigned char green;
  787.    unsigned char blue;
  788. } Palettetype;
  789.  
  790. extern int colors;
  791. extern int xdots,ydots;
  792. extern Palettetype dacbox[256];
  793. extern int cpu, daccount, cyclelimit;
  794. static int plasma_check;            /* to limit kbd checking */
  795.  
  796. void adjust(int xa,int ya,int x,int y,int xb,int yb)
  797. {
  798.    long pseudorandom;
  799.    if(getcolor(x,y) != 0)
  800.       return;
  801.  
  802.    pseudorandom = ((long)iparmx)*((rand()-16383));
  803.    pseudorandom = pseudorandom*(abs(xa-xb)+abs(ya-yb));
  804.    pseudorandom = pseudorandom >> shiftvalue;
  805.    pseudorandom = ((getcolor(xa,ya)+getcolor(xb,yb)+1)>>1)+pseudorandom;
  806.    if (pseudorandom <   1) pseudorandom =   1;
  807.    if (pseudorandom >= colors) pseudorandom = colors-1;
  808.    putcolor(x,y,(int)pseudorandom);
  809. }
  810.  
  811. void subDivide(int x1,int y1,int x2,int y2)
  812. {
  813.    int x,y;
  814.    int v;
  815.    if ((++plasma_check & 0x7f) == 1)
  816.       if(check_key()) 
  817.       {
  818.          plasma_check--;
  819.          return;
  820.       }
  821.    if(x2-x1<2 && y2-y1<2)
  822.       return;
  823.    x = (x1+x2)>>1;
  824.    y = (y1+y2)>>1;
  825.    adjust(x1,y1,x ,y1,x2,y1);
  826.    adjust(x2,y1,x2,y ,x2,y2);
  827.    adjust(x1,y2,x ,y2,x2,y2);
  828.    adjust(x1,y1,x1,y ,x1,y2);
  829.    if(getcolor(x,y) == 0)
  830.    {
  831.       v = (getcolor(x1,y1)+getcolor(x2,y1)+getcolor(x2,y2)+
  832.           getcolor(x1,y2)+2)>>2;
  833.       putcolor(x,y,v);
  834.    }
  835.    subDivide(x1,y1,x ,y );
  836.    subDivide(x ,y1,x2,y );
  837.    subDivide(x ,y ,x2,y2);
  838.    subDivide(x1,y ,x ,y2);
  839. }
  840.  
  841. plasma()
  842. {
  843.    
  844.    if(colors < 4 || diskvideo) {
  845.     setvideomode(3,0,0,0);
  846.     buzzer(2);
  847.     helpmessage(plasmamessage);
  848.     return(-1);
  849.     }
  850.    iparmx = param[0] * 8;
  851.    if (parm.x <= 0.0) iparmx = 16;
  852.    if (parm.x >= 100) iparmx = 800;
  853.  
  854.    if (!rflag) rseed = (int)time(NULL);
  855.       srand(rseed);
  856.  
  857.    if (colors == 256)            /* set the (256-color) palette */
  858.       set_Plasma_palette();        /* skip this if < 256 colors */
  859.  
  860.    if (colors > 16)            /* adjust based on # of colors */
  861.       shiftvalue = 18;
  862.    else {
  863.       if (colors > 4)
  864.          shiftvalue = 22;
  865.       else {
  866.          if (colors > 2)
  867.             shiftvalue = 24;
  868.          else
  869.             shiftvalue = 25;
  870.          }
  871.       }
  872.  
  873.    putcolor(      0,      0,1+(((rand()/colors)*(colors-1))>>(shiftvalue-11)));
  874.    putcolor(xdots-1,      0,1+(((rand()/colors)*(colors-1))>>(shiftvalue-11)));
  875.    putcolor(xdots-1,ydots-1,1+(((rand()/colors)*(colors-1))>>(shiftvalue-11)));
  876.    putcolor(      0,ydots-1,1+(((rand()/colors)*(colors-1))>>(shiftvalue-11)));
  877.  
  878.    subDivide(0,0,xdots-1,ydots-1);
  879.    if (! check_key())
  880.     return(0);
  881.    else
  882.     return(1);
  883. }
  884.  
  885. set_Plasma_palette()
  886. {
  887.    static Palettetype Red    = {63, 0, 0};
  888.    static Palettetype Green  = { 0,63, 0};
  889.    static Palettetype Blue   = { 0, 0,63};
  890.    int i;
  891.  
  892.      if( CustomLut() ) return;        /* TARGA 3 June 89 j mclain */
  893.  
  894.    dacbox[0].red  = 0 ;
  895.    dacbox[0].green= 0 ;
  896.    dacbox[0].blue = 0 ;
  897.    for(i=1;i<=85;i++)
  898.    {
  899.        dacbox[i].red   = (i*Green.red   + (86-i)*Blue.red  )/85;
  900.        dacbox[i].green = (i*Green.green + (86-i)*Blue.green)/85;
  901.        dacbox[i].blue  = (i*Green.blue  + (86-i)*Blue.blue )/85;
  902.  
  903.        dacbox[i+85].red   = (i*Red.red   + (86-i)*Green.red  )/85;
  904.        dacbox[i+85].green = (i*Red.green + (86-i)*Green.green)/85;
  905.        dacbox[i+85].blue  = (i*Red.blue  + (86-i)*Green.blue )/85;
  906.  
  907.        dacbox[i+170].red   = (i*Blue.red   + (86-i)*Red.red  )/85;
  908.        dacbox[i+170].green = (i*Blue.green + (86-i)*Red.green)/85;
  909.        dacbox[i+170].blue  = (i*Blue.blue  + (86-i)*Red.blue )/85;
  910.    }
  911.     ValidateLuts( NULL );    /* TARGA 3 June 89  j mclain */
  912.      spindac(0,1);
  913. }
  914.  
  915. check_key()
  916. {
  917. int key;
  918.   if((key = keypressed()) != 0) {
  919.     if(key != 'o' && key != 'O') 
  920.        return(-1);
  921.    getakey();
  922.    show_orbit = 1 - show_orbit;
  923.    }
  924.    return(0);
  925. }
  926.  
  927. iplot_orbit(ix, iy, color)
  928. long ix, iy;
  929. int color;
  930. {
  931. int i, j;
  932.   if (ix < lx0[0] || ix > lx0[xdots-1] || iy < ly0[ydots-1] || iy > ly0[0]
  933.     || diskvideo)
  934.     return(0);
  935.   i = (ix - lx0[0]) / delx;
  936.   j = ydots - (iy - ly0[ydots-1]) / dely;
  937.  
  938.   /* save orbit value */
  939.   if( 0 <= i && i < xdots && 0 <= j && j < ydots && orbit_ptr < 1000)
  940.   {
  941.      if(color == -1)
  942.      {
  943.         *(save_orbit + orbit_ptr++) = i; 
  944.         *(save_orbit + orbit_ptr++) = j; 
  945.         putcolor(i,j,getcolor(i,j)^orbit_color);
  946.      }
  947.      else
  948.         putcolor(i,j,color);
  949.   }
  950. }
  951.  
  952. plot_orbit(real,imag,color)
  953. double real,imag;
  954. int color;
  955. {
  956.   long ix,iy;
  957.   if((real < xxmin) || (real > xxmax) || (imag < yymin) || (imag > yymax))
  958.      return(0);
  959.   /* convert to longs */
  960.   ix = real * fudge;
  961.   iy = imag * fudge;
  962.   iplot_orbit(ix,iy,color);
  963. }  
  964. scrub_orbit()
  965. {
  966.    int i,j;
  967.    int last_orbit; 
  968.    last_orbit = orbit_ptr;
  969.    orbit_ptr = 0;
  970.    while(orbit_ptr < last_orbit)
  971.    {
  972.      i = *(save_orbit + orbit_ptr++);
  973.      j = *(save_orbit + orbit_ptr++); 
  974.      putcolor(i,j,getcolor(i,j)^orbit_color);
  975.    }
  976. }
  977. /* timer function. Assumes #include <time.h> */
  978. timer(int (*fractal)(),int argc,...)
  979. {
  980.    va_list arg_marker;    /* variable arg list */
  981.    int argcount;    /* how many optional args */
  982.    char *timestring;
  983.    time_t ltime;
  984.    FILE *fp;
  985.    int out;
  986.    long t1,t2;
  987.    int args[4];
  988.    int i;
  989.  
  990.    /* argc = number of optional arguments */   
  991.    va_start(arg_marker,argc);
  992.  
  993.    i=0;
  994.    while(i < argc)
  995.       args[i++] = va_arg(arg_marker,int);
  996.          
  997.    if(timerflag)
  998.       fp=fopen("bench","a");
  999.    t1 = clock();
  1000.    switch(argc)
  1001.    {
  1002.    case 0:
  1003.       out = (*fractal)();
  1004.       break;
  1005.    case 1:
  1006.       out = (*fractal)(args[0]);
  1007.       break;
  1008.    case 2:
  1009.       out = (*fractal)(args[0],args[1]);
  1010.       break;
  1011.    case 3:
  1012.       out = (*fractal)(args[0],args[1],args[2]);
  1013.       break;
  1014.    case 4:
  1015.       out = (*fractal)(args[0],args[1],args[2],args[3]);
  1016.       break;
  1017.    default:
  1018.       return(-1);
  1019.       break;
  1020.    }
  1021.    t2 = clock();
  1022.    if(timerflag)
  1023.    {
  1024.       time(<ime);
  1025.       timestring = ctime(<ime);
  1026.       timestring[24] = 0; /*clobber newline in time string */
  1027.  
  1028.       /* at the moment only decode timer takes a parameter */
  1029.       if(argc)
  1030.          fprintf(fp,"decode ");
  1031.       fprintf(fp,"%s type=%s resolution = %dx%d maxiter=%d",
  1032.            timestring,
  1033.            fractalspecific[fractype].name,
  1034.            xdots,
  1035.            ydots,
  1036.            maxit);
  1037.       fprintf(fp," time=%8.4f secs\n",
  1038.            ((float)(t2-t1))/CLK_TCK);
  1039.       if(fp != NULL)
  1040.           fclose(fp);
  1041.    }
  1042.    return(out); 
  1043. }
  1044.  
  1045. /* these variables are used to tune algorithm */
  1046. int block_size;      /* multiple of grid used for block */
  1047.                      /* smaller blocks catch more blocks but also have */
  1048.                      /* more overhead. */     
  1049. int xblock,yblock;   /* pixel dimension of solid blocks */
  1050. int xres,yres;       /* pixel dimension of solid-guessing grid */  
  1051.  
  1052. solidguess()
  1053. {
  1054.    int oldnumpasses;
  1055.    int i,j;
  1056.  
  1057.    /* save numpasses - this logic effectively bypasses it */ 
  1058.    oldnumpasses = numpasses; 
  1059.    numpasses = 1;
  1060.  
  1061.    /* initial tuning values */
  1062.    block_size = 1;      /* look for block size == grid size */
  1063.  
  1064.    /* set up grid and block */
  1065.    xres = 2;
  1066.    yres = 2;
  1067.  
  1068.    xblock = xres*block_size;
  1069.    yblock = yres*block_size;
  1070.  
  1071.    guessing = 1;            /* solid-guessing - first pass */
  1072.  
  1073.    (*calctype)();            /* generate the overall grid */
  1074.      if (check_key()) {        /* TARGA 2 June 89 j mclain */
  1075.         numpasses = oldnumpasses;
  1076.          return(-1);
  1077.       }                                                     
  1078.  
  1079.    guessing=2;              /* solid guessing - second pass */
  1080.  
  1081.    for(j=0;j<YYdots;j += yblock)
  1082.    {
  1083.       iystart = j;  
  1084. /*      iystop  = min(j+yblock-1,YYdots-1);  */
  1085.       iystop  = j+yblock-1;  
  1086.       for(i=0;i<XXdots;i += xblock)
  1087.       {
  1088.          ixstart = i;
  1089. /*         ixstop  = min(i+xblock-1,XXdots-1); */
  1090.          ixstop  = i+xblock-1;
  1091.          if(checkblock())
  1092.          {
  1093.             (*calctype)();
  1094.             if (check_key()) 
  1095.             {
  1096.                numpasses = oldnumpasses;
  1097.                return(-1);
  1098.             }
  1099.          }
  1100.       }
  1101.    }
  1102.    numpasses = oldnumpasses;
  1103.    guessing  = 0;
  1104. }
  1105.  
  1106. checkblock()
  1107. {
  1108.    /* check all the grid points in this block */
  1109.    /* to see if the same color */
  1110.    int same;
  1111.    int thecolor;
  1112.    int i,j; 
  1113.    same=1;
  1114.    for(j=iystart;j<=iystop+1;j += yres) {
  1115.       if (j >= ydots) continue;
  1116.       for(i=ixstart;i<=ixstop+1;i += xres)
  1117.       {
  1118.          if (i >= xdots) continue;
  1119.          if(j == iystart && i == ixstart)
  1120.             thecolor=getcolor(i,j);
  1121.          else
  1122.          {
  1123.             if(getcolor(i,j) != thecolor)
  1124.                return(1); /* darn! all colors not the same! */
  1125.          }   
  1126.       }
  1127.    }
  1128.    return(0); /* yup, all colors the same! Guess its a solid block! */
  1129. }      
  1130.  
  1131. /******************************************************************/
  1132. /* Continuous potential calculation for Mandelbrot and Julia      */
  1133. /* Reference: Science of Fractal Images p. 190.                   */
  1134. /* Special thanks to Mark Peterson for his "MtMand" program that  */
  1135. /* beautifully approximates plate 25 (same reference) and spurred */
  1136. /* on the inclusion of similar capabilities in FRACTINT.          */
  1137. /*                                                                */
  1138. /* The purpose of this function is to calculate a color value     */
  1139. /* for a fractal that varies continuously with the screen pixels  */
  1140. /* locations for better rendering in 3D.                          */
  1141. /*                                                                */
  1142. /* Here "magnitude" is the modulus of the orbit value at          */
  1143. /* "iterations". The potparms[] are user-entered paramters        */
  1144. /* controlling the level and slope of the continuous potential    */
  1145. /* surface. Returns color.  - Tim Wegner 6/25/89                  */
  1146. /*                                                                */
  1147. /*                     -- Change history --                       */
  1148. /*                                                                */
  1149. /* 09/12/89   - added floatflag support and fixed float underflow */
  1150. /*                                                                */
  1151. /******************************************************************/
  1152.  
  1153. int potential(double mag, int iterations)
  1154. {
  1155.    extern char potfile[];     /* name of pot save file */
  1156.    extern struct fractal_info save_info;    /*  for saving data */
  1157.  
  1158.    static int x,y;            /* keep track of position in image */
  1159.    float f_mag,f_tmp,pot;
  1160.    double d_tmp;
  1161.    int i_pot;
  1162.    unsigned int intbuf;
  1163.    FILE *t16_create();
  1164.    
  1165.    extern int maxit;          /* maximum iteration bailout limit */
  1166.    extern double potparam[];     /* continuous potential parameters */
  1167.    extern char potfile[];     /* pot savename */
  1168.    extern unsigned int decoderline[];
  1169.    extern char floatflag;
  1170.  
  1171.    if(top) 
  1172.    {
  1173.       top = 0;
  1174.       x   = 0;
  1175.       y   = 0;
  1176.       /* make sure file is closed */
  1177.       if(fp_pot)
  1178.          fclose(fp_pot);
  1179.       if(potfile[0])   
  1180.          fp_pot = t16_create(potfile, xdots, ydots, 
  1181.              sizeof(struct fractal_info), &save_info);
  1182.    }
  1183.    if(iterations < maxit)
  1184.    {
  1185.       i_pot = pot = iterations+2;
  1186.       if(i_pot <= 0 || mag <= 1.0)
  1187.          pot = 0.0;
  1188.       else /* pot = log(mag) / pow(2.0, (double)pot); */ 
  1189.       {
  1190.          i_pot = pot;
  1191.          if(i_pot < 120 && !floatflag) /* empirically determined limit of fShift */
  1192.          {
  1193.             f_mag = mag;
  1194.             fLog14(f_mag,f_tmp); /* this SHOULD be non-negative */
  1195.             fShift(f_tmp,-i_pot,pot);
  1196.          }   
  1197.          else
  1198.          {
  1199.             d_tmp = log(mag)/(double)pow(2.0,(double)pot);
  1200.             if(d_tmp > FLT_MIN) /* prevent float type underflow */
  1201.                pot = d_tmp;
  1202.             else
  1203.                pot = 0.0;
  1204.          }
  1205.       }
  1206.       /* following transformation strictly for aesthetic reasons */
  1207.       /* meaning of parameters: 
  1208.       potparam[0] -- zero potential level - highest color - 
  1209.       potparam[1] -- slope multiplier -- higher is steeper
  1210.       potparam[2] -- rqlim value if changeable (bailout for modulus) */
  1211.  
  1212.       if(pot > 0.0)
  1213.       {
  1214.          if(floatflag)
  1215.             pot = (float)sqrt((double)pot);
  1216.          else
  1217.          {    
  1218.             fSqrt14(pot,f_tmp);
  1219.             pot = f_tmp; 
  1220.          }   
  1221.          pot = potparam[0] - pot*potparam[1] - 1.0;
  1222.       }
  1223.       else
  1224.          pot = potparam[0] - 1.0;
  1225.    }
  1226.    else if(inside > 0)
  1227.       pot = inside; /* negative inside turns off inside replacement */
  1228.    else
  1229.       pot = potparam[0];
  1230.  
  1231.    if(pot > colors)
  1232.       pot = colors -1;
  1233.    if(pot < 1.0)
  1234.       pot = 1.0; /* avoid color 0 */     
  1235.     
  1236.     intbuf = pot*(1<<8); /* shift 8 bits */
  1237.     decoderline[x] = intbuf;
  1238.     if(++x == xdots)
  1239.     {
  1240.        /* new row */
  1241.        if(fp_pot)
  1242.           t16_putline(fp_pot, xdots, decoderline);
  1243.        x = 0;
  1244.        if(++y == ydots)
  1245.        {
  1246.           /* we're done */
  1247.           t16_flush(fp_pot);
  1248.           fclose(fp_pot);
  1249.           fp_pot = NULL;
  1250.           potfile[0] = NULL;
  1251.        }
  1252.     }
  1253.     return((int)pot);
  1254. }      
  1255.  
  1256. /*
  1257.    "intpotential" is called by the MANDEL/JULIA routines with scaled long
  1258.    integer magnitudes.  The routine just calls "potential".  Bert
  1259. */
  1260.  
  1261. int intpotential(unsigned long longmagnitude, int iterations)
  1262. {                 /* this version called from calcmand() */
  1263.    double magnititude;
  1264.  
  1265.    magnitude = ((double)longmagnitude) / fudge;
  1266.    magnitude = 256*magnitude;
  1267.    return(potential(magnitude, iterations));
  1268. }
  1269. setsymmetry(int sym)    /* set up proper symmetrical plot functions */
  1270. {
  1271.    extern int forcesymmetry;
  1272.    plot = putcolor;     /* assume no symmetry */
  1273.    XXdots = xdots;
  1274.    YYdots = ydots;
  1275.    iystop = YYdots -1;
  1276.    ixstop = XXdots -1;
  1277.    symmetry = 1;
  1278.    if(sym == NOPLOT && forcesymmetry == 999)
  1279.    {
  1280.       plot = noplot;
  1281.       return;
  1282.    }   
  1283.    /* NOTE: either diskvideo or 16-bit potential disables symmetry */
  1284.    /* also ant decomp= option and any inversion not about the origin */
  1285.    if ((!potfile[0]) && (!diskvideo) && (!invert || inversion[2] == 0.0)
  1286.    && decomp[0] == 0) 
  1287.    {
  1288.       if(sym != XAXIS && sym != XAXIS_NOPARM && inversion[1] != 0.0 && forcesymmetry == 999)
  1289.          return;
  1290.       else if(forcesymmetry != 999)
  1291.          sym = forcesymmetry;
  1292.       switch(sym)     /* symmetry switch */
  1293.       {
  1294.       case XAXIS_NOPARM:            /* X-axis Symmetry  (no params)*/
  1295.          if (parm.x != 0.0 || parm.y != 0.0)
  1296.             break;
  1297.       case XAXIS:            /* X-axis Symmetry */
  1298.          if (fabs(yymax + yymin) >= fabs(yymax - yymin)*.01 )
  1299.             if(forcesymmetry==999)
  1300.                break;
  1301.          symmetry = 0;
  1302.          if(basin)
  1303.             plot = symplot2basin;
  1304.          else   
  1305.             plot = symplot2;
  1306.          YYdots /= 2;    
  1307.          iystop = YYdots -1;
  1308.          break;
  1309.       case YAXIS_NOPARM:            /* Y-axis Symmetry (No Parms)*/
  1310.          if (parm.x != 0.0 || parm.y != 0.0)
  1311.             break;
  1312.       case YAXIS:            /* Y-axis Symmetry */
  1313.          if (fabs(xxmax + xxmin) >= fabs(xxmax - xxmin)*.01)
  1314.             if(forcesymmetry==999)
  1315.                break;
  1316.          plot = symplot2Y;
  1317.          XXdots /= 2;
  1318.          ixstop = XXdots -1;
  1319.          break;
  1320.       case XYAXIS_NOPARM:            /* X-axis AND Y-axis Symmetry (no parms)*/
  1321.          if(parm.x != 0.0 || parm.y != 0.0)
  1322.             break;
  1323.       case XYAXIS:            /* X-axis AND Y-axis Symmetry */
  1324.          if(forcesymmetry!=999)
  1325.          {
  1326.             if(basin)
  1327.                plot = symplot4basin;
  1328.             else   
  1329.                plot = symplot4 ;
  1330.             XXdots /= 2 ;
  1331.             YYdots /= 2 ;
  1332.             iystop = YYdots -1;
  1333.             ixstop = XXdots -1;
  1334.             break;
  1335.          }
  1336.          if ((fabs(yymax + yymin) < fabs(yymax - yymin)*.01 &&
  1337.              fabs(xxmax + xxmin) < fabs(xxmax - xxmin)*.01))
  1338.          {
  1339.             if(basin)
  1340.                plot = symplot4basin;
  1341.             else   
  1342.                plot = symplot4 ;
  1343.             XXdots /= 2 ;
  1344.             YYdots /= 2 ;
  1345.             iystop = YYdots -1;
  1346.             ixstop = XXdots -1;
  1347.          }
  1348.          else if (fabs(yymax + yymin) < fabs(yymax - yymin)*.01)
  1349.          {
  1350.             if(basin)
  1351.                plot = symplot2basin;
  1352.             else   
  1353.                plot = symplot2 ;
  1354.             YYdots /= 2 ;
  1355.             iystop = YYdots -1;
  1356.          }
  1357.          break;
  1358.       case ORIGIN_NOPARM:            /* Origin Symmetry (no parms)*/
  1359.          if (parm.x != 0.0 || parm.y != 0.0)
  1360.             break;
  1361.       case ORIGIN:            /* Origin Symmetry */
  1362.          if (fabs(yymax + yymin) >= fabs(yymax - yymin)*.01 ||
  1363.              fabs(xxmax + xxmin) >= fabs(xxmax - xxmin)*.01)
  1364.             if(forcesymmetry==999)
  1365.                break;
  1366.          symmetry = 0;
  1367.          plot = symplot2J;
  1368.          YYdots /= 2 ;
  1369.          iystop = YYdots -1;
  1370.          break;
  1371.       case PI_SYM_NOPARM:
  1372.          if (parm.x != 0.0 || parm.y != 0.0)
  1373.             break;
  1374.       case PI_SYM:            /* PI symmetry */
  1375.          if(invert && forcesymmetry == 999) 
  1376.          {
  1377.             if (fabs(yymax + yymin) >= fabs(yymax - yymin)*.01 ||
  1378.                 fabs(xxmax + xxmin) >= fabs(xxmax - xxmin)*.01)
  1379.                break;
  1380.             symmetry = 0;
  1381.             plot = symplot2J;
  1382.             YYdots /= 2 ;
  1383.             iystop = YYdots -1;
  1384.             break;
  1385.          }
  1386.          pixelpi = (PI/fabs(xxmax-xxmin))*xdots; /* PI in pixels */
  1387.          if(pixelpi < xdots)
  1388.          {
  1389.             XXdots = pixelpi;
  1390.             ixstop = XXdots-1;
  1391.          }
  1392.    
  1393.          if (fabs(yymax + yymin) < fabs(yymax - yymin)*.01 &&
  1394.              fabs(xxmax + xxmin) < fabs(xxmax - xxmin)*.01)
  1395.          {
  1396.             if(parm.y == 0.0)
  1397.             {
  1398.                plot = symPIplot4J ;
  1399.                YYdots /= 2 ;
  1400.                iystop = YYdots -1;
  1401.             }
  1402.             else
  1403.             {
  1404.                plot = symPIplot2J;
  1405.                YYdots /= 2 ;
  1406.                iystop = YYdots -1;
  1407.             }
  1408.          }                
  1409.          else if(!invert || forcesymmetry != 999)
  1410.             plot = symPIplot ;
  1411.          break;
  1412.       default:            /* no symmetry */
  1413.          break;
  1414.       }
  1415.    }
  1416. }
  1417.