/*
   This is a simple but rather interesting program that graphs equations.
   You enter equations in terms of x and y, something like one of the
   following:
  
    y = sin[x]
  
    x^2 + y^2 = 81
  
    y cos[x] = x sin[y]
  
   This version of the program can also graph INEQUALITIES, which have
   less-than or greater-than symbols instead of just equals.
  
   For example, try
  
   abs[y^2 + x^4 - 1] < cos[x]
   
   This uses a recursive method to subdivide and test rectangles.
  
   This version of the program adds automatic generation of grid lines and
   axis labels using the Grid.frink sample library.
*/
   
use Grid.frink

class Plot2D
{
   // Boundaries of the plot.  These can be changed by calling setBounds.
   var xmin = -10
   var xmax =  10
   var ymin = -10
   var ymax =  10

   // Change the doublings to vary the number of pixels.  This is the number
   // of doublings, so if the number is 10 we have 2^10=1024 doublings for
   // a resolution of 1024x1024.  (That is over a million pixels!  Don't
   // be surprised that graphing at that resolution takes a long time!)
   // Be warned that increasing the doublings by 1 makes 4 times as many pixels!
   var doublings = 10

   /** Plots the function (specified as a string) and returns a VoxelArray. */
   plot[func] :=
   {
      hasInequality = false
      certEq = undef
      lasteq = certFunc = func

      g = new graphics

      // If there's an inequality, let's make a test equation to see if we can
      // plot an entire rectangle using the "CERTAINLY" comparators.
      if func =~ %r/([<>]|!=)/
      {
         hasInequality = true
         g.antialiased[false]
         certFunc =~ %s/<=/ CLE /g  // Replace <= with certainly less than or equals
         certFunc =~ %s/>=/ CGE /g  // Replace >= with certainly greater than or equals
         certFunc =~ %s/</ CLT /g   // Replace <  with certainly less than
         certFunc =~ %s/>/ CGT /g   // Replace >  with certainly greater than
         certFunc =~ %s/!=/ CNE /g   // Replace = with certainly not equals
         certFunc =~ %s/=/ CEQ /g   // Replace = with certainly equals
         certEq = parseToExpression[certFunc]
      }

      // These replacements turn normal comparator and equality tests into
      // "POSSIBLY EQUALS" tests.
      func =~ %s/<=/ PLE /g  // Replace <= with possibly less than or equals
      func =~ %s/>=/ PGE /g  // Replace >= with possibly greater than or equals
      func =~ %s/</ PLT /g   // Replace <  with possibly less than
      func =~ %s/>/ PGT /g   // Replace >  with possibly greater than
      func =~ %s/!=/ PNE /g   // Replace = with possibly not equals
      func =~ %s/=/ PEQ /g   // Replace = with possibly equals
      eq = parseToExpression[func]

      // Change the last number to vary the resolution.  This is the number
      // of doublings, so if the number is 10 we have 2^10=1024 doublings for
      // a resolution of 1024x1024.
      testRect[xmin, xmax, ymin, ymax, g, eq, certEq, doublings]

      grid = new Grid
      grid.auto[g]
      g.add[grid.getGrid[]]

      return g
   }

   /* Recursive function to test an interval containing the specified bounds.
      If no possible solution exists, the recursion halts.  If a possible
      solution exists, this breaks it down into 4 sub-rectangles and tests each
      of them recursively.  level is the maximum number of levels to split, so
      the total resolution of the final graph will be 2^level.
   */
   testRect[x1, x2, y1, y2, g, eq, certEq, level] :=
   {
      nextLevel = level - 1
      x = new interval[x1, x2]
      y = new interval[y1, y2]
      
      // Test the rectangle.  If it possibly contains solutions, recursively
      // subdivide.
      res = eval[eq]
      
      if res or res==undef
      {
         if (nextLevel >= 0)
         {
            // Do we have inequalities and a CERTAINLY test?
            if (certEq != undef)
               certRes = eval[certEq]

            if certRes == true
            {
               // If the entire rectangle is a solution, fill the rectangle
               // and stop further recursion on this rectangle.
               g.fillRectSides[x1, -y1, x2, -y2]
               return
            }

            // Further subdivide the rectangle into 4 quadrants and
            // recursively test them all
            cx = (x1 + x2)/2
            cy = (y1 + y2)/2
            testRect[x1, cx, y1, cy, g, eq, certEq, nextLevel]
            testRect[cx, x2, y1, cy, g, eq, certEq, nextLevel]
            testRect[x1, cx, cy, y2, g, eq, certEq, nextLevel]
            testRect[cx, x2, cy, y2, g, eq, certEq, nextLevel]
         } else
         if (res)             // Valid point
            g.fillRectSides[x1, -y1, x2, -y2]
         else
         {
            // Error in evaluating point, plot in red.
            g.color[1,0,0]
            g.fillRectSides[x1, -y1, x2, -y2]
            g.color[0,0,0]
         }
      }            
   }

   /** Shows an already-rendered graphics on the screen. */
   show[g] :=
   {
      g.show[]
   }

   /** Sets the bounds to be plotted. */
   setBounds[xmin, xmax, ymin, ymax] :=
   {
      this.xmin = xmin
      this.xmax = xmax
      this.ymin = ymin
      this.ymax = ymax
   }

   /** Sets the number of doublings.  The number of voxels rendered will be
       2^doublings on each side. */
   setDoublings[d] :=
   {
      doublings = d
   }
}

"Plot2D included okay."