Initial population.

1 files changed, 508 insertions, 0 deletions

diff --git a/src/algebra/tube.spad.pamphlet b/src/algebra/tube.spad.pamphlet
new file mode 100644
index 00000000..9d74406f
--- /dev/null
+++ b/src/algebra/tube.spad.pamphlet
@@ -0,0 +1,508 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra tube.spad}
+\author{Clifton J. Williamson}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{domain TUBE TubePlot}
+<<domain TUBE TubePlot>>=
+)abbrev domain TUBE TubePlot
+++ Author: Clifton J. Williamson
+++ Date Created: Bastille Day 1989
+++ Date Last Updated: 5 June 1990
+++ Keywords:
+++ Examples:
+++ Description: 
+++   Package for constructing tubes around 3-dimensional parametric curves.
+++ Domain of tubes around 3-dimensional parametric curves.
+TubePlot(Curve): Exports == Implementation where
+  Curve : PlottableSpaceCurveCategory
+  B   ==> Boolean
+  L   ==> List
+  Pt  ==> Point DoubleFloat
+ 
+  Exports ==> with
+    getCurve: % -> Curve
+      ++ getCurve(t) returns the \spadtype{PlottableSpaceCurveCategory}
+      ++ representing the parametric curve of the given tube plot t.
+    listLoops: % -> L L Pt
+      ++ listLoops(t) returns the list of lists of points, or the 'loops',
+      ++ of the given tube plot t.
+    closed?: % -> B
+      ++ closed?(t) tests whether the given tube plot t is closed. 
+    open?: % -> B
+      ++ open?(t) tests whether the given tube plot t is open. 
+    setClosed: (%,B) -> B
+      ++ setClosed(t,b) declares the given tube plot t to be closed if
+      ++ b is true, or if b is false, t is set to be open.
+    tube: (Curve,L L Pt,B) -> %
+      ++ tube(c,ll,b) creates a tube of the domain \spadtype{TubePlot} from a
+      ++ space curve c of the category \spadtype{PlottableSpaceCurveCategory},
+      ++ a list of lists of points (loops) ll and a boolean b which if
+      ++ true indicates a closed tube, or if false an open tube.
+ 
+  Implementation ==> add
+ 
+--% representation
+ 
+    Rep := Record(parCurve:Curve,loops:L L Pt,closedTube?:B)
+ 
+    getCurve plot == plot.parCurve
+ 
+    listLoops plot == plot.loops
+ 
+    closed? plot == plot.closedTube?
+    open? plot   == not plot.closedTube?
+ 
+    setClosed(plot,flag) == plot.closedTube? := flag
+ 
+    tube(curve,ll,b) == [curve,ll,b]
+
+@
+\section{package TUBETOOL TubePlotTools}
+<<package TUBETOOL TubePlotTools>>=
+)abbrev package TUBETOOL TubePlotTools
+++ Author: Clifton J. Williamson
+++ Date Created: Bastille Day 1989
+++ Date Last Updated: 5 June 1990
+++ Keywords:
+++ Examples:
+++ Description: 
+++   Tools for constructing tubes around 3-dimensional parametric curves.
+TubePlotTools(): Exports == Implementation where
+  I   ==> Integer
+  SF  ==> DoubleFloat
+  L   ==> List
+  Pt  ==> Point SF
+ 
+  Exports ==> with
+    point: (SF,SF,SF,SF) -> Pt
+      ++ point(x1,x2,x3,c) creates and returns a point from the three 
+      ++ specified coordinates \spad{x1}, \spad{x2}, \spad{x3}, and also
+      ++ a fourth coordinate, c, which is generally used to specify the 
+      ++ color of the point.
+    "*" : (SF,Pt) -> Pt
+      ++ s * p returns a point whose coordinates are the scalar multiple
+      ++ of the point p by the scalar s, preserving the color, or fourth
+      ++ coordinate, of p.
+    "+" : (Pt,Pt) -> Pt
+      ++ p + q computes and returns a point whose coordinates are the sums
+      ++ of the coordinates of the two points \spad{p} and \spad{q}, using
+      ++ the color, or fourth coordinate, of the first point \spad{p}
+      ++ as the color also of the point \spad{q}.
+    "-" : (Pt,Pt) -> Pt
+      ++ p - q computes and returns a point whose coordinates are the
+      ++ differences of the coordinates of two points \spad{p} and \spad{q},
+      ++ using the color, or fourth coordinate, of the first point \spad{p}
+      ++ as the color also of the point \spad{q}.
+    dot : (Pt,Pt) -> SF
+      ++ dot(p,q) computes the dot product of the two points \spad{p}
+      ++ and \spad{q} using only the first three coordinates, and returns
+      ++ the resulting \spadtype{DoubleFloat}.
+    cross : (Pt,Pt) -> Pt
+      ++ cross(p,q) computes the cross product of the two points \spad{p}
+      ++ and \spad{q} using only the first three coordinates, and keeping
+      ++ the color of the first point p.  The result is returned as a point.
+    unitVector: Pt -> Pt
+      ++ unitVector(p) creates the unit vector of the point p and returns
+      ++ the result as a point.  Note: \spad{unitVector(p) = p/|p|}.
+    cosSinInfo: I -> L L SF
+      ++ cosSinInfo(n) returns the list of lists of values for n, in the
+      ++ form:  \spad{[[cos(n - 1) a,sin(n - 1) a],...,[cos 2 a,sin 2 a],[cos a,sin a]]}
+      ++ where \spad{a = 2 pi/n}.  Note: n should be greater than 2.
+    loopPoints: (Pt,Pt,Pt,SF,L L SF) -> L Pt
+      ++ loopPoints(p,n,b,r,lls) creates and returns a list of points 
+      ++ which form the loop with radius r, around the center point
+      ++ indicated by the point p, with the principal normal vector of
+      ++ the space curve at point p given by the point(vector) n, and the 
+      ++ binormal vector given by the point(vector) b, and a list of lists,
+      ++ lls, which is the \spadfun{cosSinInfo} of the number of points
+      ++ defining the loop.
+ 
+  Implementation ==> add
+    import PointPackage(SF)
+ 
+    point(x,y,z,c) == point(l : L SF := [x,y,z,c])
+ 
+    getColor: Pt -> SF
+    getColor pt == (maxIndex pt > 3 => color pt; 0)
+ 
+    getColor2: (Pt,Pt) -> SF
+    getColor2(p0,p1) ==
+      maxIndex p0 > 3 => color p0
+      maxIndex p1 > 3 => color p1
+      0
+ 
+    a * p ==
+      l : L SF := [a * xCoord p,a * yCoord p,a * zCoord p,getColor p]
+      point l
+ 
+    p0 + p1 ==
+      l : L SF := [xCoord p0 + xCoord p1,yCoord p0 + yCoord p1,_
+                   zCoord p0 + zCoord p1,getColor2(p0,p1)]
+      point l
+ 
+    p0 - p1 ==
+      l : L SF := [xCoord p0 - xCoord p1,yCoord p0 - yCoord p1,_
+                   zCoord p0 - zCoord p1,getColor2(p0,p1)]
+      point l
+ 
+    dot(p0,p1) ==
+      (xCoord p0 * xCoord p1) + (yCoord p0 * yCoord p1) +_
+        (zCoord p0 * zCoord p1)
+ 
+    cross(p0,p1) ==
+      x0 := xCoord p0; y0 := yCoord p0; z0 := zCoord p0;
+      x1 := xCoord p1; y1 := yCoord p1; z1 := zCoord p1;
+      l : L SF := [y0 * z1 - y1 * z0,z0 * x1 - z1 * x0,_
+                   x0 * y1 - x1 * y0,getColor2(p0,p1)]
+      point l
+ 
+    unitVector p == (inv sqrt dot(p,p)) * p
+ 
+    cosSinInfo n ==
+      ans : L L SF := nil()
+      theta : SF := 2 * pi()/n
+      for i in 1..(n-1) repeat             --!! make more efficient
+        angle := i * theta
+        ans := concat([cos angle,sin angle],ans)
+      ans
+ 
+    loopPoints(ctr,pNorm,bNorm,rad,cosSin) ==
+      ans : L Pt := nil()
+      while not null cosSin repeat
+        cossin := first cosSin; cos := first cossin; sin := second cossin
+        ans := cons(ctr + rad * (cos * pNorm + sin * bNorm),ans)
+        cosSin := rest cosSin
+      pt := ctr + rad * pNorm
+      concat(pt,concat(ans,pt))
+
+@
+\section{package EXPRTUBE ExpressionTubePlot}
+<<package EXPRTUBE ExpressionTubePlot>>=
+)abbrev package EXPRTUBE ExpressionTubePlot
+++ Author: Clifton J. Williamson
+++ Date Created: Bastille Day 1989
+++ Date Last Updated: 5 June 1990
+++ Keywords:
+++ Examples:
+++ Package for constructing tubes around 3-dimensional parametric curves.
+ExpressionTubePlot(): Exports == Implementation where
+  B   ==> Boolean
+  I   ==> Integer
+  FE  ==> Expression Integer
+  SY  ==> Symbol
+  SF  ==> DoubleFloat
+  L   ==> List
+  S   ==> String
+  SEG ==> Segment
+  F2F ==> MakeFloatCompiledFunction(FE)
+  Pt  ==> Point SF
+  PLOT3 ==> Plot3D
+  TUBE  ==> TubePlot Plot3D
+ 
+  Exports ==> with
+    constantToUnaryFunction: SF -> (SF -> SF)
+      ++ constantToUnaryFunction(s) is a local function which takes the
+      ++ value of s, which may be a function of a constant, and returns
+      ++ a function which always returns the value \spadtype{DoubleFloat} s.
+    tubePlot: (FE,FE,FE,SF -> SF,SEG SF,SF -> SF,I) -> TUBE
+      ++ tubePlot(f,g,h,colorFcn,a..b,r,n) puts a tube of radius r(t) with
+      ++ n points on each circle about the curve \spad{x = f(t)}, 
+      ++ \spad{y = g(t)}, \spad{z = h(t)} for t in \spad{[a,b]}. 
+      ++ The tube is considered to be open.
+    tubePlot: (FE,FE,FE,SF -> SF,SEG SF,SF -> SF,I,S) -> TUBE
+      ++ tubePlot(f,g,h,colorFcn,a..b,r,n,s) puts a tube of radius \spad{r(t)}
+      ++ with n points on each circle about the curve \spad{x = f(t)}, 
+      ++ \spad{y = g(t)},
+      ++ \spad{z = h(t)} for t in \spad{[a,b]}. If s = "closed", the tube is
+      ++ considered to be closed; if s = "open", the tube is considered
+      ++ to be open.
+    tubePlot: (FE,FE,FE,SF -> SF,SEG SF,SF,I) -> TUBE
+      ++ tubePlot(f,g,h,colorFcn,a..b,r,n) puts a tube of radius r with
+      ++ n points on each circle about the curve \spad{x = f(t)}, 
+      ++ \spad{y = g(t)}, \spad{z = h(t)} for t in \spad{[a,b]}. 
+      ++ The tube is considered to be open.
+    tubePlot: (FE,FE,FE,SF -> SF,SEG SF,SF,I,S) -> TUBE
+      ++ tubePlot(f,g,h,colorFcn,a..b,r,n,s) puts a tube of radius r with
+      ++ n points on each circle about the curve \spad{x = f(t)}, 
+      ++ \spad{y = g(t)}, \spad{z = h(t)} for t in \spad{[a,b]}. 
+      ++ If s = "closed", the tube is
+      ++ considered to be closed; if s = "open", the tube is considered
+      ++ to be open.
+ 
+  Implementation ==> add
+    import Plot3D
+    import F2F
+    import TubePlotTools
+ 
+--% variables
+ 
+    getVariable: (FE,FE,FE) -> SY
+    getVariable(x,y,z) ==
+      varList1 := variables x
+      varList2 := variables y
+      varList3 := variables z
+      (not (# varList1 <= 1)) or (not (# varList2 <= 1)) or _
+       (not (# varList3 <= 1)) =>
+        error "tubePlot: only one variable may be used"
+      null varList1 =>
+        null varList2 =>
+          null varList3 =>
+            error "tubePlot: a variable must appear in functions"
+          first varList3
+        t2 := first varList2
+        null varList3 => t2
+        not (first varList3 = t2) =>
+          error "tubePlot: only one variable may be used"
+      t1 := first varList1
+      null varList2 =>
+        null varList3 => t1
+        not (first varList3 = t1) =>
+          error "tubePlot: only one variable may be used"
+        t1
+      t2 := first varList2
+      null varList3 =>
+        not (t1 = t2) =>
+          error "tubePlot: only one variable may be used"
+        t1
+      not (first varList3 = t1) or not (t2 = t1) =>
+        error "tubePlot: only one variable may be used"
+      t1
+ 
+--% tubes: variable radius
+ 
+    tubePlot(x:FE,y:FE,z:FE,colorFcn:SF -> SF,_
+             tRange:SEG SF,radFcn:SF -> SF,n:I,string:S) ==
+      -- check value of n
+      n < 3 => error "tubePlot: n should be at least 3"
+      -- check string
+      flag : B :=
+        string = "closed" => true
+        string = "open" => false
+        error "tubePlot: last argument should be open or closed"
+      -- check variables
+      t := getVariable(x,y,z)
+      -- coordinate functions
+      xFunc := makeFloatFunction(x,t)
+      yFunc := makeFloatFunction(y,t)
+      zFunc := makeFloatFunction(z,t)
+      -- derivatives of coordinate functions
+      xp := differentiate(x,t)
+      yp := differentiate(y,t)
+      zp := differentiate(z,t)
+      -- derivative of arc length
+      sp := sqrt(xp ** 2 + yp ** 2 + zp ** 2)
+      -- coordinates of unit tangent vector
+      Tx := xp/sp; Ty := yp/sp; Tz := zp/sp
+      -- derivatives of coordinates of unit tangent vector
+      Txp := differentiate(Tx,t)
+      Typ := differentiate(Ty,t)
+      Tzp := differentiate(Tz,t)
+      -- K = curvature = length of curvature vector
+      K := sqrt(Txp ** 2 + Typ ** 2 + Tzp ** 2)
+      -- coordinates of principal normal vector
+      Nx := Txp / K; Ny := Typ / K; Nz := Tzp / K
+      -- functions SF->SF giving coordinates of principal normal vector
+      NxFunc := makeFloatFunction(Nx,t);
+      NyFunc := makeFloatFunction(Ny,t);
+      NzFunc := makeFloatFunction(Nz,t);
+      -- coordinates of binormal vector
+      Bx := Ty * Nz - Tz * Ny
+      By := Tz * Nx - Tx * Nz
+      Bz := Tx * Ny - Ty * Nx
+      -- functions SF -> SF giving coordinates of binormal vector
+      BxFunc := makeFloatFunction(Bx,t);
+      ByFunc := makeFloatFunction(By,t);
+      BzFunc := makeFloatFunction(Bz,t);
+      -- create Plot3D
+      parPlot := plot(xFunc,yFunc,zFunc,colorFcn,tRange)
+      tvals := first tValues parPlot
+      curvePts := first listBranches parPlot
+      cosSin := cosSinInfo n
+      loopList : L L Pt := nil()
+      while not null tvals repeat
+        -- note: tvals and curvePts have the same number of elements
+        tval := first tvals; tvals := rest tvals
+        ctr := first curvePts; curvePts := rest curvePts
+        pNormList : L SF :=
+          [NxFunc tval,NyFunc tval,NzFunc tval,colorFcn tval]
+        pNorm : Pt := point pNormList
+        bNormList : L SF :=
+          [BxFunc tval,ByFunc tval,BzFunc tval,colorFcn tval]
+        bNorm : Pt := point bNormList
+        lps := loopPoints(ctr,pNorm,bNorm,radFcn tval,cosSin)
+        loopList := cons(lps,loopList)
+      tube(parPlot,reverse_! loopList,flag)
+ 
+    tubePlot(x:FE,y:FE,z:FE,colorFcn:SF -> SF,_
+             tRange:SEG SF,radFcn:SF -> SF,n:I) ==
+      tubePlot(x,y,z,colorFcn,tRange,radFcn,n,"open")
+ 
+--% tubes: constant radius
+ 
+    project: (SF,SF) -> SF
+    project(x,y) == x
+ 
+    constantToUnaryFunction x == project(x,#1)
+ 
+    tubePlot(x:FE,y:FE,z:FE,colorFcn:SF -> SF,_
+             tRange:SEG SF,rad:SF,n:I,s:S) ==
+      tubePlot(x,y,z,colorFcn,tRange,constantToUnaryFunction rad,n,s)
+ 
+    tubePlot(x:FE,y:FE,z:FE,colorFcn:SF -> SF,_
+             tRange:SEG SF,rad:SF,n:I) ==
+      tubePlot(x,y,z,colorFcn,tRange,rad,n,"open")
+
+@
+\section{package NUMTUBE NumericTubePlot}
+<<package NUMTUBE NumericTubePlot>>=
+)abbrev package NUMTUBE NumericTubePlot
+++ Author: Clifton J. Williamson
+++ Date Created: Bastille Day 1989
+++ Date Last Updated: 5 June 1990
+++ Keywords:
+++ Examples:
+++ Package for constructing tubes around 3-dimensional parametric curves.
+NumericTubePlot(Curve): Exports == Implementation where
+  Curve : PlottableSpaceCurveCategory
+  B   ==> Boolean
+  I   ==> Integer
+  SF  ==> DoubleFloat
+  L   ==> List
+  S   ==> String
+  SEG ==> Segment
+  Pt  ==> Point SF
+  TUBE ==> TubePlot Curve
+  Triad ==> Record(tang:Pt,norm:Pt,bin:Pt)
+ 
+  Exports ==> with
+    tube: (Curve,SF,I) -> TUBE
+      ++ tube(c,r,n) creates a tube of radius r around the curve c.
+ 
+  Implementation ==> add
+    import TubePlotTools
+ 
+    LINMAX  := convert(0.995)@SF
+    XHAT := point(1,0,0,0)
+    YHAT := point(0,1,0,0)
+    PREV0 := point(1,1,0,0)
+    PREV := PREV0
+ 
+    colinearity: (Pt,Pt) -> SF
+    colinearity(x,y) == dot(x,y)**2/(dot(x,x) * dot(y,y))
+ 
+    orthog: (Pt,Pt) -> Pt
+    orthog(x,y) ==
+      if colinearity(x,y) > LINMAX then y := PREV
+      if colinearity(x,y) > LINMAX then
+        y := (colinearity(x,XHAT) < LINMAX => XHAT; YHAT)
+      a := -dot(x,y)/dot(x,x)
+      PREV := a*x + y
+ 
+    poTriad:(Pt,Pt,Pt) -> Triad
+    poTriad(pl,po,pr) ==
+      -- use divided difference for t.
+      t := unitVector(pr - pl)
+      -- compute n as orthogonal to t in plane containing po.
+      pol := pl - po
+      n   := unitVector orthog(t,pol)
+      [t,n,cross(t,n)]
+ 
+    curveTriads: L Pt -> L Triad
+    curveTriads l ==
+      (k := #l) < 2 => error "Need at least 2 points to specify a curve"
+      PREV := PREV0
+      k = 2 =>
+        t := unitVector(second l - first l)
+        n := unitVector(t - XHAT)
+        b := cross(t,n)
+        triad : Triad := [t,n,b]
+        [triad,triad]
+      -- compute interior triads using divided differences
+      midtriads : L Triad :=
+        [poTriad(pl,po,pr) for pl in l for po in rest l _
+               for pr in rest rest l]
+      -- compute first triad using a forward difference
+      x := first midtriads
+      t := unitVector(second l - first l)
+      n := unitVector orthog(t,x.norm)
+      begtriad : Triad := [t,n,cross(t,n)]
+      -- compute last triad using a backward difference
+      x := last midtriads
+      -- efficiency!!
+      t := unitVector(l.k - l.(k-1))
+      n := unitVector orthog(t,x.norm)
+      endtriad : Triad := [t,n,cross(t,n)]
+      concat(begtriad,concat(midtriads,endtriad))
+ 
+    curveLoops: (L Pt,SF,I) -> L L Pt
+    curveLoops(pts,r,nn) ==
+      triads := curveTriads pts
+      cosSin := cosSinInfo nn
+      loops : L L Pt := nil()
+      for pt in pts for triad in triads repeat
+        n := triad.norm; b := triad.bin
+        loops := concat(loopPoints(pt,n,b,r,cosSin),loops)
+      reverse_! loops
+ 
+    tube(curve,r,n) ==
+      n < 3 => error "tube: n should be at least 3"
+      brans := listBranches curve
+      loops : L L Pt := nil()
+      for bran in brans repeat
+        loops := concat(loops,curveLoops(bran,r,n))
+      tube(curve,loops,false)
+
+@
+\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>>
+ 
+<<domain TUBE TubePlot>>
+<<package TUBETOOL TubePlotTools>>
+<<package EXPRTUBE ExpressionTubePlot>>
+<<package NUMTUBE NumericTubePlot>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}