home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #31 / NN_1992_31.iso / spool / sci / math / symbolic / 3336 < prev    next >
Encoding:
Internet Message Format  |  1992-12-25  |  5.9 KB
  1. Path: sparky!uunet!spool.mu.edu!agate!soda.berkeley.edu!phr
  2. From: phr@soda.berkeley.edu (Paul Rubin)
  3. Newsgroups: sci.math.symbolic
  4. Subject: cute topology demo in Mathematica
  5. Date: 25 Dec 1992 20:54:40 GMT
  6. Organization: University of California, Berkeley
  7. Lines: 138
  8. Message-ID: <1hfsegINN6ja@agate.berkeley.edu>
  9. NNTP-Posting-Host: soda.berkeley.edu
  10.  
  11. Here's something I wrote for a geometry class last spring.
  12. It shows how you can continuously deform two-holed donut into
  13. one where the holes are linked together, by generating a series
  14. of still pictures showing different stages of the transformation.
  15. (The first frame is a "before and after" picture).  I wanted to
  16. expand it to make an animation but never got around to upgrading
  17. to a Windows version of Mathematica, so I'm posting it now.
  18.  
  19. Note: it takes a couple of minutes to generate each frame on a 486/33
  20. (MSDOS Mathematica 2.0), so be patient.  Also, you may need to change
  21. the hardcopy function (3rd-to-last line in the program) to get
  22. hardcopy on non-MSDOS systems.  To get screen display, just
  23. delete that line.
  24.  
  25. ================================================================
  26. (* torus.m -- Paul Rubin (phr@soda.berkeley.edu), April-May 1992 
  27.    Copyright (c) 1992, Paul Rubin.  Permission granted to copy and
  28.    distribute per terms of GNU General Public License.  To distribute
  29.    with different permissions, ask me first. *)
  30.  
  31. (* Acknowledgement: I thank Profs. Silvio Levy and Bill Thurston
  32.    for helpful advice and encouragement for this project. *)
  33.  
  34. norm[vec_] := Sqrt [vec . vec]
  35. normalize[vec_] := vec / norm[vec]
  36.  
  37. linearCombination[vecs_,weights_]:=Inner[Times,weights,vecs,Plus]
  38.  
  39. (* This generates a parametrized tube of radius `r' around a curve.
  40.    Explanation: D[curve,t] is a vector giving the direction of the curve;
  41.    NullSpace[D[curve,t]] is a pair of vectors that are perpendicular
  42.    to the curve and to each other.  By normalizing these and generating
  43.    the equation for a circle around the curve in the plane they form
  44.    for any given t, we get the equation for the tube around the curve.
  45.  
  46.    We simplify this equation and make the substitution (a^b)^c == a^(bc)
  47.    in order to change expressions in it like Sqrt[x^2] to x.
  48.    This gets rid of some singularities in in the equations this
  49.    generates for tori which keep the tori from getting drawn
  50.    correctly.  It's surprising that Simplify doesn't make that
  51.    substitution already! *)
  52.  
  53. tube[curve_] :=
  54.      (Simplify [curve + 
  55.                r * linearCombination [Map [normalize,
  56.                         NullSpace [{D [curve, t]}]],
  57.                     {Cos[u], Sin[u]}]]
  58.     /. (a_^b_)^c_ -> a^(b c))
  59.  
  60. (* draw1d draws the specified curves in 3-space.  This mainly exists
  61.    to make sure the shapes of the curves are right, since it displays
  62.    the curves much faster than draw2d does. *)
  63.    
  64. draw1d[curves__]:=ParametricPlot3D[Evaluate[{curves}],{t,0,2Pi},Boxed->False,
  65.                 PlotRange->All,Axes->False]
  66.  
  67. (* draw2d draws the 2-dimensional surface made by putting a tube
  68.    around the specified curves.  It's quite slow. *)
  69. draw2d[curves__]:=ParametricPlot3D[Evaluate[Map[tube,{curves}]],
  70.         {t,0,2Pi},{u,0,2Pi},Boxed->False,PlotRange->All,
  71.                     Axes->False]
  72.  
  73. (* Set up shapes for some of the actual pictures in the
  74.     torus-transformation sequence *)
  75. r=1/4
  76. circle = {Cos[t],Sin[t],0}
  77. vcircle = {{1,0,0},{0,0,-1},{0,1,0}} . circle (* vertically oriented circle *)
  78. cyl:={t/(2Pi),0,0}    (* unit straight line along x axis *)
  79.  
  80. circle2 = circle + {3, 0, 0}
  81. vcircle2 = vcircle + {3, 0, 0}
  82. cyl2:=cyl + {1,0,0}
  83.  
  84. vcircle3 = vcircle + {2, 0, 0}  (* to draw with no cylinder *)
  85. vcircle4 = vcircle + {1, 0, 0}  (* draw with inward cylinder *)
  86.  
  87. cyl5 = {{0,-1,0},{0,1,0},{-1,0,0}} . cyl + {1,0,0} (* downward cylinder *)
  88.  
  89. (* define a "draw" function so we can easily change between 1d and 2d *)
  90. (* draw = draw1d *)
  91. draw = draw2d
  92.  
  93. (* xarch generates parametrization for an arch that goes through
  94.    specified points.  Probably not a good idea to use this for
  95.    much more than 3 points, as would get high-degree polynomials.
  96.    What's really needed is a way to generate this arch as a bunch
  97.    of cubic patches.  Interpolate[] does this, but returns an
  98.    object that cannot be differentiated by D[]  :-(  *)
  99. xarch[pts__] :=
  100.     Table [InterpolatingPolynomial [
  101.               Table [{2Pi (k-1)/(Length[{pts}]-1), {pts}[[k,i]]},
  102.                {k, 1, Length[{pts}]}], t],
  103.        {i, 1, Length[{pts}[[1]]]}]
  104.  
  105. arch1 := xarch[{0,-1,0},{3/2,-10,0},{3,-1,0}]
  106. arch1a := xarch[{0,-1,0},{3/2,-5/2,0},{3,-1,0}]
  107.  
  108. (* now make expressions for the actual pictures in the sequence *)
  109. (* d0 is the starting figure (two unlinked tori connected by a long arch).
  110.    d0a is the same figure with a shorter arch.
  111.    d1 is two unlinked tori in the X-Y plane connected by a straight cylinder.
  112.    d2 is d1 with one torus rotated into the Z-plane
  113.    d3 has the two tori touching
  114.    d4 pushes them through each other, connecting them with a cylinder
  115.       that is drawn in the "intersection"
  116.    d5 rotates one torus 90 degrees so the cylinder is no longer in
  117.       the "intersection"
  118.    d6 rotates it a further 90 degrees so the cylinder is in the X-Y plane
  119.       and completely outside one of the tori
  120.    d7 bends the cylinder so it becomes an arch that is outside both tori
  121.    d8 lengthens the arch so it looks like the one from d0. 
  122.    dstart has the starting and ending figures alongside each other
  123.       in one drawing. *)
  124.  
  125. d0 := draw[circle,circle2,arch1]
  126. d0a := draw[circle,circle2,arch1a]
  127. d1:=draw[circle,circle2,cyl2] 
  128. d2:=draw[circle,vcircle2,cyl2]
  129. d3:=draw[circle,vcircle3]
  130. d4:=draw[circle,vcircle4,cyl]
  131. d5:=draw[circle,vcircle4,cyl5]
  132. d6:=draw[circle,vcircle4,cyl2]
  133.  
  134. arch2 := xarch[{0,-1,0},{1,-2,0},{2,0,0}]
  135. arch3 := xarch[{0,-1,0},{1,-10,0},{2,0,0}]
  136. d7 := draw[circle,vcircle4,arch2]
  137. d8 := draw[circle,vcircle4,arch3]
  138.  
  139. dx={-7,0,0}
  140. dstart:=draw[circle+dx,circle2+dx,arch1+dx,circle,vcircle4,arch3]
  141.  
  142. (* to print all frames, uncomment the following code.  Watch out because
  143.    it takes several minutes to do each frame. *)
  144.  
  145.  
  146.     $DisplayFunction=Hardcopy
  147.     draw=draw2d
  148.     dstart; d0; d0a; d1; d2; d3; d4; d5; d6; d7; d8
  149.