path: root/src/algebra/tube.spad.pamphlet

\documentclass{article}
\usepackage{open-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}