Tubagraphics 
(June 1999) by M A Berger, Mathematics, University College London

Print[
      "Tuba package written by Prof M A Berger, Mathematics, University \
College London. Last revised June 1999.
		Please send comments and suggestions to m.berger@ucl.ac.uk"]

BeginPackage["Topographics`",  "Graphics`Colors`"];  
    Off[General::spell1];

Usage Messages

normalize, rotate, curveinterp, todata, slidebraid

\(normalize::usage\  = \ 
\ "normalize[v_] := N[v/ Sqrt[v.v]]
		"; 
rotate::usage\  = \ 
\ 
      "rotate[v, theta] rotates the 2d vector v anticlockwise by an angle theta
		"; 
curveinterp::usage\  = \ 
\ 
      "curveinterp[curve, ncurve] changes the number of points evaluated along a curve. It takes a curve through 3d space stored as an(N,3) list where N need not be known. It then uses (cubic) Interpolation to return the same curve as an (ncurve,3) list.
		"; 
todata::usage\  = 

      "todata[fun, np] takes a function and returns a list of the values
N[fun[0]], N[fun[1/(np-1)]], ..., fun[1]."; 

    slidebraid::usage = 
      "slidebraid[b, {20,0,0}] moves braid b by {20,0,0} "; \)

tube, quicktube,  checkedcurve

\(\(tube::usage\  = 
\ 
      \*""tube[curve, color, shape] takes a curve represented by either a function or a list (a list of {x,y,z} coordinates of arbitrary length), and draws a tube around it. It returns a Graphics3D object; use Show[[tube[curve, color,shape]]
		   to display (or try quicktube[curve]). Tube will interpolate nalong points along the curve, and then find naround points in the meridional direction (see defaults below).
		If curve is a function, it has domain and range curve:[0,1] -> (ℝ^3). 
		
		Let the variable s measure position the short way around, going from 0 to 1. Also let t measure position the long way around (along the axis), also varying from 0 to 1.
		Both shape and color can be constants or functions of [s,t].
		If shape is a number, then the radius of the tube will be shape. Otherwise, the radius is given by shape[s,t].
		
		Color can be a single graphics directive or list of graphics directives. 
		An assignment like color[s_,t_] := {EdgeForm[],Hue[s]} 
		will automatically be translated to {EdgeForm[],SurfaceColor[Hue[s]]}. 
		This works for RGBColor, Hue, and CMYKColor, as well as the
	   color names listed in AllColors (part of Graphics`Colors package automatically loaded by tuba.m).
		
		For a frenettube (the default), numerical differentiation is used to find normal and binormals, with avoidance of flips at inflection points. Other algorithms are available. If the first and last points in curve are identical then tube assumes a closed curve. 
	
			There are several options:
		1. naroundtube -> na gives na points in azimuthal direction (default = 12);
		2. nalongtube  -> nc gives nc points in axial direction along curve (default 40);
		3. tuberange   -> myfun will only plot parts of the surface inside the region myfun(x,y,z) ≥ 0. 
		4. tubehelicity-> hh will add hh units of twist (in units of 2π radians). The default is hh = 0; this gives zero framing (zero helicity).
 	  5. tubealgorithm -> xytube will plot circles perpendicular to curve in xy plane. tubealgorithm -> {a,b,c} gets normal by crossing tangent with fixed vector {a,b,c}. Default
         is frenettube, which uses normal and binormal for circles.""; \)\)

tubedefaults::usage = "default tube options";
tuberange::usage= "tube option";
tubehelicity::usage ="tube option, default = 0";
tubealgorithm::usage = "tube option, default = frenettube = 1.(xytube = 2)";
naroundtube::usage = "tube option, default = 12";
nalongtube::usage = "tube option, default = 40 ";
frenettube::usage = "tubalgorithm option (=1)";
xytube::usage = "tubealgorithm option (=2)";

\(quicktube::usage\  = \ 
      "quicktube[curve]; shows radius 1 tube around curve with no boxes and white surface color."; 

    checkedcurve::usage\  = \ 
      "checkedcurve[curve, n] does a moderate amount of checking to see that curve is either a 3-vector valued function of t (at least curve[t=0] is a 3-vector), or that curve is a list of 3-vectors. It then ouputs an interpolated list of 3-vector positions of length n"; \)

fluxlines

fluxlines::usage  =  
      "Suppose you want to draw one set of lines from a tube surrounding a \
curve. 
		fluxlines[curve, radius, nalong, naround] creates a list of two Graphics3D \
objects: each consisting of one set of lines (parallel or perp to axis). 
		There are two options:
		1 fluxthickness  -> th gives thickness of lines th (default 0.005);
		2 fluxcolor  ->c determines color c of lines (default GrayLevel[0]).";

\(fluxthickness::usage\  = \ \ "Option for fluxlines"; 

    fluxcolor::usage\  = \ \ "Option for fluxlines"; \)

arrow, sidearrow

\(arrowcolor::usage\  = \ \ \ "arrow and sidearrow option"; 

    arrowthickness::usage\  = \ \ \ "arrow and sidearrow option"; 

    arrow::usage\  = \ 
\ 
      "arrow[p,q,n] draws an arrow with tail at p and tip at q and normal to its surface along n (actually n perp to q-p). 
				A graphics3D object is output. To show, use Show[arrow[p,q,n]].
		There are two options:
		
		1. arrowcolor -> color gives surface color;
		
		2. arrowthickness  -> th The thickness of the arrow is th|q-p| (default th = 0.1);
		"; 
sidearrow::usage\  = \ 
\ 
      "sidearrow[fun, t, eps, radius, theta] draws a red arrow with tail at fun[t] and tip at
		fun[t+eps]. The arrow is displaced by radius in the default viewing direction (ViewPoint). Use, e.g.
		SetOptions[Graphics3D, ViewPoint -> {-1,1,3}] to change ViewPoint so that sidearrow knows about it (changing from inside graphics3D will leave the default ViewPoint unchanged).
		Sidearrow is meant to stick onto a tube plot. You can rotate it around the tube by the angle theta."; \)

tress

braidcolors::usage =   "tress option";
pedestal::usage =   "tress option: pedestal -> {4,0.75} (the default), makes the margin between the side of the pedestal and the tubes equal to 4 * radius (of a tube). Also the thickness of the pedestal becomes 0.75 radius";
nicecolors::usage = "list of common colors";
tresswidth::usage = "returns width of bounding box for tress";
tress::usage =
 "tress[braid] takes a set of curves and draws tubes for each curve with a pedestal at bottom and top. Options:
		braidcolours->default selects colours for you (with repeats for more than 10 colours). 
		braidcolours->list colours tube1 with first colour in list, e.g. RGBColor[1,1,0], etc. 
		pedestal -> {4,0.75} (the default), makes the margin between the side of the pedestal and the tubes 
		equal to 4 * radius (of a tube). Also the thickness of the pedestal becomes 0.75 radius";

Code

Begin["`Private`"];

Basic Data Type: curvedata

curvedata  =  x : {{_?NumericQ, _?NumericQ, _?NumericQ} .. };

User functions: normalize, rotate, curveinterp, todata

normalize[v_]  :=  N[v/ Sqrt[v . v]]\(rotate[v_, \ theta_]\  := 
      \ {{Cos[theta], \ Sin[theta]}, {\(-Sin[theta]\), \ Cos[theta]}}\  . \ 
        v; 
curveinterp[curve_, \ ncurve_]\  := \ 
      Module[
		{tran, \ inter, \ t, \ morepoints, \ npoints, \ 
          interporder, \ i}, \ tran\  = \ Transpose[curve]; \ 
		
        npoints\  = \ Length[curve]; 
        interporder\  = \ If[npoints > 3, 3, npoints - 1]; 
		\ 
        inter\  = \ 
          Map[Interpolation[#, \ InterpolationOrder -> interporder]&, \ 
            tran]; \ 
		
        t[i_]\  := \ 
          N[\((i - 1)\) \((npoints\  - \ 1)\)/\((ncurve\  - \ 1)\)\  + \ 1]; 
        
		morepoints[ifun_]\  := \ 
		
          Append[Prepend[\ \ 
              Table[N[ifun[t[i]]], \ {i, \ 2, \ ncurve\  - 1}], \ 
					
              N[ifun[1]]], N[ifun[npoints]]\ ]; 
\ \ \ \ \ \ 
        Transpose[\ Map[morepoints, \ inter]\ ]
		]\)\(todata[fun_Symbol, \ nx_]\  := \ 
      Table[N[fun[\((i - 1)\)/\((nx - 1)\)]], \ {i, \ 1, \ nx}]; 

    todata[fun_Symbol, \ nx_, \ ny_] := \ 
		\ 
      Table[fun[\((i - 1)\)/\((nx - 1)\)\ , \((j - 1)\)/\((ny - 1)\)], 
        \ {i, nx}, \ {j, \ ny}]; \)

Derivative routines

weirdvector  :=   normalize[{Pi +  1, 3.7, 1.9}];
	\(deriv[c_]\  := 
          0.5\ *\ If[First[c]\  == \ Last[c], \ 
				\ 
              Module[{ccut\  = \ Drop[c, \ \(-1\)], \ derivcut}, \ \ 
                derivcut\  = \ RotateLeft[ccut]\  - \ RotateRight[ccut]; 

                					\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Append[derivcut, \ First[derivcut]]\ ], 
			\ \ \ 
              Module[{np\  = \ Length[c], \ i}, \ 

                						\ \ \ \ \ \ \ \ \ \ \ 
                N[Append[
                    Prepend[
                      Table[\((
                          c[[i + 1]]\  - 
                            \ c[[i - 1]])
                          \), \ {i, \ 2, \ np\  - 1}], 

                      									\ \ \ \ 
                      2 \((\ c[[3]] - 
                              c[[1]])
                            \)\  - \ 
                        \((c[[4]] - 
                            c[[2]])\)], 
                    \ 2 \((c[[np]]\  - \ 
                            c[[np - 2]])
                          \)\  - \ 
                      \((c[[np - 1
                              ]]\ \  - \ \ 
                          c[[np - 3]])
                        \)\ ]]]\ ]\);
    
  
                
     \(normalvector[tangent_]\  := \ 
          Module[{np\  = \ Length[tangent], \ n1, \ vv, \ vw, \ i}, 

            				\ n1\  = \ deriv[tangent]; 
		
            n1[[1]]\  = \ 
              Module[{v2\  = \ 
                    N[\ n1[[1]] . \ 
                        n1[[1]]]}, \ 
                If[v2\  != \ 0, 
                  normalize[\ n1[[1]]], 
                  
		normalize[
                    Cross[tangent[[1]], 
                      \ weirdvector]]]\ ]\ ; \ 
		
            Do[
			{
                vv\  = \ 
                  n1[[i]] . 
                    n1[[i]]; \ 
                If[vv\  > \ 10^\(-6\), 
                  n1[[i]]\  = \ 
                    n1[[i]]/Sqrt[vv], \ 
                  n1[[i]]\  = \ 
                    n1[[i - 1]]]; \ 

                						
                vw\  = \ 
                  n1[[i]] . 
                    n1[[i - 1]]; \ 
                If[vw\  < \ 0, \ 
                  n1[[i]]\  = \ 
                    \(-n1[[i]]\)]}, 
              \ {i, \ 2, \ np}
			]; \ n1]\);

Utility functions for restricting tubes to region specified by tuberange

\(interpoint[x_, \ y_, \ i_]\  := \ 
        Module[{j\  = \ Mod[i, \ 4]\  + \ 1}, 
		
          x[[j]]\  - \ 
            \((y[[j]]/
                  \((y[[j]]\  - \ 
                      y[[i]])\)\ )\) 
              \((x[[j]]\  - \ 
                  x[[i]])\)]\)
		\(slice[poly_, \ funvalues_]\  := \ 
        Module[{test, \ signs, \ edges, \ newpoly = {}, \ i}, \ 
          signs\  = \ Sign[funvalues]; \ 
          edges\  = 
			signs\ *\ RotateLeft[signs]; 
          Do[{\ If[signs[[i]]\  == \ 1, 
                \ newpoly\  = \ Append[newpoly, Extract[poly, i]]], 

              				\ \ \ \ \ 
              If[edges[[i]]\  == \ 
                  \(-1\), \ 
                newpoly\  = \ 
                  Append[newpoly, \ interpoint[poly, \ funvalues, \ i]\ ]
                    \ ]}, \ {i, \ 4}]; newpoly
		]\);

\(zap[poly_, \ funtest_]\  := \ 
        Module[{funvalues, howmanybadpoints}, \ 
          funvalues\  = \ Map[funtest, \ poly]; \ 
          howmanybadpoints\  = \ Count[funvalues, x_ /; x < 0]; \ 
          Which[howmanybadpoints\  == \ 0, \ poly, \ 

            	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ howmanybadpoints\  == \ \ 4, \ 
            Null, 
	\ \ \ \ \ \ \ \ \ \ \ True, \ slice[poly, \ funvalues]]
            
	\ \ \ \ \ \ \ \ \ \ \ ]\);

  
crop[polylist_,  funtest_]  :=  
        DeleteCases[Map[zap[#,  funtest]&,  polylist],  Null];
  
colorcrop[polylist_,  
          funtest_] := \(Module[{biglist}, 
          biglist\  = \ 
            Map[MapAt[zap[#, funtest]&, \ #, \ \(-1\)]&, \ polylist]; 
	
          DeleteCases[biglist, 
            x_ /; x[[\(-1\)]] == Null]]
          \);

Algorithms for getting normals and constructing lattice

getnormals[curve_,  algo_] := Module[{ ncurve  =  Length[curve],  tangent,  normal,  binormal, 
                 i},
	Switch[algo,
			x : {_, _, _}, \({
				\ \ 
                      tangent\ \ \  = \ Map[normalize, \ deriv[curve]], 

                      				\ \ 
                      normal\ \ \ \ \  = \ 
                        Map[normalize, \ \ Map[Cross[#, algo]&, \ tangent]], 
                      
			\ \ \ \ 
                      binormal\  = MapThread[Cross, \ {tangent, \ normal}]\ }
                      \), 
		      xytube,\({
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                      normal\  = \ Table[{0.6, \ 0.8, 0}, \ {i, \ ncurve}], 
                      binormal\  = \ 
                        Table[{0.8, \ \(-0.6\), \ 0}, \ {i, \ ncurve}]\ }\),
		  _,\({
\ \ \ \ \ \ \ \ \ \ \ \ 
                      tangent\ \ \  = \ Map[normalize, \ deriv[curve]], 

                      \ \ \ \ \ \ \ \ \ \ \ \ 
                      normal\ \ \ \ \  = \ normalvector[tangent], 

                      		\ \ \ \ \ \ 
                      binormal\  = MapThread[Cross, \ {tangent, \ normal}]\ }
                      \) 
		]; Transpose[{normal,  binormal}]
	];

getlattice[curve_,  shapefun_,  naround_,  nalong_,  helicity_, 
             algo_]:= Module[{costable, sintable, bothnormals, i, j, offset,  lattice}, 
   
		 bothnormals  =  getnormals[curve,  algo];
  
		\(If[helicity\  != \ 0, \ 
                  Module[{angle, \ newvecs}, \ 
                    angle[t_]\  := \ 
                      N[2\ Pi\ \((t - 1)\)\ helicity/\ \((nalong - 1)\)]; 

                    			bothnormals\  = \ 
                      Table[rotate[\ 
                          bothnormals[[i
                            ]]\ , \ \ angle[i]]\ , 
                        \ {i, \ 1, \ nalong}]
				]\ ]\);
  
   		costable  =  
                  Table[N[Cos[2 Pi j/naround]],  {j,  1,  naround}];
	       sintable  =  
                  Table[N[Sin[2 Pi j/naround]],  {j,  1,  naround}];
  	
	\(offset\  = \ 
                  If[NumberQ[shapefun], 
                    shapefun\ *\ 
                      Table[{costable[[j
                              ]], \ 
                            sintable[[j
                              ]]} . \ 
                          bothnormals[[i
                            ]], {j, \ 1, \ naround}, {i, \ 
                          1, \ nalong}], 			
	
                    Module[{radiustable}, \ 
                      radiustable\  = todata[shapefun, \ naround, \ nalong]; 
                      
					
                      Table[\ radiustable[[j, i
                            ]] {
                              costable[[j
                                ]], \ 
                              sintable[[j
                                ]]} . \ 
                            bothnormals[[i
                              ]], {j, \ 1, \ naround}, {i, 
                          \ 1, \ nalong}]]
					\ \ \ \ \ \ \ ]\);
  
  lattice =  Map[#  +  curve&, offset];
  AppendTo[lattice,  
                  lattice[[1]]]; lattice 
		];

Main functions gettube and getcolortube

tube::bad = "You have fed something weird to tube";

 polygondata[lattice_,  j_, i_] := 
         {lattice[[j, i]],  
          lattice[[j + 1, i]],  
          lattice[[j + 1,  i + 1]],  
          lattice[[j,  i + 1]]};

  
tosurf[x_] :=  
        If[Head[x]  ===  RGBColor  ||  Head[x]  ===  Hue  ||  
            Head[x]  ===  CMYKColor,  SurfaceColor[x],  x];
resolvecolors[directive_] := Map[tosurf,   Flatten[{directive}]];

gettube[curve_,  directive_,  shapefun_,  naround_,  range_,  
          helicity_,  algo_] := If[MatchQ[curve,  curvedata], 
		Module[{nalong  =  Length[curve],  polylist,  i,  j,  lattice, 
                   sheaf}, 
			lattice  =  
                    getlattice[curve,  shapefun,  naround,  nalong,  
                      helicity,  algo];
           sheaf        =  
                    Flatten[
                      Table[polygondata[ lattice,  j, i], 
                         {i,  1,  nalong - 1},  {j,  1,  naround}], 1]
                      ;
  
                      	
			polylist  =  
                    If[NumberQ[range[{0, 0, 0}]],   crop[sheaf,  range],  
                      sheaf];
  	
			polylist  =  Map[Polygon,  polylist];
			directiveplusreflectivity  =  
                    Prepend[resolvecolors[directive], 
                      SurfaceColor[Eggshell,  GrayLevel[1], 15]];
			Graphics3D[
                    FlattenAt[{directiveplusreflectivity,  polylist},  1]]
                    
			], Message[tube::bad]];

Main program getcolortube for variable color directives

(* This version of polygondata combines a list of graphics directives with a set of lattice points .  *)
colorpolygondata[directivetable_,  lattice_,  j_, i_]:= 
FlattenAt[{ 
              resolvecolors[
                directivetable[[j, i
                  ]]], {
                lattice[[j, i]],  
                lattice[[j + 1, i]],  
                lattice[[j + 1,  
                  i + 1]],  
                lattice[[j,  i + 1]]}}, 
             1];

getcolortube[curve_,  directive_,  shapefun_,  naround_,  
            range_,  helicity_,  algo_] := If[MatchQ[curve,  curvedata], 
		Module[{ nalong  =  Length[curve],  polylist,  i,  j,  
                    lattice,  sheaf,  colortable}, 
		lattice  =  
                      getlattice[curve,  shapefun,  naround,  nalong,  
                        helicity,  algo];
       colortable  = 
                      todata[directive,  naround,  nalong];			
	   \(sheaf\ \ \ \ \ \ \ \ \ \  = 	
                      Flatten[
                        Table[colorpolygondata[colortable, \ lattice, \ j, 
                            i], \ {i, \ 1, \ nalong - 1}, 
                          \ {j, \ 1, \ naround}], 1]\);
  	
			polylist  =  
                      If[NumberQ[range[{0, 0, 0}]],   
                        colorcrop[sheaf,  range],   sheaf];
  	
	polylist  =  
                      Map[ MapAt[Polygon,  #,  -1]&,  polylist]; 
	Graphics3D[
                      Join[{SurfaceColor[Eggshell,  GrayLevel[1], 15] },  
                        polylist]]
			], Message[tube::bad]];

User functions: tube, quicktube, checkedcurve

\(frenettube\  = \ 1; \ xytube\  = \ 2; 

    Options[tube] = {nalongtube -> 40, naroundtube -> 12, 
        tuberange -> All, tubehelicity -> 0, 
        tubealgorithm\  -> \ frenettube}; 
tubedefaults\  = \ Options[tube]; 
    

checkedcurve[curve_, \ n_]\  := \ \ 
      Which[\ MatchQ[curve, \ curvedata], curveinterp[curve, \ n], 

        \ \ \ \ \ \ \ VectorQ[N[curve[0]], \ NumberQ], todata[curve, n], 

        		True, \ Null]; 


    tube[curve_, \ directive_, \ shapefun_, \ opts___] := \ 
      gettube[checkedcurve[curve, 
          \(nalongtube /. {opts}\) /. Options[tube]\ ], directive, shapefun, 
        
		\ \ \(naroundtube\ \ \ \  /. {opts}\) /. Options[tube], 

        \ \ \ \ \ \ \ \ \ 
        \(tuberange\ \ \ \ \ \ \ \ \  /. {opts}\) /. Options[tube], 
	\ \ 	
        \(tubehelicity\ \ \  /. {opts}\) /. Options[tube], 

        \ \ \ \ \ \ \ \ \ \(tubealgorithm\  /. {opts}\) /. Options[tube]

        \ \ \ \ \ ]; 


    tube[curve_, \ directive_Symbol, \ shapefun_, \ opts___] := \ 
      getcolortube[
        checkedcurve[curve, \(nalongtube /. {opts}\) /. Options[tube]\ ], 
        directive, shapefun, 
		\ \ 
        \(naroundtube\ \ \ \  /. {opts}\) /. Options[tube], 

        \ \ \ \ \ \ \ \ \ 
        \(tuberange\ \ \ \ \ \ \ \ \  /. {opts}\) /. Options[tube], 
	\ \ 	
        \(tubehelicity\ \ \  /. {opts}\) /. Options[tube], 

        \ \ \ \ \ \ \ \ \ \(tubealgorithm\  /. {opts}\) /. Options[tube]

        \ \ \ \ \ ]; 

\ 
    quicktube[curve_] := 
      Show[tube[curve, Eggshell, \ 1], \ Boxed -> \ False]; \)

User functions: fluxlines

\(\(getlonglines[tubeobject_, \ nlonglines_, \ naround_, \ nalong_] := 

      Module[{polylist, \ toomanylines, \ i, \ j}, \ 
        polylist\  = \(tubeobject\  // First\) // Last; 
						\ 
        toomanylines\  = \ 
          polylist /. \ 
            Polygon[{a_, \ b_, \ c_, \ d_}]\ \  -> \ Line[{a, d}]; 

        										
        Flatten[\ 
          Table[toomanylines
              [[\((i\ \ naround + 1)\)\  + \ 
                j\ naround/\ nlonglines]], 
            \ {i, 0, \ nalong - 2}, {j, \ 0, \ \ nlonglines - 1}], 1]
	]; 

    
getshortlines[tubeobject_, \ nshortlines_, \ \ naround_, \ nalong_] := 
      
	Module[{polylist, \ toomanylines, \ i, \ j}, \ 
        polylist\  = \(tubeobject // First\) // Last\ ; 
						\ 
        toomanylines\  = \ 
          polylist /. \ 
            Polygon[{a_, \ b_, \ c_, \ d_}]\ \  -> \ Line[{a, b}]; 

        										
        Flatten[\ 
          Table[toomanylines
              [[i\ \  + \ 
                j\ \((nalong - 1)\)\ 
                  naround/nshortlines]], 
            \ {i, \ 1, \ naround}, {j, \ 0, \ nshortlines - 1}], 1]
	]; 

    										

    getflux[curve_, \ radius_, \ nlonglines_, \ nshortlines_, \ 
        fluxthickness_, \ fluxcolor_] := 
      Module[{newtube, \ lines, naround, nalong}, 
				
        naround\  = naroundtube /. Options[tube, \ naroundtube]\ ; 
        nalong\  = nalongtube /. \ Options[tube, \ nalongtube]; \ 
			\ 
        naround\  = \ Ceiling[naround/nlonglines]\ nlonglines; 
			
        nalong\  = \ 1 + Ceiling[\((nalong - 1)\)/nshortlines]\ nshortlines; 
        
		newtube = 
          tube[curve, \ White, \ radius, \ naroundtube\  -> naround, \ 
            nalongtube -> nalong]; 

        							{
          Graphics3D[
            \ {Thickness[fluxthickness], fluxcolor, \ 
              getlonglines[newtube, \ \ \ nlonglines, \ naround, \ nalong]
                \ }], 
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Graphics3D[
            \ {Thickness[fluxthickness], fluxcolor, 
              getshortlines[newtube, nshortlines, \ naround, \ nalong]\ }\ ]

          		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }]; 			


    Options[fluxlines] = {fluxthickness -> 0.005, 
        fluxcolor -> \ \ GrayLevel[0]}; 


    fluxlines[curve1_, \ radius_, \ nlonglines_, \ nshortlines_, \ 
        opts___] := 
\ \ \ \ \ \ 
      getflux[\ curve1, \ radius, \ nlonglines, \ nshortlines, \ 
		
        \(fluxthickness /. {opts}\) /. Options[fluxlines], 
		
        \(fluxcolor /. {opts}\) /. Options[fluxlines]\ ]; 			\ \)\)

Main function getarrow

\(\(\(getarrow[pp_, \ qq_, \ normal_, \ thickness_]\  := \ 
      Module[{a, b, c, d, e, f, g, nhat, len, \ thick, \ sticklength, \ 
          trilength, hwidth, 
			top, \ bot, \ bottom, \ side1, side2, \ 
          wing1, \ wing2, \ back1, \ back2}, \ 
		w\ \ \ \ \  = qq - pp; \ 
        
		len\  = \ \ N[Sqrt[w . w]]; 
		\ w\ \ \ \  = \ N[w/len]; 

        		\ nhat\  = \ 
          N[Cross[w, \ Cross[normal, \ w]]/Sqrt[normal . normal]]; \ 

        	\ \ \ v\  = \ Cross[nhat, \ w]; 
		sticklength\  = \ len/2; 

        		trilength\  = \ len\  - \ sticklength; 
		
        hwidth\  = \ sticklength/6; 
		thick\ \  = \ thickness\ len; 

        		 (*\ top\ \ of\ arrow\ *) 
		a\  = \ pp\  + \ hwidth\ v; 
        b\  = \ pp\  - \ hwidth\ v; c\  = \ b\  + \ sticklength\ w; 
        g\  = \ a\  + \ sticklength\ w; 
		d = \ c\  - \ \ \ hwidth\ v; 
        e\  = \ qq; 
		f = \ g\  + hwidth\ v; 

        		 (*\ bottom\ of\ arrow\ *) 
	\ \ \ 
        bot[{x_, y_, z_}]\  := \ {x, y, z}\  - \ thick\ nhat; 
		
        top\  = \ {a, b, c, g}; arrowtop\  = \ {d, \ e, f}; \ 
        bottom\  = \ Map[bot, \ top]; \ 
        arrowbottom\  = \ Map[bot, \ arrowtop]; \ 
        end\  = \ {a, b, bot[b], \ bot[a]}; 
		
        side1\ \  = \ {b, c, \ bot[c], \ bot[b]}; \ 	
        side2\ \  = \ {g, a, \ bot[a], \ bot[g]}; 
		
        wing1\ \  = \ {d, e, \ bot[e], \ bot[d]}; \ 	
        wing2\ \  = \ {e, f, \ bot[f], \ bot[e]}; \ 
		
        back1\ \  = \ {d, c, \ bot[c], \ bot[d]}; \ 	
        back2\ \  = \ {g, f, \ bot[f], \ bot[g]}; 
		\ 
        polly\  = 
          \ {\ top, \ arrowtop, \ bottom, \ arrowbottom, \ end, \ side1, \ 
            side2, \ wing1, \ wing2, \ back1, \ back2}; 
		
        Map[Polygon, polly]]; \)
		\)\)

User functions: arrow, sidearrow

\(Options[arrow]\  = 
      \ {arrowcolor\  -> \ Red, arrowthickness\ \  -> 0.1}; 
\ 
    showarrow[p_, q_, n_, color_, \ thickness_]\  := 
		\ 
      Graphics3D[{EdgeForm[], SurfaceColor[color, GrayLevel[1.0], \ 15], 
          getarrow[p, q, n, thickness]}]; 


    arrow[p_, q_, n_, opts___]\  := \ 
      showarrow[p, q, n, \(arrowcolor\  /. {opts}\) /. Options[arrow], \ 

        		\(arrowthickness\ \ \ \  /. {opts}\) /. Options[arrow]]; \)

\(Options[sidearrow]\  = 
      \ {arrowcolor\  -> Red, arrowthickness\ \  -> 0.1}; 

    sidearrow[fun_, \ t_, \ eps_, \ radius_, \ theta_, opts___]\  := \ 
      Module[{ppp, \ qqq, \ norm, \ a, \ b, \ qp, \ woof, view, off}, 
		
        a\  = \ fun[t]; \ 
		b\  = \ fun[Mod[t + eps, \ 1]]; 
		
        len\  = \ Sqrt[\((b - a)\) . \((b - a)\)]; 
		
        qp\  = \ \((b - a)\)/len; 
		
        woof\  = \ Options[Graphics3D, \ ViewPoint]; 
		
        view\  = \ ViewPoint\  /. \ woof; 
		norm\  = \ normalize[view]; 

        		norm\  = Cross[qp, \ Cross[norm, \ qp]]; 
		
        binorm\  = \ Cross[norm, \ qp]; 	\ \ \ 
        off\  = \((radius\  + \ 0.1\ len)\)*\ 
            \(rotate[{norm, \ binorm}, \ theta]\)[[1
              ]]; 
		ppp\  = \ a\  + \ off; 
		
        qqq\  = \ b\  + \ off; 
		
        showarrow[ppp, qqq, off, \ 
          \(arrowcolor\  /. {opts}\) /. Options[sidearrow], \ 
		
          \(arrowthickness\ \ \ \  /. {opts}\) /. Options[sidearrow]]]; \)

User Functions: tress

\(boundingbox[braid_] := \ 
      Module[{pts, \ nz, \ nc, \ i}, nc = Length[braid]; 
			
        nz = Length[braid[[1]]]; \ \ 

        		pts = 
          Transpose[
            Join[Table[braid[[i, 1]], 
                \ {i, nc}], 
              Table[braid[[i, nz]], 
                \ {i, nc}]]]; 
		{Map[Min, \ pts], \ Map[Max, pts]}]; 


    pedestal[box_, \ \ thick_, \ margin_] := \ \ 
      Module[{left, \ right, \ front, \ back, \ bot, \ top}, \ 
			
        left = \ box[[1, 1]] - margin; \ 
        right = box[[2, 1]] + margin; \ 
        bot = box[[1, 3]]; 
        tresswidth = right; 
		
        front = \ box[[1, 2]] - margin; 
        \ back = box[[2, 2]] + margin; \ 
        top\  = \ box[[2, 3]]; \ 

        \ \ \ \ \ 	{
          Graphics3D[
            Cuboid[{left, \ front, bot\  - \ thick/2}, 
              \ {right, \ back, \ bot + thick/2}]], 
          Graphics3D[
            Cuboid[{left, \ front, top\  - \ thick/2}, 
              \ {right, \ back, top + \ thick/2}]]}]; \)

\(\(braidtubes[braid_, \ colors_, \ radius_, \ braidrange_] := 
      Table[
		
        tube[braid[[i]], 
          	{EdgeForm[], \ 
            colors[[Mod[i - 1, \ \ Length[colors]] + 1
                ]]}, radius, \ tuberange\  -> \ braidrange]
          \ , 
		\ \ \ {i, \ 1, \ Length[braid]}]; \)\)

\(\(nicecolors = {RGBColor[1, \ 0, \ 0], RGBColor[0, 1, \ 1], 
        RGBColor[1, \ 1, \ 0], RGBColor[0, 1, \ 0], 
		
        RGBColor[1, \ 0, \ 1], \ RGBColor[0, \ 0, \ 1], \ 
        RGBColor[1, \ 0.5, \ 0.5], RGBColor[1, \ 0.5, \ 0], 
        RGBColor[0, \ 1, \ 0.5], RGBColor[1, \ 0, \ 0.5], 
        RGBColor[0.5, \ 0, 1], RGBColor[0.25, \ 0.5, \ 0.5]}; \)\)

\(Options[tress]\  = 
      \ {braidcolors\  -> nicecolors, \ pedestal\  -> \ {4, \ 0.75}}; 


    tress[curvelist_, radius_, \ opts___]\  := 
      Module[{colors, braidlist, ped, \ margin, \ thick, \ box, braidrange, 
          i}, \ 
	\ 
        colors\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  = \ 
          \(braidcolors\  /. \ {opts}\)\  /. \ Options[tress]; 
	\ 
        braidlist = 
          Table[checkedcurve[
              curvelist[[i]], 
              nalongtube /. Options[tube]], {i, 1, Length[curvelist]}]; 
	\ 
        ped\  = \ \(pedestal\  /. \ {opts}\)\  /. \ Options[tress]; 
	\ 
        ped\  = \ ped\ *\ If[NumberQ[radius], \ radius, \ radius[0, 0]]; 

        		margin\  = \ ped[[1]]\ ; 

        		thick\ \ \ \ \  = \ 
          ped[[2]]; 
		
        box = boundingbox[braidlist]; \ 
        tresswidth\  = \ 
          box\ [[2, 1]]\  - \ 
            box\ [[1, 1]]; 
		
        braidrange[vec_]\  := \ 
          If[vec[[3]] >= \ 
                box\ [[1, 3]]\  && \ 
              vec[[3]] <= \ 
                box\ [[2, 3]], \ 1, \ 
            \(-1\)]; 
		
        graphicslist = \ 
          Flatten[{pedestal[box, \ thick, \ margin], 
              braidtubes[braidlist, colors, \ radius, \ braidrange]}]; 
		
        Graphics3D[
          Table[graphicslist[[i, 1]], 
            \ {i, \ 1, \ Length[graphicslist]}], \ PlotRange -> All\ ]
	]; 
    \)

End[ ];  EndPackage[ ];  On[General::spell1]

Tubafunctions 
(June 1999) by M A Berger, Mathematics, University College London

BeginPackage["Topofunctions`"];   Needs["Topographics`"];  
    Off[General::spell1];

Usage Messages

arch

arch::usage = 
      "arch[{x1, y1}, {x2, y2}, t] gives a halfcircle in a verticle plane \
from {x1, y1,0} to {x2, y2, 0}. For more complicated loops above z=0, \
consider drawing a full circle or ellipse, then cutting off the part below \
z=0 by creating a function range[z_] = z, then setting tuberange -> range for \
the tube graphics routine.";

artin

\(\(artin::usage\  = \ 
      \*""artin[braidword, nstrings]creates a braid  with nstrings curves (nstrings lists of points) corresponding to braidword; for example artin[{1,-2,3}, 4] has four strings and artin word (Subscript[σ, 1])((Subscript[σ, 2])^(-1))(Subscript[σ, 3]). Each element has height 1, i.e. in the example z goes from 0 to 3.
    artin[braidword, nstrings, deltaz]: each element has height deltaz.
		""; \)\)

scaledbraid

scaledbraid::usage = "
		scaledbraid[braid, xyscale, n] first converts braid to a set of lists of \
length n. 
		It then chooses z values so that the ith point in the list for each curve \
is at z=i-1. ";

readcurves, writecurves

readcurves::usage = 
"readcurves[\"tiger.dat\", ncurves]  will read a data file called tiger.dat . This file should contain nz*ncurves lines, where each line has three floating point numbers ({x,y,z} coordinates). The top nz lines
describe curve 1, and so on. readcurves[\"tiger.dat\"] assumes that the first line consists of the two numbers ncurves, nz.";
writecurves::usage = 
"writecurves[\"tiger.dat\", curves] will write out the coordinates in the list curves to the file 
tiger.dat. The first lines gives  the number of curves and the number of points per curve nz. After this will each line gives a point as three floating point numbers ({x,y,z} coordinates). The top nz points describe curve 1, and so on.";

torusknot

\(\(torusknot::usage\  = \ 
\ 
      "torusknot[p,q, majorradius,minorradius, t] is a function of t which gives a torusknot curve. The curve is periodic in t with period 1. It is imbedded in a torus described by majorradius and minorradius, and goes p times the long way around and q times the short way around.
		"; \)\)

Code

Begin["`Private`"];

arch

\(\(arch[{x1_, \ y1_}, \ {x2_, \ y2_}, \ t_]\  := \ 
      Module[{u, \ uhat, \ mu, \ distance, \ \ midpoint}, 
		
        u\  = \ {x2 - x1, \ y2 - y1, \ 0}; 
		distance\  = \ Sqrt[u . u]; 
        
		uhat\  = \ u/distance; 
		
        midpoint\  = \ \(({x1, \ y1, \ 0}\  + \ {x2, \ y2, \ 0})\)/2; 
		
        mu\  = \ Pi\ \((t - 1/2)\); 
		
        N[midpoint\  + \ 
            distance\ \((Sin[mu]\ uhat\  + \ Cos[mu]\ {0, 0, 1})\)/2]
]; 
    \)\)

artin

Clear[artin]

\(artin[braidword_, \ nstrings_, \ deltaz_]\  := \ 
      If[VectorQ[braidword, \ IntegerQ], \ 
        Module[{nsigmas\  = \ Length[braidword], \ stringposn, \ sig, \ mu, \ 
            h, \ braid, \ j, \ iz}, \ stringposn\  = \ Range[nstrings]; \ \ 
          braid\  = \ Table[{{j, 0, 0}}, \ {j, nstrings}]; \ 
			
          Do[{
		
              sig\  = Abs[
                  braidword[[iz]]], 
              mu\  = \ Sign[
                  braidword[[iz]]], 

              					
              If[0 < sig < nstrings, {
					
                  h = \ \((iz - 0.5)\) deltaz, \ 
                  braid = \ 
                    braidappend[braid, \ 
						
                      simplebraid[nstrings, \ stringposn, \ sig, \ mu, 

                        							\ h]], 
                  stringposn = flip[stringposn, \ sig], 
                  If[iz < nsigmas, \ 
                    If[braidword[[iz
                          ]] == 
                        braidword
                          [[iz + 1]], 

                      									
                      braid = \ 
                        braidappend[braid, \ \ 
                          Table[{
                              stringposn[[j
                                ]], 0, iz\ deltaz}, 
                            \ {j, nstrings}]]
								]]
\ }]

              				\ \ \ }, {iz, \ nsigmas}]; 
			
          If[nsigmas\  == 0, h = 1, \ h = nsigmas*deltaz]; \ 
          braidappend[braid, \ \ 
            Table[{stringposn[[j]], 0, 
                h}, \ {j, nstrings}]]
		], \ 
        Print["Bad input for artin function\!"]]; 

    artin[braidword_, \ nstrings_] := \ artin[braidword, \ nstrings, \ 1]; 
    \)

\(flip[perm_, \ s_]\  := 
      If[0 < s < Length[perm], \ \ 
        Module[{new, \ leftstring, \ rightstring}, \ 
			
          leftstring\ \  = 
            \(Position[perm, \ s]\)[[1, 1
              ]]; 
			
          rightstring = 
            \(Position[perm, \ s + 1]\)[[1, 1
              ]]; 
			
          new\  = \ ReplacePart[perm, \ s + 1, \ leftstring]; 
			
          ReplacePart[new, \ s, \ rightstring]]
	]; \ 

    braidcat[botbraid_, \ topbraid_]\  := \ 
      MapThread[Join, \ {botbraid, \ topbraid}]; 

    braidappend[braid_, \ points_]\  := \ 
      MapThread[Append, \ {braid, \ points}]; \)

\(\(\(simplebraid[nstrings_, perm_, \ s_, \ mu_, h_] := 
      If[0 < s < Length[perm], \ \ 
        Module[{nextlevel, \ leftstring, \ rightstring, \ j}, \ 
			
          leftstring\ \  = 
            \(Position[perm, \ s]\)[[1, 1
              ]]; 
			
          rightstring = 
            \(Position[perm, \ s + 1]\)[[1, 1
              ]]; 
		
          nextlevel\ \ \ \ \  = \ 
            Table[{perm[[j]], 0, h}, 
              \ {j, nstrings}]; 
		
          nextlevel\ \ \ \ \  = \ 
            ReplacePart[nextlevel, {s + 0.5, \ \(-0.5\) mu, \ h}\ , 
              leftstring]; \ 	
		
          ReplacePart[nextlevel, \ {s + 0.5, \ \ 0.5  mu, \ h}, rightstring]
            \ ], 
		]; \)
		\)\)

scaledbraid

scaledbraid[braid_, xyscale_, n_] :=  Module[{ncurves = Length[braid], tempbraid, jcurve, ipoint}, tempbraid = Map[checkedcurve[#, n]&, braid]; 
Table[{xyscale tempbraid[[jcurve, ipoint, 1]],
       xyscale tempbraid[[jcurve, ipoint, 2]], ipoint-1}, 
			{jcurve, ncurves}, {ipoint,n}]];

readcurves, writecurves

readcurves[file_, ncurves_] := Module[{nz, intnz, data}, data = ReadList[file, Number, RecordLists -> True]; nz = Length[data]/ncurves; intnz = IntegerPart[nz]; If[Not[nz==intnz], {nz=intnz, Print["Warning! Length of data set not an integer multiple of number of curves."]}];  Table[Take[data, {(i-1)nz + 1, i*nz}], {i, ncurves}] ];

readcurves[file_] := Module[{ncurves, nz, data, i}, data = ReadList[file, Number, RecordLists -> True]; 
		ncurves = data[[1, 1]];
		nz = data[[1, 2]]; 
		ntest = (Length[data]-1)/ncurves; If[Not[nz+1==ntest], Print["Warning! Length of data set not equal to (nz+1)*ncurves from first line of file! "]];  Table[Take[data, {(i-1)(nz+1) + 2, i*(nz+1) + 1}], {i, ncurves}] ];

writecurves[file_,curves_]:= Module[{t, ncurves=Length[curves], 
nz = Length[ curves[[1]]]-1, icurve, ipoint}, t = OpenWrite[file, FormatType -> OutputForm];
		Write[t, ncurves,"  ", nz];
Do[Do[Write[t,curves[[icurve, ipoint, 1]]," ",
curves[[icurve, ipoint, 2]]," ",
curves[[icurve, ipoint, 3]]], {ipoint, 1, nz+1}], {icurve,ncurves}];
Close[t];
];

torusknot

\(\(torusknot[p_, \ q_, \ majorradius_, \ minorradius_, \ t_]\  := \ 
      Module[{ax, \ az, \ circ, \ offset}, 
		ax\  = \ 2\ Pi\ p\ t; \ 

        		az\  = \ 2\ Pi\ q\ t; 
		
        circ\  = \ majorradius\ {Cos[ax], \ Sin[ax], \ 0}; 
			
        offset\  = \ 
          minorradius 
            \(({Cos[ax] Cos[az], Sin[ax] Cos[az], \ \(-Sin[az]\)})\); \ 
		
        circ\  + \ offset]; \)\)

End[ ];  EndPackage[ ];  On[General::spell1];