axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / EXPRTUBE.spad

)abbrev package EXPRTUBE ExpressionTubePlot
++ Author: Clifton J. Williamson
++ Date Created: Bastille Day 1989
++ Date Last Updated: 5 June 1990
++ Description:
++ Package for constructing tubes around 3-dimensional parametric curves.

ExpressionTubePlot() : SIG == CODE 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
 
  SIG ==> 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.
 
  CODE ==> 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 that 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 == s +-> project(x,s)
 
    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")