From: phr@soda.berkeley.edu (Paul Rubin)
Newsgroups: sci.math.symbolic
Subject: cute topology demo in Mathematica
Date: 25 Dec 1992 20:54:40 GMT
Organization: University of California, Berkeley
NNTP-Posting-Host: soda.berkeley.edu

Here's something I wrote for a geometry class last spring.
It shows how you can continuously deform two-holed donut into
one where the holes are linked together, by generating a series
of still pictures showing different stages of the transformation.
(The first frame is a "before and after" picture).  I wanted to
expand it to make an animation but never got around to upgrading
to a Windows version of Mathematica, so I'm posting it now.

Note: it takes a couple of minutes to generate each frame on a 486/33
(MSDOS Mathematica 2.0), so be patient.  Also, you may need to change
the hardcopy function (3rd-to-last line in the program) to get
hardcopy on non-MSDOS systems.  To get screen display, just
delete that line.

================================================================
(* torus.m -- Paul Rubin (phr@soda.berkeley.edu), April-May 1992 
   Copyright (c) 1992, Paul Rubin.  Permission granted to copy and
   distribute per terms of GNU General Public License.  To distribute
   with different permissions, ask me first. *)

(* Acknowledgement: I thank Profs. Silvio Levy and Bill Thurston
   for helpful advice and encouragement for this project. *)

norm[vec_] := Sqrt [vec . vec]
normalize[vec_] := vec / norm[vec]

linearCombination[vecs_,weights_]:=Inner[Times,weights,vecs,Plus]

(* This generates a parametrized tube of radius `r' around a curve.
   Explanation: D[curve,t] is a vector giving the direction of the curve;
   NullSpace[D[curve,t]] is a pair of vectors that are perpendicular
   to the curve and to each other.  By normalizing these and generating
   the equation for a circle around the curve in the plane they form
   for any given t, we get the equation for the tube around the curve.

   We simplify this equation and make the substitution (a^b)^c == a^(bc)
   in order to change expressions in it like Sqrt[x^2] to x.
   This gets rid of some singularities in in the equations this
   generates for tori which keep the tori from getting drawn
   correctly.  It's surprising that Simplify doesn't make that
   substitution already! *)

tube[curve_] :=
     (Simplify [curve + 
       		r * linearCombination [Map [normalize,
						NullSpace [{D [curve, t]}]],
					{Cos[u], Sin[u]}]]
    /. (a_^b_)^c_ -> a^(b c))

(* draw1d draws the specified curves in 3-space.  This mainly exists
   to make sure the shapes of the curves are right, since it displays
   the curves much faster than draw2d does. *)
   
draw1d[curves__]:=ParametricPlot3D[Evaluate[{curves}],{t,0,2Pi},Boxed->False,
				PlotRange->All,Axes->False]

(* draw2d draws the 2-dimensional surface made by putting a tube
   around the specified curves.  It's quite slow. *)
draw2d[curves__]:=ParametricPlot3D[Evaluate[Map[tube,{curves}]],
		{t,0,2Pi},{u,0,2Pi},Boxed->False,PlotRange->All,
					Axes->False]

(* Set up shapes for some of the actual pictures in the
	torus-transformation sequence *)
r=1/4
circle = {Cos[t],Sin[t],0}
vcircle = {{1,0,0},{0,0,-1},{0,1,0}} . circle (* vertically oriented circle *)
cyl:={t/(2Pi),0,0}	(* unit straight line along x axis *)

circle2 = circle + {3, 0, 0}
vcircle2 = vcircle + {3, 0, 0}
cyl2:=cyl + {1,0,0}

vcircle3 = vcircle + {2, 0, 0}  (* to draw with no cylinder *)
vcircle4 = vcircle + {1, 0, 0}  (* draw with inward cylinder *)

cyl5 = {{0,-1,0},{0,1,0},{-1,0,0}} . cyl + {1,0,0} (* downward cylinder *)

(* define a "draw" function so we can easily change between 1d and 2d *)
(* draw = draw1d *)
draw = draw2d

(* xarch generates parametrization for an arch that goes through
   specified points.  Probably not a good idea to use this for
   much more than 3 points, as would get high-degree polynomials.
   What's really needed is a way to generate this arch as a bunch
   of cubic patches.  Interpolate[] does this, but returns an
   object that cannot be differentiated by D[]  :-(  *)
xarch[pts__] :=
    Table [InterpolatingPolynomial [
    		  Table [{2Pi (k-1)/(Length[{pts}]-1), {pts}[[k,i]]},
		  	 {k, 1, Length[{pts}]}], t],
	   {i, 1, Length[{pts}[[1]]]}]

arch1 := xarch[{0,-1,0},{3/2,-10,0},{3,-1,0}]
arch1a := xarch[{0,-1,0},{3/2,-5/2,0},{3,-1,0}]

(* now make expressions for the actual pictures in the sequence *)
(* d0 is the starting figure (two unlinked tori connected by a long arch).
   d0a is the same figure with a shorter arch.
   d1 is two unlinked tori in the X-Y plane connected by a straight cylinder.
   d2 is d1 with one torus rotated into the Z-plane
   d3 has the two tori touching
   d4 pushes them through each other, connecting them with a cylinder
      that is drawn in the "intersection"
   d5 rotates one torus 90 degrees so the cylinder is no longer in
      the "intersection"
   d6 rotates it a further 90 degrees so the cylinder is in the X-Y plane
      and completely outside one of the tori
   d7 bends the cylinder so it becomes an arch that is outside both tori
   d8 lengthens the arch so it looks like the one from d0. 
   dstart has the starting and ending figures alongside each other
      in one drawing. *)

d0 := draw[circle,circle2,arch1]
d0a := draw[circle,circle2,arch1a]
d1:=draw[circle,circle2,cyl2] 
d2:=draw[circle,vcircle2,cyl2]
d3:=draw[circle,vcircle3]
d4:=draw[circle,vcircle4,cyl]
d5:=draw[circle,vcircle4,cyl5]
d6:=draw[circle,vcircle4,cyl2]

arch2 := xarch[{0,-1,0},{1,-2,0},{2,0,0}]
arch3 := xarch[{0,-1,0},{1,-10,0},{2,0,0}]
d7 := draw[circle,vcircle4,arch2]
d8 := draw[circle,vcircle4,arch3]

dx={-7,0,0}
dstart:=draw[circle+dx,circle2+dx,arch1+dx,circle,vcircle4,arch3]

(* to print all frames, uncomment the following code.  Watch out because
   it takes several minutes to do each frame. *)


	$DisplayFunction=Hardcopy
	draw=draw2d
	dstart; d0; d0a; d1; d2; d3; d4; d5; d6; d7; d8