\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} <>= )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} <>= )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} <>= )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} <>= )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} <>= --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. @ <<*>>= <> <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}