引言
之前的一篇文章（行星形变理论（3））提到了可以用高阶的球谐势场来近似非质点天体的引力场，并且给出了根据球谐系数仅使用浮点加乘快速计算引力场的表达式，以及二阶球谐势场和天体的转动惯量张量的关系（也即如何从天体的质量分布积分获得二阶球谐势场）。这些结论对大天体的星历积分已经足够，但还有一些问题：
1. 很多小天体并没有足够高阶的球谐势场测量，只有形状的模型（三角网格）。如果假定小天体的内部是密度均匀的，理论上可以对三角网格包围的体积进行积分获得任意阶球谐系数。但具体如何快速计算？
2. 即使获得了小天体的球谐势系数表，但在特别接近小天体表面附近（尤其是在小天体内部，不过此种情形实用价值不大）时，也不能用球谐势场去进行引力计算。在与天体质心的距离小于天体的最大半径（也即产生引力的质量微元距质心的最大距离）时，球谐势是发散的：累加的阶数越高，结果越离谱。
接下来大概有两篇文章会解决上述两个问题，这篇文章我要先彻底解决第二个，并且提出一个惊人的结论：
一个由三角面组成的任意封闭匀质网格体，对空间中任意一点的牛顿引力场，可以使用\(O(N)\)次初等函数计算精确求出。其中\(N\)为网格的三角面数。
其实基本就是这篇论文里的内容。
第一章 高斯定理
下图是“由三角面组成的任意封闭匀质网格体”的一个例子：
这是Metis（木卫十六）的一个模型，模型顶点到质心的距离从16.3km到31.1km不等，由\(N=720\)个三角形组成。每个三角形\(T\)由三个顶点\(\vec p_1=(x_1,y_1,z_1)\)，\(\vec p_2=(x_2,y_2,z_2)\)，\(\vec p_3=(x_3,y_3,z_3)\)组成，记为\(T=\Delta\vec p_1\vec p_2\vec p_3\)。这三个点以及它们的顺序确定了三角形的法向量\(\hat n\)：法向量指向模型的外侧，从外侧看去，三角形的三个点按顺序总是逆时针排列的。据此可以写出计算法向量的表达式\[\hat n=\frac{(\vec p_2-\vec p_1)\times(\vec p_3-\vec p_1)}{|(\vec p_2-\vec p_1)\times(\vec p_3-\vec p_1)|}\]
这个表达式不唯一。因为顶点之间的轮换对称，\(\Delta\vec p_1\vec p_2\vec p_3\)、\(\Delta\vec p_2\vec p_3\vec p_1\)和\(\Delta\vec p_3\vec p_1\vec p_2\)是完全等价的，故将表达式中对应项依此替换结果不变。不过，不能直接交换一对顶点，如变成\(\Delta\vec p_3\vec p_2\vec p_1\)，这样会逆转三个点的轮换顺序，使得法向量反向。按照定义，一个封闭的多面体不仅要被三角网格无缝隙、不重叠地拼合起来，而且所有三角形的法向都需要朝外，否则不能算是“封闭”的。
虽然模型是被表面的三角形确定的，但用暴力体积分方法计算引力需要对整个模型内部进行体素划分，所以一般复杂度要远高于\(O(N)\)，而且精度还受限于划分体素的大小，想得到精确的结果更是不可能。不过，说到体积分比面积分更加复杂，我们就不由得联想到高斯定理，这可以将任意场\(\vec A\)散度的体积分转化成场的面积分：
\[\int_V\nabla\cdot \vec A\mathrm{d}V=\int_{\partial V}\vec A\cdot\mathrm{d}\vec S\]
首先直接写出模型对空间中任意一点\(\vec r=(x,y,z)\)的引力，并且注意到这个引力场恰好可以写成一个形式更加简洁的张量场的散度：
\[\begin{aligned}\vec g(\vec r)&=G\rho\int_{\vec r’\in V}\frac{-(\vec r-\vec r’)}{|\vec r-\vec r’|^{3}}\mathrm{d}V\\
&=G\rho\int_{V}\frac{-(x-x’,y-y’,z-z’)}{\left((x-x’)^2+(y-y’)^2+(z-z’)^2\right)^{3/2}}\mathrm{d}x’\mathrm{d}y’\mathrm{d}z’\\
&=G\rho\int_{V}\nabla’\cdot\left(\frac{-I}{\sqrt{(x-x’)^2+(y-y’)^2+(z-z’)^2}}\right)\mathrm{d}x’\mathrm{d}y’\mathrm{d}z’\\
\end{aligned}\]
其中\(G\)是万有引力常量；\(\rho\)是行星密度，假设为一常数；\(\vec r’=(x’,y’,z’)\)是行星内被积体积元的位置；\(I\)为二阶单位张量；注意微分算符\(\nabla’\)作用在\(\vec r’\)上。
使用高斯定理后，上述体积分成为对三角网格表面的积分，并且可以继续拆分成对各个三角面积分的和：
\[\begin{aligned}\vec g(\vec r)&=G\rho\int_{S’\in\partial V}\left(\frac{-I}{\sqrt{(x-x’)^2+(y-y’)^2+(z-z’)^2}}\right)\cdot\mathrm{d}\vec S’\\
&=G\rho\int_{S’\in\partial V}\frac{-\mathrm{d}\vec S’}{|\vec r-\vec r’|}\\
&=G\rho\sum_{i=1}^N \int_{S’\in T_i}\frac{-\hat n_i\mathrm{d}S’}{|\vec r-\vec r’|}\\
&=G\rho\sum_{i=1}^N \vec L(T_i,\vec r)
\end{aligned}\]
其中\(T_i\)为第\(i\)个三角形；\(\hat n_i\)为该三角形的单位法向量；\(\vec r’=(x’,y’,z’)\)是被积面元的位置；\(\vec L\)为一个函数，对三角形\(T\)和点\(\vec r=(x,y,z)\)，给出：
\[\vec L(T,\vec r)=\int_{S’\in T}\frac{-\hat n\mathrm{d}S’}{|\vec r-\vec r’|}\]
如此，我们就将复杂的体积分转换成了简单（？）的面积分~~~
第二章 符号积分
简单的理论推导结束了，现在开始是暴力积分时间！
目标是对任意三角形\(T=\Delta\vec p_1\vec p_2\vec p_3\)和点\(\vec r=(x,y,z)\)计算出\(\vec L(T,\vec r)\)。首先进行一次换元，令
\[\vec r’=\vec p_1+u(\vec p_2-\vec p_1)+v(\vec p_3-\vec p_2)\]
也即将\(\Delta\vec p_1\vec p_2\vec p_3\)换元到\(uv\)平面上\(\{(u,v)|0\le u\le 1, 0\le v\le u\}\)的区域。在换元后，三角形的三个顶点\(\vec p_1\vec p_2\vec p_3\)分别对应\((u,v)=(0,0), (1,0), (1,1)\)。
由于三角形是线性平面，所以雅可比行列式非常好算，就是三角形在换元前后面积的比。换元前三角形\(T\)的面积可以写为\[S_T=\frac{1}{2}|(\vec p_2-\vec p_1)\times(\vec p_3-\vec p_1)|\]
注意到这个叉积也可以用来计算三角形的法向量，故可以定义非归一的法向量
\[\vec n=(\vec p_2-\vec p_1)\times(\vec p_3-\vec p_1)\]
这样三角形的法向量可以写为\(\hat n=\vec n/|\vec n|\)，面积为\(S_T=|\vec n|/2\)。
换元后在\(uv\)平面上的三角形面积为\(S_{uv}=1/2\)。综上，原积分在换元后成为
\[\begin{aligned}\vec L(T,\vec r)&=\int_{S’\in T}\frac{-\hat n\mathrm{d}S’}{|\vec r-\vec r’|}\\
&=\int_0^1\mathrm{d}u\int_0^u\mathrm{d}v\cdot\left|\frac{S_T}{S_{uv}}\right|\cdot\frac{-\hat n}{|\vec r-\vec r’|}\\
&=\int_0^1\mathrm{d}u\int_0^u\mathrm{d}v\cdot|\vec n|\cdot\frac{-\vec n/|\vec n|}{|\vec r-\vec r’|}\\
&=\vec n\int_0^1\mathrm{d}u\int_0^u\frac{-\mathrm{d}v}{|\vec r-\vec r’|}
\end{aligned}\]
雅可比行列式与单位法向量中的\(|\vec n|\)恰好抵消，剩余一个未归一的法向量\(\vec n\)。
接下来需要将被积表达式展开成关于\(uv\)的函数：
\[\begin{aligned}
|\vec r-\vec r’|&=\sqrt{(\vec r-\vec r’)^2}\\
&=\sqrt{((\vec r-\vec p_1)-u(\vec p_2-\vec p_1)-v(\vec p_3-\vec p_2))^2}\\
&=\sqrt{b+2ju+su^2+2(a+pu)v+cv^2}
\end{aligned}\]
其中：\[
a=(\vec r-\vec p_1)\cdot(\vec p_2-\vec p_3)\\
j=(\vec r-\vec p_1)\cdot(\vec p_1-\vec p_2)\\
p=(\vec p_1-\vec p_2)\cdot(\vec p_2-\vec p_3)\\
b=(\vec r-\vec p_1)^2\\
s=(\vec p_1-\vec p_2)^2\\
c=(\vec p_2-\vec p_3)^2\\
\]
你相信微积分吗？符号推导……要开始加速了！
\[\begin{aligned}\vec L(T,\vec r)&=\vec n\int_0^1\mathrm{d}u\int_0^u\frac{-\mathrm{d}v}{\sqrt{b+2ju+su^2+2(a+pu)v+cv^2}}\\
&=\vec n\int_0^1\mathrm{d}u\left.\left[-\frac{\ln{\left(a+pu+cv+\sqrt{c}\sqrt{b+2ju+su^2+2(a+pu)v+cv^2}\right)}}{\sqrt{c}}\right]\right|_0^u\\
&=\frac{\vec n}{\sqrt{c}}\int_0^1\ln{\left(\frac{a+pu+\sqrt{c}\sqrt{b+(2j+su)u}}{a+qu+\sqrt{c}\sqrt{b+(2k+tu)u}}\right)}\mathrm{d}u\\
\end{aligned}\]
其中\[
q=p+c=(\vec p_1-\vec p_3)\cdot(\vec p_2-\vec p_3)\\
k=j+a=(\vec r-\vec p_1)\cdot(\vec p_1-\vec p_3)\\
t=s+2p+c=(\vec p_1-\vec p_3)^2\\
\]
来自于\(v=u\)时不定积分中的对应项。可以发现对\(v\)的定积分结果具有一个对称性，可以看成是两个符号形式相同的\(\ln\)函数的差，只是第一式中的\(p,j,s\)在第二式中被换为了\(q,k,t\)。
利用此对称性，接下来只需计算\(\ln\left(a+pu+\sqrt{c}\sqrt{b+(2j+su)u}\right)\)的不定积分，然后进行替换和作差即可。
\[\begin{aligned}\mathscr{I}&=\int\ln\left(a+pu+\sqrt{c}\sqrt{b+(2j+su)u}\right)\mathrm{d}u\\
&=u\ln\left(a+pu+\sqrt{c}\sqrt{b+(2j+su)u}\right)-\int u\frac{p+\sqrt{c}\frac{j+su}{\sqrt{b+(2j+su)u}}}{a+pu+\sqrt{c}\sqrt{b+(2j+su)u}}\mathrm{d}u
\end{aligned}\]
换元，令
\[w=\mathrm{arcsinh}\frac{j+su}{\sqrt{bs-j^2}}\]
可以验证\(bs-j^2=((\vec r-\vec p_1)\times(\vec p_1-\vec p_2))^2\gt 0\)。记\(h=\sqrt{bs-j^2}\)，则有
\[\begin{aligned}u&=\frac{h\sinh{w}-j}{s}\\
\sqrt{b+(2j+su)u}&=\frac{h}{\sqrt{s}}\cosh{w}\end{aligned}\]
继续计算不定积分。这是一个双曲有理式，不定积分必为初等函数：
\[\begin{aligned}\mathscr{I}=~&u\ln\left(a+pu+\sqrt{c}\sqrt{b+(2j+su)u}\right)\\
&\quad-\int\frac{h(h\sinh{w}-j)(p\cosh{w}+\sqrt{c}\sqrt{s}\sinh{w})}{s(as-jp+h(p\sinh{w}+\sqrt{c}\sqrt{s}\cosh{w}))}\mathrm{d}w\\
=~&u\ln\left(a+pu+\sqrt{c}\sqrt{b+(2j+su)u}\right)-\left.\frac{1}{s(cs-p^2)}\cdot\right(\\
~&(cs-p^2)h\sinh{w}-(as-jp)w\sqrt{c}\sqrt{s}\\
-&(cj-ap)s\ln(as-jp+h(p\sinh{w}+\sqrt{c}\sqrt{s}\cosh{w}))\\
-&\left.2\sqrt{c}\sqrt{s}\sqrt{h^2(cs-p^2)-(as-jp)^2}\arctan\frac{hp-(as-jp-h\sqrt{c}\sqrt{s})\tanh(w/2)}{\sqrt{h^2(cs-p^2)-(as-jp)^2}}\right)\\
\end{aligned}
\]
不论这一堆东西如何复杂，它总是初等函数。所以到此为止，我们已经论证了开头给出的结论：均匀封闭三角网格体对空间中任一点的引力场可以用初等函数写出。但为了实际计算这个引力场，还得继续化简这一坨式子……
第三章 符号化简
上面那坨好几行的式子只是不定积分\(\mathscr{I}\)的结果，为了计算\(\vec L(T,\vec r)\)需要首先将\(w\)代换回\(u\)，然后将那些式子重复四次：积分的上限减去下限这是两次，再将\(p,j,s\)换为\(q,k,t\)后又是两次。总之详细的化简过程此处略去，感兴趣的读者可以自行尝试，否则这篇文章里的公式实在会突破天际……
化简结果如下：
\[\begin{aligned}\vec L(T,\vec r)=\frac{\vec n}{\vec n^2}\cdot\left(\frac{\vec e_1\cdot\vec e_2~\vec e_1\cdot\vec r_1-\vec e_1^2~\vec r_1\cdot\vec e_2}{|\vec e_1|}\ln\frac{\vec e_1\cdot\vec r_1+|\vec e_1||\vec r_1|}{\vec e_1\cdot\vec r_2+|\vec e_1||\vec r_2|}\right.&\\
+\frac{\vec e_2\cdot\vec e_3~\vec e_2\cdot\vec r_2-\vec e_2^2~\vec r_2\cdot\vec e_3}{|\vec e_2|}\ln\frac{\vec e_2\cdot\vec r_2+|\vec e_2||\vec r_2|}{\vec e_2\cdot\vec r_3+|\vec e_2||\vec r_3|}&\\
+\frac{\vec e_3\cdot\vec e_1~\vec e_3\cdot\vec r_3-\vec e_3^2~\vec r_3\cdot\vec e_1}{|\vec e_3|}\ln\frac{\vec e_3\cdot\vec r_3+|\vec e_3||\vec r_3|}{\vec e_3\cdot\vec r_1+|\vec e_3||\vec r_1|}&\\
+2\vec r_1\cdot(\vec r_2\times\vec r_3)\left.\arctan\frac{\vec r_1\cdot(\vec r_2\times\vec r_3)}{|\vec r_1||\vec r_2||\vec r_3|+|\vec r_1|~\vec r_2\cdot \vec r_3+|\vec r_2|~\vec r_3\cdot \vec r_1+|\vec r_3|~\vec r_1\cdot \vec r_2}\right)&
\end{aligned}\]
其中\[\vec e_1=\vec p_2-\vec p_1,\quad\vec e_2=\vec p_3-\vec p_2,\quad\vec e_3=\vec p_1-\vec p_3\]是沿三角形的三条边的有向矢量；\[\vec r_1=\vec p_1-\vec r,\quad\vec r_2=\vec p_2-\vec r,\quad\vec r_3=\vec p_3-\vec r\]是从场点\(\vec r\)指向三角形的三个顶点的矢量；\[\vec n=\vec e_1\times\vec e_2=\vec e_2\times\vec e_3=\vec e_3\times\vec e_1\]是三角形\(T\)的未归一化的法向量。
注意此表达式中的\(\arctan\)为二元 \(\arctan\) 函数，即\(\arctan\frac{y}{x}\equiv\arg(x+\mathrm{i}y)\)。如果直接将自变量的分数视为除法然后计算一元\(\arctan\)，在极少数的情况下会得出错误结果。
之后只需对每个三角面\(T_i\)计算\(\vec L(T_i,\vec r)\)再求和即可计算出网格体的引力场\(g(\vec r)=G\rho\sum_{i=1}^N \vec L(T_i,\vec r)\)。也即对\(N\)个三角面围成的均匀网格体，只需\(3N\)次\(\ln\)计算及\(N\)次\(\arctan\)计算即可精确求得其对空间中任意一点的引力场。还可以注意到式中的3个\(\ln\)的自变量实际上只各与三角形的一条边有关，而这条边是与两个三角形共用的，所以可以对两个三角形的公用边只计算一次\(\ln\)，降低一半\(\ln\)的计算次数。
最后在写程序数值计算\(\vec L\)时还需要注意表达式的奇点：如果式中任意一个分母或者\(\ln\)的自变量趋于\(0\)，那么实际上其所在的那一项也是趋于\(0\)的（除了\(\arctan\)项的自变量外，因为二元\(\arctan\)的自变量分母为\(0\)属于正常情况，并非奇点）。这一结论不仅适用于式中的三个\(\ln\)，对整个表达式也适用：如果\(\vec n^2=0\)，那么\(\vec L\)也是\(\vec 0\)。此外还需要注意在\(\ln\)的自变量的分式中，分子分母其实都具有\(|\vec e_i||\vec r_j|(1+\cos\theta_{ij})\)的形式，其中\(\theta_{ij}\)是\(\vec e_i\)和\(\vec r_j\)的夹角。显然它们理论上都是非负的，但数值计算时，浮点误差可能导致极接近\(0\)的值被误计算为负数，使\(\ln\)爆炸。所以在计算时如果发现分子分母之一非正，直接将此项\(\ln\)赋\(0\)即可。加入这些特判后，程序可以对一切退化情形（如场点\(\vec r\)在网格体顶点、棱所在直线、表面所在平面上，或者网格体包含面积为\(0\)的三角形）正确计算引力场。
（未完待续？）
