用 Mathematica 玩转环面
环面及其变体
要玩转环面，先要构造出环面，然后才可以谈其它。本节将介绍如何从环面出发，用数学公式让它发生各种形变，以及如何变化参数，生成动画。
01
构造环面
我们都很熟悉圆的参数方程，比如对一个半径为 R 的圆心位于原点的圆，圆上的点的坐标 (x, y) 满足以下方程：
这个方程虽然准确，却不容易根据它直接画出图形。为了能构造性的生成曲线，采用所谓“参数方程”比较方便，也就是把 x、y 当作另一个参数 t 的函数。上述圆方程的一种可能的参数方程是：
参数 t 可以看作圆上的点和圆心所连成的直线和 x 轴的夹角，t 从 0 取到 2π，就完成了绕圆一周，也就画出了一个圆。有了圆的参数方程，就不难构造出环面的参数方程。直观想象一下，把一个圆 B 放在 xy 平面上，再有一个垂直于 xy 平面的立起来的圆 A，A 的圆心 O 始终落在 B 上，朝向是 B 上 O 点处的切线方向。只要圆 A 绕 B 一周，就形成了一个环面。如下图所示，红色的就是圆 A，黑色的就是圆 B。
为不失一般性，不妨假设红圆的半径是 r，其圆心在黑圆上的点 (R Cos[u], R Sin[u], 0) 处。该怎么才能画出红圆呢？我们需要知道此时红圆的"局部标架"，也就是上图中绿色部分的三个箭头代表的三个单位向量。垂直向上的箭头很好办，就是 (0, 0, 1)；垂直于红圆，也就是往里的那个箭头不影响绘制红圆可以不管；而剩下的朝右的那个箭头和圆心是共线的，所以就是 (Cos[u], Sin[u], 0)。于是假设红圆的参数为 v，可计算其参数方程如下：
圆心坐标+上箭头向量 r sin(v)+右箭头向量 r cos(v)
用 Mathematica 计算化简
于是我们定义环面参数方程函数如下。说明一下这个函数还有个可选的设置 L 用于控制环面在 z 轴方向的拉伸程度。默认就是 1，放在这里主要是为了绘制某种曲面时能复用这部分代码，省得重新定义计算。
这个函数定义有点特别，设定好 R、r、L 后返回的是 Function，一个代表曲面的向量值函数。返回的函数有两个参数 u 和 v：u 从 0 到 2\[Pi] 的话，就相当于绕大圆一圈，而 v 从 0 到 2\[Pi]，则相当于绕小圆一圈。u、v 彼此独立，则互相交织形成了环面。可以用 ParametricPlot3D 把它画出来，就是一个甜甜圈的样子：
02
变化环面
我们得到了环面的参数方程，并把它画了出来，从视觉方面验证了这个方程的正确性。下面我们就可以在之前计算的基础上，变换花样了。
改变环绕曲线
前文说了，环面是一个圆 A 绕另一个圆 B 形成的曲面，但 A 和 B 只要拓扑上还是个圆，也就是自身不相交的闭合曲线，那么这么绕圈操作之后，拓扑上得到的就还是一个环面。这里我们考虑一种特殊的曲线：内摆线。内摆线定义如下，它将是我们这一节讨论的，不同于圆的环绕曲线。
在这一小节，我们只讨论改变 A 的情形，下一节“环面上的曲线”会讨论推广的改变 B 的情形。现在我们就定义一个可以任意指定环绕曲线 A 的环面函数：
这个定义和 torusSurface 十分相似，无非是把 r Cos[v] 换成 c[[1]]，把 r Sin[v] 换成 c[[2]] 而已。接下来我们把横截面曲线和对应生成的环面对比的绘制出来，环面只绘制出了四分之三，空出一角并加了厚度，让横截面更加清晰可辨。
从定义和上面几幅图可以看到，闭合曲线 A 在环绕过程中并不发生变化。注意内摆线的形状都是某种正 k 边形，我们可以考虑让 A 在环绕过程中自身绕中心旋转，只要旋转速率适当，就可以在环绕一周后，仍然形成闭合曲面。经过一番思考，可以把这个想法写成如下函数：
这个函数里的 k 表示了截面是正 k 边形，n 则表示自身旋转扭曲的程度，相比原来多旋转了 n/k 圈。下图展示了 k 分别为 3、4、5 时的环面：
再来一张固定 k 为 3，n 分别为 1、2、3 时的环面，从顶部看：
我们构造的当 k = 3、n = 1 时的扭曲环面，和雕塑家 Helaman Ferguson 的著名作品“Umbilic Torus SC”的整体形状相当一致：
改变表面起伏
对于环面上任意一点，都存在一个平面刚好和环面相切，我们称之为环面在这一点的切平面。垂直于切平面的单位向量称之为环面在这一点的法向量，垂直于平面的法向量有两个，一个朝上一个朝下。环面是一个闭曲面，把空间划分成了内外两部分，所以环面上一点的两个法向量可以说一个朝里一个朝外。让环面沿着其法向量的方向有规律的起伏，可以形成一些有趣的形状。
从法向量的定义可以知道，要计算法向量，关键在于计算切平面，而环面上一点切平面则取决于该点上两个不共线的切向量。我们的环面参数方程可以看作二元向量函数，它的两个一阶偏导数刚好就可以用作切向量，它们叉积的方向就是该点上一个法向量的方向，改变叉积顺序就能得到相反方向。我们可以适当选择叉积顺序来得到朝外的法向量，以方便后续计算。
利用 Mathematica 的符号好计算功能，我们很容易得到环面上法向量的公式，计算很直接，求偏导数，求叉积，然后在归一化，因为法向量是单位向量。
据此我们可以定义环面上的法向量函数如下：
然后我们就可以定义不同的 “起伏函数” 来改变环面的形状了。起伏函数是一个标量函数，表示沿法向量起伏的高度就可以了。新形状参数式总是这样的：环面 + 起伏 * 法向量
下面我们举几个例子，首先可以让环面沿被环绕的大圆的方向涨缩，定义起伏函数如下，m 表示起伏高度，num 表示涨缩次数：
我们也可以让环面长出一些尖刺：
更有意思的是，我们可以给上述计算再添加一个参数 t，根据不同的 t 生成不同的图像，然后输出成动画。特别是我们可以精心挑选 t 的范围，让这个动画能够首尾平滑连接，形成无限循环的 GIF 动画效果。下面就是在尖刺圆环的起伏函数的 v 方向，即环绕的小圆方向上加上参数 t，形成一种不断吞噬的效果：
然后输出成 GIF 动图：
上图的尖刺分布有点太过“平实”了，我们还可以定义一个新的起伏函数，让尖刺能够交错分布：
然后一样我们可以生成动图：
上面定义的两个起伏函数都不小于零，可以说有起无伏，为了让起伏函数名副其实，我们再定义一个，看看效果：
环面上的曲线
横向弹簧
之前的计算里，我们都用了 u、v 两个互相独立的参数来生成曲面。如果它们不是互相独立的会怎么样呢？单个参数会生成曲线，而由于我们用环面函数做基础，那么很大可能我们得到的是环面上的曲线。这就来试试看好了：
这看起来就像是一圈弹簧，为什么是这个形状呢？这是可以解释的：前面环面函数 torusSurface[R, r][u, v] 中，u、v 彼此独立相互交织，形成环面。现在不是独立的了，u 绕大圆的同时，v 以 9 倍于 u 的速度绕小圆，这不就形成螺旋了嘛。上面把 v 固定成 9u 得到的是很均匀的"弹簧"形状，未免有点无趣，我们可以把 v 变成更加复杂的关于 u 的函数，再加上和之前类似的添加额外参数 t 的技术，就可以得到比较有意思的动画了：
动画虽然有趣，但光这样的线条还是有些单调了，我们可以把它变成管状。这就是我们之前提到的用圆 A 绕圆 B，但现在 B 不再是一个圆的情形。所用的方法和一开始我们从圆得到环面是一样的：找到沿着曲线的一个局部标架，然后就可以得到参数方程了。幸运的是，微分几何中提供了空间曲线一种局部标架：Frenet 标架。要计算这种标架，先要知道的是曲线的一阶和二阶导数，这可以说是 Mathematica 的强项：
有了一阶和二阶导数就可以得到组成 Frenet 标架的三个单位向量，进而定义上述"弹簧"管状版本的函数：
我们用和之前相似的做法，可以到的上述管状曲线的 GIF 动画：
纵向弹簧
还可以交换一下 u 和 9u 的位置，曲线变成小圆旋转一圈时，大圆旋转了九圈，这样就是一个纵向的弹簧了，而且由于环面的关系，这个弹簧是内外两层连在一起形成的。为了美观，我们拉伸了环面 Z 方向的高度。
动画效果也还不错：
然后我们仍然可以得到一个管状版本，并输出动画：
上述的计算过程非常繁琐，运行速度比较慢，我用这种方法只是为了展示"管状化"的数学原理。Mathematica 里提供了一个 Tube 函数，可以把折线、曲线"管状化"，我们用它的话速度会非常快，直接计算曲线上的 1000 个点，然后用 Tube 转换成管状：
环面扭结
前面的两个例子展现了曲线 torusSurface[R, r][u, 9u] 和 torusSurface[R, r][9u, u] 两种情况，如果最后的参数是一般性的 [n*u, m*u] 会是什么样子呢？当 n 和 m 互素时，我们得到的是一大类纽结，叫做"环面纽结"。下面就用 Tube 展示了几种环面纽结：
从顶部看这些环面纽结是这个样子的：
进一步推广
我们最后再举一个例子。著名的荷兰画家 Escher 有一幅版画 Spirals，可以看作四条螺旋带组成的，不断缩小小圆半径的环面：
我们只要让环面函数 torusSurface[R, r][u, v] 里的 r 随着 u 不断变化，同时限制 v 的变化范围，就能生成一条螺旋带，然后把四条这样的螺旋带放在一起，就可以构造出和 Spirals 里极为相似的图形了。这里我们再次用到了 ParametricPlot3D 里的 PlotTheme->"ThickSurface" 选项，让曲面带有厚度：
在环面上能做的变化还有很多很多，限于篇幅这里就不继续列举下去了。希望目前为止的这些展示能让读者有所启发，有兴趣的读者可以进一步继续探索其它可能的变化。
END