这是均匀三角网格体引力计算相关的第二部分~
在第一部分中我们已经计算出了这个引力的精确解析形式，这当然很好，但每计算一次引力就需要对每个三角面做一遍求和，这有时仍然代价过高。如果能计算出引力场的球谐系数，那么在大多数情况下，使用有限的几阶球谐场进行近似就已经足够了。这也是第一部分中提出的两个问题之一，这篇文章将彻底解决这个问题。
此外，在这篇文章的最后两章还给出了均匀三轴椭球体的各阶球谐系数的解析算式，以及如何估计非均匀天体的球谐系数。
第五章 预备知识
进一步的计算需要一些球谐函数的知识，这里不加证明地直接列出。定义球谐函数：
\[Y_n^m(\lambda,\phi)=(-1)^m\sqrt{\frac{2n+1}{4\pi}\frac{(n-m)!}{(n+m)!}}P_n^m(\sin\phi)\mathrm{e}^{\mathrm{i}m\lambda}\]
注意本文采取球坐标的经度\(\lambda\)与纬度\(\phi\)作为球谐函数的自变量，这与数学中常用的极角\(\theta\)和经度\(\phi\)不同；\(P_n^m\)是连带勒让德函数：
\[P_n^m(x)=\frac{1}{2^nn!}(1-x^2)^{m/2}\frac{\mathrm d^{n+m}}{\mathrm d x^{n+m}}(x^2-1)^n\]
注意本文采用的约定并未将球谐函数的Condon–Shortley相位\((-1)^m\)放入连带勒让德函数的定义。这与Mathematica中不同，但与c++标准一致。
当\(m=0\)时，连带勒让德函数\(P_n^m\)成为勒让德多项式：
\[P_n(x)=P_n^0(x)\]
球谐函数有正交归一关系：
\[\int Y_n^m(\lambda,\phi)Y_{n’}^{m’}(\lambda,\phi)^*\mathrm{d}\Omega=\delta_{nn’}\delta_{mm’}\]
其中\(\mathrm{d}\Omega=\cos\phi\mathrm{d}\phi\mathrm{d}\lambda\)为立体角微元，积分取遍\(4\pi\)立体角；\(^*\)为取复共轭；\(\delta\)为Kronecker符号。
\(Y_n^{-m}\)可用\(Y_n^m\)表示：
\[Y_n^{-m}(\lambda,\phi)=(-1)^mY_n^m(\lambda,\phi)^*\]
最后，当两矢量\(\vec r\)、\(\vec r’\)满足\(|\vec r’|\lt|\vec r|\)时，如下恒等式成立：
\[\frac{1}{|\vec r-\vec r’|}=4\pi\sum_{n=0}^\infty\frac{1}{2n+1}\frac{r’^{n}}{r^{n+1}}\sum_{m=-n}^n Y_n^m(\lambda,\phi)Y_n^m(\lambda’,\phi’)^*\]
其中\((r,\lambda,\phi)\)和\((r’,\lambda’,\phi’)\)分别为\(\vec r\)和\(\vec r’\)的模长、经度和纬度。这些内容的详情可见这篇教案中的式(7~11)和(32)。
第六章 任意天体引力势的球谐展开
这一章先不考虑均匀三角网格体……而是先考虑一个任意天体，给出尽可能通用的公式。
首先按定义写出天体的引力势场去除质点项的部分：
\[\begin{aligned}\Delta\phi(\vec r)&=\left(-G\int_{\vec r’\in V}\frac{\rho(\vec r’)\mathrm{d}V}{|\vec r-\vec r’|}\right)-\left(-\frac{GM}{r}\right)\\
&=-G\int_{\vec r’\in V}\left(\frac{1}{|\vec r-\vec r’|}-\frac{1}{r}\right)\rho(\vec r’)\mathrm{d}V
\end{aligned}\]
其中\(V\)为天体体积；\(M=\int\rho\mathrm{d}V\)为天体质量；\(\vec r’\)为天体上质量微元的位置；\(\rho(\vec r’)\)为天体密度；\(r=|\vec r|\)为场点\(\vec r\)与原点的距离；原点设在天体的质心处。
利用上一章中提到的恒等式，将\(1/{|\vec r-\vec r’|}\)展开：
\[\begin{aligned}
\Delta\phi(\vec r)=&-G\int_{\vec r’}\rho(\vec r’)\mathrm{d}V\left(-\frac{1}{r}+4\pi\sum_{n=0}^\infty\frac{1}{2n+1}\frac{r’^{n}}{r^{n+1}}\sum_{m=-n}^n Y_n^m(\lambda,\phi)Y_n^m(\lambda’,\phi’)^*\right)\\
=&-G\int_{\vec r’}\rho(\vec r’)\mathrm{d}V\left(-\frac{1}{r}+4\pi\frac{1}{r}Y_0^0(\lambda,\phi)Y_0^0(\lambda’,\phi’)^*\right.\\
&\quad+\left.4\pi\sum_{n=1}^\infty\frac{1}{2n+1}\frac{r’^{n}}{r^{n+1}}\sum_{m=-n}^n Y_n^m(\lambda,\phi)Y_n^m(\lambda’,\phi’)^*\right)\\
=&-G\int_{\vec r’}\rho(\vec r’)\mathrm{d}V\left(-\frac{1}{r}+4\pi\frac{1}{r}\sqrt{\frac{1}{4\pi}}\sqrt{\frac{1}{4\pi}}\right.\\
&\quad+\left.4\pi\sum_{n=1}^\infty\frac{1}{2n+1}\frac{r’^{n}}{r^{n+1}}\right(Y_n^0(\lambda,\phi)Y_n^0(\lambda’,\phi’)^*\\
&\quad\quad+\left.\left.\sum_{m=1}^n (Y_n^{-m}(\lambda,\phi)Y_n^{-m}(\lambda’,\phi’)^*+Y_n^m(\lambda,\phi)Y_n^m(\lambda’,\phi’)^*)\right)\right)\\
=&-G\left.\int_{\vec r’}\rho(\vec r’)\mathrm{d}V\right(\quad 0\\
&\quad+\left.4\pi\sum_{n=1}^\infty\frac{1}{2n+1}\frac{r’^{n}}{r^{n+1}}\right(\sqrt{\frac{2n+1}{4\pi}}P_n(\sin\phi)\sqrt{\frac{2n+1}{4\pi}}P_n(\sin\phi’)\\
&\quad\quad+\left.\left.\sum_{m=1}^n (Y_n^m(\lambda,\phi)^*Y_n^m(\lambda’,\phi’)+Y_n^m(\lambda,\phi)Y_n^m(\lambda’,\phi’)^*)\right)\right)\\
=&-G\left.\int_{\vec r’}\rho(\vec r’)\mathrm{d}V\sum_{n=1}^\infty\frac{r’^{n}}{r^{n+1}}\right(P_n(\sin\phi)P_n(\sin\phi’)\\
&\quad+\left.\frac{4\pi}{2n+1}\sum_{m=1}^n 2\mathrm{Re}(Y_n^m(\lambda,\phi)Y_n^m(\lambda’,\phi’)^*)\right)
\end{aligned}\]
为了避免一屏写不下，先歇会~~
以上推导过程先将\(n=0\)的项抽出来，由于\(Y_0^0=1/\sqrt{4\pi}\)，此项刚好和\(1/r\)抵消。这代表零阶球谐势场等价于天体的质点势场；再将对\(m\)的求和拆分，利用\(Y_n^{-m}(\lambda,\phi)=(-1)^mY_n^m(\lambda,\phi)^*\)将对\(\pm m\)的求和配对相加，最后利用\(z_1^*z_2+z_1z_2^*=2\mathrm{Re}(z_1z_2^*)\)继续简化。
继续把球谐函数用连带勒让德函数写出来，其实已经没几步了……
\[\begin{aligned}
\Delta\phi(\vec r)=&-G\left.\int_{\vec r’}\rho(\vec r’)\mathrm{d}V\sum_{n=1}^\infty\frac{r’^{n}}{r^{n+1}}\right(P_n(\sin\phi)P_n(\sin\phi’)\\
&+\left.\sum_{m=1}^n 2\frac{(n-m)!}{(n+m)!}P_n^m(\sin\phi)P_n^m(\sin\phi’)\mathrm{Re}\left(\mathrm{e}^{\mathrm i m(\lambda-\lambda’)}\right)\right)\\
=&-G\left.\int_{\vec r’}\rho(\vec r’)\mathrm{d}V\sum_{n=1}^\infty\frac{r’^{n}}{r^{n+1}}\right(P_n(\sin\phi)P_n(\sin\phi’)\\
&+\left.\sum_{m=1}^n 2\frac{(n-m)!}{(n+m)!}P_n^m(\sin\phi)P_n^m(\sin\phi’)(\cos(m\lambda)\cos(m\lambda’)+\sin(m\lambda)\sin(m\lambda’))\right)
\end{aligned}\]
算到这里就是成功！现在翻出行星形变理论（3）里以球谐系数\(J_n, C_{nm}, S_{nm}\)表达的\(\Delta\phi\)的式子：
\[\begin{aligned}\Delta\phi(\vec r)=&\frac{GM}{r}\left(\sum_{n=2}^{\infty} J_n\left(\frac{R}{r}\right)^nP_n(\sin\phi)\right.\\
&\left.-\sum_{n=2}^{\infty}\left(\frac{R}{r}\right)^n \sum_{m=1}^{n}P_n^m(\sin\phi)(C_{nm}\cos m\lambda+S_{nm}\sin m\lambda)\right)
\end{aligned}\]
两式分别对照\(P_n(\sin\phi)\)、\(P_n^m(\sin\phi)\cos(m\lambda)\)和\(P_n^m(\sin\phi)\sin(m\lambda)\)的系数，不难看出：
\[J_n=-\frac{1}{MR^n}\int_{\vec r’}r’^{n}P_n(\sin\phi’)\rho(\vec r’)\mathrm{d}V\\
\begin{pmatrix}C_{nm}\\S_{nm}\end{pmatrix}=2\frac{(n-m)!}{(n+m)!}\frac{1}{MR^n}\int_{\vec r’}r’^{n}P_n^m(\sin\phi’)\begin{pmatrix}\cos(m\lambda’)\\\sin(m\lambda’)\end{pmatrix}\rho(\vec r’)\mathrm{d}V\\
\]
终于~~我们有了用体积分表示的任意天体球谐系数的计算式了~~
不过还有一点遗留问题：球谐系数的下标\(n\)是从\(2\)开始的，那么\(n=1\)的项哪去了？这是因为计算球谐系数的积分式中被积函数（除质量微元\(\rho(\vec r’)\mathrm{d}V\)外的部分）在完全展开后是坐标\((x’,y’,z’)=\vec r’\)的\(n\)次齐次多项式。当\(n=1\)时，被积函数就是线性函数，对质量微元的体积分可以仅用总质量和质心位置表达。实际上，天体的质心\(\vec c\)可以使用一阶球谐系数表示为\(\vec c=(C_{11}R, S_{11}R, -J_{1}R)\)，此式的证明留作习题。通过选取坐标系原点与天体质心重合，可以消除掉这些系数。这样做只有好处，而且符合直觉，一般都约定这样做，所以从不需要讨论或计算一阶球谐系数。
第七章 均匀三角网格体的球谐系数
呃……
剩下的内容好像没有什么可以说的了……由于已经有了能算球谐系数的积分式，那只需要对三角网格体进行这些积分就行了。积分所用的方法和第一部分一样，都是先找到一个矢量函数，其散度等于被积函数，然后利用高斯定理将体积分转换为面积分，再利用三角面参数化换元对每个三角面解析计算这个面积分。
每个球谐系数\(h\)的计算可以写成如下形式：
\[h=\int_{\vec r’} p(\vec r’)\mathrm{d}V\]
被积函数\(p(\vec r’)\)都是坐标的\(n\)次齐次多项式，故有多个矢量场\(\vec q^\alpha\)满足\(\nabla’\cdot\vec q^\alpha(\vec r’)=p(\vec r’)\)可供选择:
\[\vec q^x=\hat x’\int p(\vec r’)\mathrm{d}x’\\
\vec q^y=\hat y’\int p(\vec r’)\mathrm{d}y’\\
\vec q^z=\hat z’\int p(\vec r’)\mathrm{d}z’
\]
对\(\alpha\in{\{x,y,z\}}\)，\(\vec q^\alpha\)均为坐标的\(n+1\)次齐次多项式。之后可利用高斯定理：
\[\begin{aligned}
h&=\int_{\vec r’} p(\vec r’)\mathrm{d}V\\
&=\int_{\vec r’} \nabla’\cdot\vec q^\alpha(\vec r’)\mathrm{d}V\\
&=\int_{\partial V} \vec q^\alpha(\vec r’)\cdot\mathrm{d}\vec S’\\
&=\sum_{i=1}^N\int_{T_i} \vec q^\alpha(\vec r’)\cdot\hat n\mathrm{d}S’\\
&=\sum_{i=1}^N L_h^\alpha(T_i)
\end{aligned}
\]
其中，
\[L_h^\alpha(T)=\int_{T}\vec q^\alpha(\vec r’)\cdot\hat n\mathrm{d}S’\]
可以继续使用如第一部分的换元法计算。由于被积函数均为多项式，结果也将为三角形各顶点坐标位置\(T=\Delta\vec p_1\vec p_2\vec p_3\)的多项式。实际计算中，可选取使\(L_h^\alpha(T)\)表达式最简单的一个\(\alpha\)。
这些计算没有技术上的难度，因为都是多项式积分，但十分繁琐。从这里可以下载一份Mathematica代码。代码中实现了第一、第二部分中描述的所有公式与算法：首先将球谐系数计算式中的被积函数展开为多项式，再对这些多项式积分并编译结果以供快速计算三角网格体的体积、质心及球谐系数，最后在一个测试三角网格模型上比较球谐场和第一部分中描述的精确解析算法的差异。
最后，本文描述的算法已由c++实现并整合进入开源仓库GrayscaleGenerator，下载即可计算任意封闭均匀三角网格的球谐势场。不过这个仓库的Geopotential输出是供KSP/Principia使用的，而Principia使用的归一化约定与我们不同，故其输出的各个系数\(\tilde C_{nm}, \tilde S_{nm}\)都与本文中的定义相差了一个因子。具体而言，它输出的是
\[\begin{aligned}
m=0:\quad &\tilde C_{n0}=-\frac{J_n}{\sqrt{2n+1}},\quad \tilde S_{n0}=0\\
m\gt 0:\quad &\begin{pmatrix}\tilde C_{nm}\\\tilde S_{nm}\end{pmatrix}=\begin{pmatrix}C_{nm}\\S_{nm}\end{pmatrix}\sqrt{\frac{1}{2(2n+1)}\frac{(n+m)!}{(n-m)!}}
\end{aligned}
\]
第八章 均匀三轴椭球体的球谐系数
考虑一个均匀的三轴椭球体，三半长轴分别为\(a,b,c\)，且分别与\(xyz\)坐标轴对齐，坐标原点为椭球体中心：
\[V=\left\{(x,y,z)\left|\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}\right.\le 1\right\}\]
这样的椭球体的球谐系数是可以直接使用第六章中的体积分公式算出的：因为被积函数是坐标的齐次多项式，而被积区域可以进行经过缩放的球坐标换元。换元后被积函数成为各项可直接分离变量的三角函数，没有任何难度，但同样较为繁琐。这里直接给出最初几阶非\(0\)的球谐系数的计算式；所有的\(S_{nm}\)或下标不全为偶数的\(J_n,C_{nm}\)均为\(0\)：
\[\begin{aligned}
J_2&=\frac{2}{5}t& C_{22}&=\frac{2}{5}s&\\
J_4&=-\frac{12}{35}(2s^2+t^2)& C_{42}&=-\frac{4}{35}st\\ & & C_{44}&=\frac{1}{35}s^2\\
J_6&=\frac{8}{21}(6s^2t+t^3)& C_{62}&=\frac{8}{105}(s^3+st^2)\\ & & C_{64}&=-\frac{2}{315}s^2t\\ & & C_{66}&=\frac{1}{945}s^3\\
J_8&=-\frac{16}{33}(6s^4+12s^2t^2+t^4)& C_{82}&=-\frac{16}{231}(3s^3t+st^3)\\ & & C_{84}&=\frac{4}{3465}(2s^4+3s^2t^2)\\ & & C_{86}&=-\frac{2}{10395}s^3t\\ & & C_{88}&=\frac{1}{41580}s^4
\end{aligned}\]
其中，\(s,t\)为与\(a,b,c\)有关的参数：
\[s=\frac{a^2-b^2}{8R^2}\quad\quad t=\frac{a^2+b^2-2c^2}{4R^2}\]
\(R\)为参考半径，一般可取为体积平均半径\(R=(abc)^{1/3}\)。其具体取值会影响球谐系数的数值，但不会影响引力场的计算。
经观察，我们注意到，这些系数的计算式可以写为统一的形式：
\[
J_n=-\frac{1}{2}Q(n,0)\quad\quad C_{nm}=Q(n,m)\quad\quad S_{nm}=0\\
Q(n,m)=\left\{\begin{aligned}
&\frac{6(-1)^{\frac{n+m}{2}}(n-m)!}{(n+3)!!}\sum_{j=0}^{\lfloor\frac{n-m}{2}\rfloor}\frac{s^{2j+\frac{m}{2}}t^{\frac{n-m}{2}-2j}}{j!\left(j+\frac{m}{2}\right)!\left(\frac{n-m}{2}-2j\right)!},& &(n,m)~\text{are even,}\\
&0,& &\text{else.}
\end{aligned}\right.
\]
感兴趣的读者可自行验证，有能力的读者可尝试证明~~~
第⑨章 非均匀情形的球谐系数
啊……
虽然第七章中给出的方法以及代码已经可以计算均匀三角网格体的球谐系数了，但其实用性还是有些受限：实际天体不可能是密度均匀的。一般来说天体中心密度会更高，而表面密度较低。所以本章就讨论一个简单的一维密度分布模型：\[\rho(\vec r)=f\left(\frac{|\vec r|}{|\vec r_S|}\right)=f(\xi)\]
其中\(\vec r_S\)为从行星中心（即坐标原点）出发，过\(\vec r\)的射线与行星表面的交点；无量纲参数\(\xi=|\vec r|/|\vec r_S|\in [0,1]\)描述\(\vec r\)位于行星内部的比例。这个模型假定了行星密度可以由一元函数\(f(\xi)\)描述，在大多数情况下这是一个较好的近似。
首先考虑在三轴椭球体情形下加入这么一个密度分布会造成什么影响。已知\(n\)阶球谐系数\(h_n\)的计算式是对一个关于坐标的\(n\)次多项式\(p_n(\vec r)\)的体积分。\(p_n\)中还包含一些与\(M\)、\(\rho(\vec r)\)有关的项，我们可以将这些项显式地提出来，把剩下的仅与坐标相关的多项式记为\(\tilde p_n\)：
\[\begin{aligned}
h_n&=\int_{\vec r’} p_n(\vec r’)dV\\
&=\frac{1}{MR^n}\int_{\vec r’} \tilde p_n(\vec r’)\rho(\vec r’)dV\\
&=\frac{1}{R^n\int_{\vec r’} \rho(\vec r’)dV}\int_{\vec r’} \tilde p_n(\vec r’)\rho(\vec r’)dV
\end{aligned}\]
现在进行椭球坐标换元：
\[\vec r’=\begin{pmatrix}x’\\y’\\z’\end{pmatrix}=\begin{pmatrix}a\xi\cos\phi\cos\lambda\\b\xi\cos\phi\sin\lambda\\c\xi\sin\phi\end{pmatrix},\quad\quad
dV=abc\xi^2\mathrm{d}\xi\mathrm{d}\Omega=R^3\xi^2\mathrm{d}\xi\mathrm{d}\Omega\]
由于换元前\(\tilde p_n\)是关于\(\vec r’=(x’,y’,z’)\)的\(n\)次齐次多项式，故换元后它的每一项中都有一个公共的\(\xi^n\)，即可以拆写为\(\tilde p_n(\vec r’)=\xi^n g(\lambda,\phi)\)，其中\(g\)为仅与\(\lambda,\phi\)有关的函数（实际上是一个三角多项式）。而\(\rho(\vec r’)\)根据定义可直接由\(f(\xi)\)代替。将这些代入\(h_n\)可得：
\[\begin{aligned}
h_n&=\frac{1}{R^n\int f(\xi)R^3\xi^2\mathrm{d}\xi\mathrm{d}\Omega}\int_{\vec r’} \xi^n g(\lambda,\phi)f(\xi)R^3\xi^2\mathrm{d}\xi\mathrm{d}\Omega\\
&=\frac{1}{R^n\int f(\xi)\xi^2\mathrm{d}\xi\int\mathrm{d}\Omega}\int f(\xi)\xi^{n+2}\mathrm{d}\xi\int g(\lambda,\phi)\mathrm{d}\Omega\\
&=\frac{\int f(\xi)\xi^{n+2}\mathrm{d}\xi}{\int f(\xi)\xi^2\mathrm{d}\xi}\cdot\frac{1}{R^n\int\mathrm{d}\Omega}\int g(\lambda,\phi)\mathrm{d}\Omega\\
&=\frac{\int f(\xi)\xi^{n+2}\mathrm{d}\xi}{\int f(\xi)\xi^2\mathrm{d}\xi}\cdot\frac{\int \rho_0\xi^2\mathrm{d}\xi}{\int \rho_0\xi^{n+2}\mathrm{d}\xi}\cdot\left(\frac{1}{R^n\int \rho_0\xi^2\mathrm{d}\xi\int\mathrm{d}\Omega}\int \rho_0\xi^{n+2}\mathrm{d}\xi\int g(\lambda,\phi)\mathrm{d}\Omega\right)\\
&=\frac{\int f(\xi)\xi^{n+2}\mathrm{d}\xi}{\int f(\xi)\xi^2\mathrm{d}\xi}\cdot\frac{n+3}{3}\cdot (h_n^*)
\end{aligned}
\]
其中\(h_n^*\)为将\(h_n\)（在上式第2个等号后的形式）中的\(f(\xi)=\rho(\vec r’)\)替换为一个定常数\(\rho_0\)后的结果。
定义\(I_n\)：
\[I_n=\int_0^1 (n+1)f(\xi)\xi^n\mathrm{d}\xi=\int_0^1 f(\xi)\mathrm{d}\xi^{n+1}=\int_0^1 f\left(\eta^{\frac{1}{n+1}}\right)\mathrm{d}\eta\]
则上式可写为
\[h_n=\frac{I_{n+2}}{I_2}h_n^*\]
这意味着，考虑行星内部具有从中心到表面的密度分布函数\(f(\xi)\)后，行星的各阶球谐系数\(h_n\)可由假定其均匀计算出的球谐系数\(h_n^*\)简单地乘上一个与行星内部密度分布函数\(f\)和阶数\(n\)有关的系数得出。
而对非三轴椭球的天体，比如任意三角网格体，也可以做同样的推导，只是推导中的换元会更加复杂：雅可比行列式与\(\lambda,\phi\)的关系会和网格的具体形状有关。但无论如何，这一关系仍然可以适用，只需要使用一个合适的\(f(\xi)\)，就能估计出具有从内到外非均匀密度的天体的球谐系数。
完结撒花~~~~
