Initial population.

1 files changed, 341 insertions, 0 deletions

diff --git a/src/algebra/clip.spad.pamphlet b/src/algebra/clip.spad.pamphlet
new file mode 100644
index 00000000..04f9f42a
--- /dev/null
+++ b/src/algebra/clip.spad.pamphlet
@@ -0,0 +1,341 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra clip.spad}
+\author{Clifton J. Williamson}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{package CLIP TwoDimensionalPlotClipping}
+<<package CLIP TwoDimensionalPlotClipping>>=
+)abbrev package CLIP TwoDimensionalPlotClipping
+++ Automatic clipping for 2-dimensional plots
+++ Author: Clifton J. Williamson
+++ Date Created: 22 December 1989
+++ Date Last Updated: 10 July 1990
+++ Keywords: plot, singularity
+++ Examples:
+++ References:
+ 
+TwoDimensionalPlotClipping(): Exports == Implementation where
+  ++ The purpose of this package is to provide reasonable plots of
+  ++ functions with singularities.
+  B      ==> Boolean
+  L      ==> List
+  SEG    ==> Segment
+  RN     ==> Fraction Integer
+  SF     ==> DoubleFloat
+  Pt     ==> Point DoubleFloat
+  PLOT   ==> Plot
+  CLIPPED ==> Record(brans: L L Pt,xValues: SEG SF,yValues: SEG SF)
+ 
+  Exports ==> with
+    clip: PLOT -> CLIPPED
+      ++ clip(p) performs two-dimensional clipping on a plot, p, from
+      ++ the domain \spadtype{Plot} for the graph of one variable,
+      ++ \spad{y = f(x)}; the default parameters \spad{1/4} for the fraction 
+      ++ and \spad{5/1} for the scale are used in the \spadfun{clip} function.
+    clip: (PLOT,RN,RN) -> CLIPPED
+      ++ clip(p,frac,sc) performs two-dimensional clipping on a plot, p, 
+      ++ from the domain \spadtype{Plot} for the graph of one variable
+      ++ \spad{y = f(x)}; the fraction parameter is specified by \spad{frac} 
+      ++ and the scale parameter is specified by \spad{sc} for use in the 
+      ++ \spadfun{clip} function.
+    clipParametric: PLOT -> CLIPPED
+      ++ clipParametric(p) performs two-dimensional clipping on a plot,
+      ++ p, from the domain \spadtype{Plot} for the parametric curve
+      ++ \spad{x = f(t)}, \spad{y = g(t)}; the default parameters \spad{1/2} 
+      ++ for the fraction and \spad{5/1} for the scale are used in the 
+      ++ \fakeAxiomFun{iClipParametric} subroutine, which is called by this
+      ++ function.
+    clipParametric: (PLOT,RN,RN) -> CLIPPED
+      ++ clipParametric(p,frac,sc) performs two-dimensional clipping on a 
+      ++ plot, p, from the domain \spadtype{Plot} for the parametric curve
+      ++ \spad{x = f(t)}, \spad{y = g(t)}; the fraction parameter is 
+      ++ specified by \spad{frac} and the scale parameter is specified 
+      ++ by \spad{sc} for use in the \fakeAxiomFun{iClipParametric} subroutine,
+      ++ which is called by this function.
+    clipWithRanges: (L L Pt,SF,SF,SF,SF) -> CLIPPED
+      ++ clipWithRanges(pointLists,xMin,xMax,yMin,yMax) performs clipping
+      ++ on a list of lists of points, \spad{pointLists}.  Clipping is 
+      ++ done within the specified ranges of \spad{xMin}, \spad{xMax} and 
+      ++ \spad{yMin}, \spad{yMax}.  This function is used internally by 
+      ++ the \fakeAxiomFun{iClipParametric} subroutine in this package.
+    clip: L Pt -> CLIPPED
+      ++ clip(l) performs two-dimensional clipping on a curve l, which is
+      ++ a list of points; the default parameters \spad{1/2} for the 
+      ++ fraction and \spad{5/1} for the scale are used in the 
+      ++ \fakeAxiomFun{iClipParametric} subroutine, which is called by this 
+      ++ function.
+    clip: L L Pt -> CLIPPED
+      ++ clip(ll) performs two-dimensional clipping on a list of lists 
+      ++ of points, \spad{ll}; the default parameters \spad{1/2} for
+      ++ the fraction and \spad{5/1} for the scale are used in the 
+      ++ \fakeAxiomFun{iClipParametric} subroutine, which is called by this 
+      ++ function.
+ 
+  Implementation ==> add
+    import PointPackage(DoubleFloat)
+    import ListFunctions2(Point DoubleFloat,DoubleFloat)
+ 
+    point:(SF,SF) -> Pt
+    intersectWithHorizLine:(SF,SF,SF,SF,SF) -> Pt
+    intersectWithVertLine:(SF,SF,SF,SF,SF) -> Pt
+    intersectWithBdry:(SF,SF,SF,SF,Pt,Pt) -> Pt
+    discardAndSplit: (L Pt,Pt -> B,SF,SF,SF,SF) -> L L Pt
+    norm: Pt -> SF
+    iClipParametric: (L L Pt,RN,RN) -> CLIPPED
+    findPt: L L Pt -> Union(Pt,"failed")
+    Fnan?: SF ->Boolean
+    Pnan?:Pt ->Boolean
+
+    Fnan? x == x~=x
+    Pnan? p == any?(Fnan?,p)
+   
+    iClipParametric(pointLists,fraction,scale) ==
+      -- error checks and special cases
+      (fraction < 0) or (fraction > 1) =>
+        error "clipDraw: fraction should be between 0 and 1"
+      empty? pointLists => [nil(),segment(0,0),segment(0,0)]
+      -- put all points together , sort them according to norm
+      sortedList := sort(norm(#1) < norm(#2),select(not Pnan? #1,concat pointLists))
+      empty? sortedList => [nil(),segment(0,0),segment(0,0)]
+      n := # sortedList 
+      num := numer fraction
+      den := denom fraction
+      clipNum := (n * num) quo den
+      lastN := n - 1 - clipNum
+      firstPt := first sortedList
+      xMin : SF := xCoord firstPt
+      xMax : SF := xCoord firstPt
+      yMin : SF := yCoord firstPt 
+      yMax : SF := yCoord firstPt
+      -- calculate min/max for the first (1-fraction)*N points
+      -- this contracts the range
+      -- this unnecessarily clips monotonic functions (step-function, x^(high power),etc.)
+      for k in 0..lastN  for pt in rest sortedList repeat
+        xMin := min(xMin,xCoord pt)
+        xMax := max(xMax,xCoord pt)
+        yMin := min(yMin,yCoord pt)
+        yMax := max(yMax,yCoord pt)
+      xDiff := xMax - xMin; yDiff := yMax - yMin
+      xDiff = 0 =>
+        yDiff = 0 =>
+          [pointLists,segment(xMin-1,xMax+1),segment(yMin-1,yMax+1)]
+        [pointLists,segment(xMin-1,xMax+1),segment(yMin,yMax)]
+      yDiff = 0 =>
+        [pointLists,segment(xMin,xMax),segment(yMin-1,yMax+1)]
+      numm := numer scale; denn := denom scale
+      -- now expand the range by scale
+      xMin := xMin - (numm :: SF) * xDiff / (denn :: SF)
+      xMax := xMax + (numm :: SF) * xDiff / (denn :: SF)
+      yMin := yMin - (numm :: SF) * yDiff / (denn :: SF)
+      yMax := yMax + (numm :: SF) * yDiff / (denn :: SF)
+      -- clip with the calculated range
+      newclip:=clipWithRanges(pointLists,xMin,xMax,yMin,yMax)
+      -- if we split the lists use the new clip
+      # (newclip.brans) > # pointLists   => newclip
+      -- calculate extents
+      xs :L SF:= map (xCoord,sortedList)
+      ys :L SF:= map (yCoord,sortedList)
+      xMin :SF :=reduce (min,xs)
+      yMin :SF :=reduce (min,ys)
+      xMax :SF :=reduce (max,xs)
+      yMax :SF :=reduce (max,ys) 
+      xseg:SEG SF :=xMin..xMax
+      yseg:SEG SF :=yMin..yMax
+      -- return original
+      [pointLists,xseg,yseg]@CLIPPED
+      
+
+
+ 
+    point(xx,yy) == point(l : L SF := [xx,yy])
+ 
+    intersectWithHorizLine(x1,y1,x2,y2,yy) ==
+      x1 = x2 => point(x1,yy)
+      point(x1 + (x2 - x1)*(yy - y1)/(y2 - y1),yy)
+ 
+    intersectWithVertLine(x1,y1,x2,y2,xx) ==
+      y1 = y2 => point(xx,y1)
+      point(xx,y1 + (y2 - y1)*(xx - x1)/(x2 - x1))
+ 
+    intersectWithBdry(xMin,xMax,yMin,yMax,pt1,pt2) ==
+      -- pt1 is in rectangle, pt2 is not
+      x1 := xCoord pt1; y1 := yCoord pt1
+      x2 := xCoord pt2; y2 := yCoord pt2
+      if y2 > yMax then
+        pt2 := intersectWithHorizLine(x1,y1,x2,y2,yMax)
+        x2 := xCoord pt2; y2 := yCoord pt2
+      if y2 < yMin then
+        pt2 := intersectWithHorizLine(x1,y1,x2,y2,yMin)
+        x2 := xCoord pt2; y2 := yCoord pt2
+      if x2 > xMax then
+        pt2 := intersectWithVertLine(x1,y1,x2,y2,xMax)
+        x2 := xCoord pt2; y2 := yCoord pt2
+      if x2 < xMin then
+        pt2 := intersectWithVertLine(x1,y1,x2,y2,xMin)
+      pt2
+ 
+    discardAndSplit(pointList,pred,xMin,xMax,yMin,yMax) ==
+      ans : L L Pt := nil()
+      list : L Pt  := nil()
+      lastPt? : B  := false
+      lastPt : Pt  := point(0,0)
+      while not empty? pointList repeat
+        pt := first pointList
+        pointList := rest pointList
+        pred(pt) =>
+          if (empty? list) and lastPt? then
+            bdryPt := intersectWithBdry(xMin,xMax,yMin,yMax,pt,lastPt)
+            -- print bracket [ coerce bdryPt ,coerce pt ]	
+            --list := cons(bdryPt,list)
+          list := cons(pt,list)
+        if not empty? list then
+          bdryPt := intersectWithBdry(xMin,xMax,yMin,yMax,first list,pt)
+	  -- print bracket [ coerce bdryPt,coerce first list]	
+          --list := cons(bdryPt,list)
+          ans := cons( list,ans)
+        lastPt := pt 
+	lastPt? := true
+	list := nil()
+      empty? list => ans
+      reverse_! cons(reverse_! list,ans)
+ 
+    clip(plot,fraction,scale) ==
+--      sayBrightly(["   clip: "::OutputForm]$List(OutputForm))$Lisp
+      (fraction < 0) or (fraction > 1/2) =>
+        error "clipDraw: fraction should be between 0 and 1/2"
+      xVals := xRange plot
+      empty?(pointLists := listBranches plot) =>
+        [nil(),xVals,segment(0,0)]
+      more?(pointLists := listBranches plot,1) =>
+        error "clipDraw: plot has more than one branch"
+      empty?(pointList := first pointLists) =>
+        [nil(),xVals,segment(0,0)]
+      sortedList := sort(yCoord(#1) < yCoord(#2),pointList)
+      n := # sortedList; num := numer fraction; den := denom fraction
+      clipNum := (n * num) quo den
+      -- throw out points with large and small y-coordinates
+      yMin := yCoord(sortedList.clipNum)
+      yMax := yCoord(sortedList.(n - 1 - clipNum))
+      if Fnan? yMin then yMin : SF := 0
+      if Fnan? yMax then yMax : SF := 0
+      (yDiff := yMax - yMin) = 0 =>
+        [pointLists,xRange plot,segment(yMin - 1,yMax + 1)]
+      numm := numer scale; denn := denom scale
+      xMin := lo xVals; xMax := hi xVals
+      yMin := yMin - (numm :: SF) * yDiff / (denn :: SF)
+      yMax := yMax + (numm :: SF) * yDiff / (denn :: SF)
+      lists := discardAndSplit(pointList,_
+         (yCoord(#1) < yMax) and (yCoord(#1) > yMin),xMin,xMax,yMin,yMax)
+      yMin := yCoord(sortedList.clipNum)
+      yMax := yCoord(sortedList.(n - 1 - clipNum))
+      if Fnan? yMin then yMin : SF := 0
+      if Fnan? yMax then yMax : SF := 0
+      for list in lists repeat
+        for pt in list repeat
+          if not Fnan?(yCoord pt) then
+            yMin := min(yMin,yCoord pt)
+            yMax := max(yMax,yCoord pt)
+      [lists,xVals,segment(yMin,yMax)]
+ 
+    clip(plot:PLOT) == clip(plot,1/4,5/1)
+ 
+    norm(pt) == 
+      x := xCoord(pt); y := yCoord(pt)
+      if Fnan? x then
+        if Fnan? y then
+          r:SF := 0
+        else
+          r:SF := y**2
+      else
+        if Fnan? y then
+          r:SF := x**2
+        else
+          r:SF := x**2 + y**2
+      r
+ 
+    findPt lists ==
+      for list in lists repeat
+        not empty? list => 
+	     for p in list repeat 
+               not Pnan? p => return p
+      "failed"
+
+    clipWithRanges(pointLists,xMin,xMax,yMin,yMax) ==
+      lists : L L Pt := nil()
+      for pointList in pointLists repeat
+        lists := concat(lists,discardAndSplit(pointList,_
+           (xCoord(#1) <= xMax) and (xCoord(#1) >= xMin) and _
+           (yCoord(#1) <= yMax) and (yCoord(#1) >= yMin), _
+           xMin,xMax,yMin,yMax))
+      (pt := findPt lists) case "failed" =>
+        [nil(),segment(0,0),segment(0,0)]
+      firstPt := pt :: Pt
+      xMin : SF := xCoord firstPt; xMax : SF := xCoord firstPt
+      yMin : SF := yCoord firstPt; yMax : SF := yCoord firstPt
+      for list in lists repeat
+        for pt in list repeat
+          if not Pnan? pt then
+            xMin := min(xMin,xCoord pt)
+            xMax := max(xMax,xCoord pt)
+            yMin := min(yMin,yCoord pt)
+            yMax := max(yMax,yCoord pt)
+      [lists,segment(xMin,xMax),segment(yMin,yMax)]
+ 
+    clipParametric(plot,fraction,scale) ==
+      iClipParametric(listBranches plot,fraction,scale)
+ 
+    clipParametric plot == clipParametric(plot,1/2,5/1)
+ 
+    clip(l: L Pt)   == iClipParametric(list l,1/2,5/1)
+    clip(l: L L Pt) == iClipParametric(l,1/2,5/1)
+
+@
+\section{License}
+<<license>>=
+--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
+--All rights reserved.
+--
+--Redistribution and use in source and binary forms, with or without
+--modification, are permitted provided that the following conditions are
+--met:
+--
+--    - Redistributions of source code must retain the above copyright
+--      notice, this list of conditions and the following disclaimer.
+--
+--    - Redistributions in binary form must reproduce the above copyright
+--      notice, this list of conditions and the following disclaimer in
+--      the documentation and/or other materials provided with the
+--      distribution.
+--
+--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
+--      names of its contributors may be used to endorse or promote products
+--      derived from this software without specific prior written permission.
+--
+--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+@
+<<*>>=
+<<license>>
+
+<<package CLIP TwoDimensionalPlotClipping>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}