这是行星形变理论的第三部分。主要描述一下形变理论在实际星历积分中遇到的问题，以及实现上的解决方案。关于理论部分的描述，可见第一部分和第二部分。
啊……
之前讨论的形变理论，都假定行星初始是完美的球形，一切非质点引力项都来自于自转和潮汐引起的形变。但实际天体中，行星在不受影响时仍可能具有非球对称的形状（尤其是小行星）。且对于某些自转速度过快的天体（如妊神星），其平衡形态为三轴椭球体，也是本理论的自转形变矩阵的算法无法描述的。本篇将探讨并解决这些问题。
A1.以球谐系数表示的行星重力场非质点引力项
对一颗实际天体，描述其非质点引力项的物理量为引力势场的各阶球谐系数，即重力场模型。这些系数一般在行星的表面坐标系\(x_sy_sz_s\)中给出，其中\(x_s\)轴由行星中心指向本初子午线与赤道的交点，\(z_s\)为行星极轴，\(y_s\)轴由右手坐标系确定。注意，对同一个重力场，这些系数的归一化约定有很多种，不同约定间的同一个系数可差一常数倍。我们采用NASA/JPL在DE星历积分时使用的约定。在此约定下，行星因自身非质点引力项与外界质点\(m^*\)作用产生的额外加速度\(\vec a^*\)可在\(\xi\eta\zeta\)坐标系中给出。其中\(\xi\)轴由行星中心指向质点\(m^*\)，\(\xi\zeta\)平面包含行星极轴，\(\eta\)轴由右手坐标系确定。在这一坐标系中，有：
\[
\begin{aligned}
\begin{pmatrix}
\vec a^*\cdot\hat\xi\\
\vec a^*\cdot\hat\eta\\
\vec a^*\cdot\hat\zeta
\end{pmatrix}=&-\frac{Gm^*}{r^2}\left(\sum_{n=2}^{N_z} J_n\left(\frac{R}{r}\right)^n\begin{pmatrix}
(n+1)P_n(\sin\phi)\\
0\\
-\cos\phi P_n'(\sin\phi)
\end{pmatrix}\right.\\
&\left.+\sum_{n=2}^{N_t}\left(\frac{R}{r}\right)^n \sum_{m=1}^{n}\begin{pmatrix}
-(n+1)P_n^m(\sin\phi)(+C_{nm}\cos m\lambda+S_{nm}\sin m\lambda)\\
m\sec\phi P_n^m(\sin\phi)(-C_{nm}\sin m\lambda+S_{nm}\cos m\lambda)\\
\cos\phi {P_n^m}'(\sin\phi)(+C_{nm}\cos m\lambda+S_{nm}\sin m\lambda)
\end{pmatrix}\right)
\end{aligned}
\]
其中\(m^*\)是质点的质量；\(r\)是行星与\(m^*\)的质心距离；\(R\)是行星半径；\(J_n\)、\(C_{nm}\)与\(S_{nm}\)是球谐展开的系数，也即重力场模型系数；\(N_z\)和\(N_t\)是球谐展开的阶数上限；\(P_n\)是勒让德多项式；\(P_n^m\)是连带勒让德函数，注意此处的连带勒让德函数的定义与Mathematica中不同，不包含Condon–Shortley相位\((-1)^m\)；\(\lambda\)与\(\phi\)分别是质点\(m^*\)在行星表面坐标系\(x_sy_sz_s\)中的经度与纬度；撇号代表求导。
上式的直接数值计算极为繁杂，且在质点接近行星极轴时会产生\(0\times\infty\)型不定式。实际计算中，对固定的阶数\(n\)，可以对上式预先进行符号计算展开，结果在行星表面坐标系\(x_sy_sz_s\)中是质点坐标的齐次多项式函数。球谐展开的最低阶数是\(n=2\)，因为一阶系数可以通过将坐标原点设为行星质心而消除。更高阶的球谐势场随距离衰减更快，且一般对最终结果贡献更小。对行星星历的精度来说，计算到\(n=6\)已经绰绰有余。
计算得到行星自身加速度\(\vec a^*\)后，由反作用力即可得质点\(m^*\)受到的重力加速度，也即行星的非质点引力场。对该引力场进行一次积分，即可得到行星在自身外一点\(\vec r\)处的非质点引力势：\[\Delta\phi(\vec r)=-\int_\infty^{\vec r} \left(-\frac{M}{m^*}\vec a^*(\vec r’)\right)\cdot\mathrm d \vec r’\]
其中\(M\)是行星的质量。对\(n=2\)的情形，如上计算（过程留做习题）后可得\[
\Delta\phi_2(\vec r)=-\frac{GMR^2}{r^5}\vec r\cdot H_2\vec r\]
其中\(H_2\)为一张量。在行星表面坐标系\(x_sy_sz_s\)（即球谐系数所在的坐标系）中，\(H_2\)的分量如下：
\[
H_2=\frac{1}{2}\begin{pmatrix}
J_2+6C_{22} & 6S_{22} & 3C_{21}\\
6S_{22} & J_2-6C_{22} & 3S_{21}\\
3C_{21} & 3S_{21} & -2J_2
\end{pmatrix}
\]
注意到\(\Delta\phi_2\)的形式与我们之前所得的形变附加势非常相似，而且\(H_2\)也是零迹对称矩阵。比较后可发现形变理论的结论可写为\(H_2=(k-1)C\)，其中\(k\)与形变矩阵\(C\)的定义可见前两章。
从这里可以下载一份利用上述一堆勒让德和三角函数计算加速度\(\vec a_m\)、将结果展开成为关于坐标的多项式、再积分得到势场\(\Delta\phi\)的Mathematica代码。这些计算非常耗时，所以随代码还附带了直到32阶的展开结果，可直接加载使用。使用展开的多项式计算非质点引力只需浮点乘法，远远快于直接计算勒让德与三角函数。但需要注意，虽然对2到6阶的球谐引力场这样做是很合适的，但对高阶球谐场（比如32阶），这些多项式的每一项系数都巨大无比，求和后结果会大部分抵消，这会造成严重的浮点误差。而且时空复杂度均达到\(O(N^3)\)，其中\(N=\max(N_z,N_t)\)。使用勒让德和三角函数计算虽慢，但可以避免这些误差。（不过最好的方法是根本不去计算那么高阶的势场项，反正对结果贡献不大emmm……→_→）
————————更新（2024.1.28）
实际上那些勒让德和三角函数可以使用递推式逐级计算，巧妙实现后计算一次球谐引力也仅需要进行浮点乘除和有限次数的平方根，不需要直接调用超越函数。此方法在阶数\(n\)稍大时会快于完全展开成浮点乘法的形式，时间复杂度仅为\(O(N_z+N_t^2)\)，而使用的临时空间甚至可以降低至\(O(1+N_t)\)。
————————更新完毕
A2.转动惯量
直接对行星的引力势进行多极展开，将展开式里的几个积分对比转动惯量张量\(J\)的定义所用的几个积分，可得以下严格关系（读者自证不难）：\[\frac{J}{MR^2}=\frac{2}{3}(A-H_2)\]
其中\[A=\frac{1}{MR^2}\int_{\vec r\in V} r^2\rho(\vec r)dV\]
正是行星形变理论中定义的无量纲量\(A\)在非球对称行星中的扩展定义，即\((r/R)^2\)对行星质量的加权平均；\(M\)是行星质量；\(R\)是行星半径；\(H_2\)为二阶球谐势张量，定义见上一章。
行星的形变矩阵是零迹的，故形变不对\(A\)造成影响，\(A\)是一个由行星内部物质分布决定的常量。于是只需考虑形变对\(H_2\)的影响，求出\(H_2\)后即可直接算出转动惯量。
A3.行星的固有非球对称性
在本章我们将修改此前讨论的形变理论\(H_2=(k-1)C\)，使得修改后的理论能够描述行星固有的非球对称性。
首先展开\(H_2\)：
\[\begin{aligned}
H_2&=(k-1)C\\
&=(k-1)C_\omega+(k-1)k_tC_m\\
&=k_\omega C_\omega+k_mC_m
\end{aligned}
\]
其中，\(C_\omega=C(\vec\omega,f_\omega)\)与\(C_m=\sum_i C(\vec r_i,f_{m_i})\)分别为自转与潮汐力所对应的形变矩阵，定义与计算方法可见形变理论第二部分。实际计算中，\(C_m\)可以使用经过潮汐延迟的算法。\(k_\omega=k-1\)与\(k_m=(k-1)k_t\)为行星对自转和潮汐影响的形变反应系数（即Love Number）。
之前的形变理论认为行星总是对自转拥有完全的形变反应，即\(k_\omega\)无需类似\(k_m\)乘上一个\(k_t\)因子。这在实际星历积分中会导致问题：虽然类似地球的行星在长时间的稳定自转速度下，整体形变确实已接近完全平衡，但行星对这一平衡点附近自转速度的微小变化的响应不会是完全的。也就是说，自转导致的形变中有一部分应该被视为固定形变，剩下的部分才是随自转速度变化的。根据这一想法，可以将上式修改为\[H_2=H_{20}+k_\omega C_\omega+k_mC_m\]
此即形变理论的修改版本。其中\(H_{20}\)描述了行星的固有非球对称性；\(k_\omega\)不再等于\(k-1\)，而是与\(k_m\)类似，描述行星对自转速度的微小变化的形变响应。
修改后的形变理论中，描述行星形变性质的自由参数不再是\(k\)、\(A_r\)和\(k_t\)，而是\(A\)、\(k_\omega\)、\(k_m\)与\(H_{20}\)。\(H_{20}\)的具体取值自由度较高，一般不作为输入参数。而是根据天体实测的二阶球谐势\(H_2\)以及参数\(k_\omega\)、\(k_m\)的取值，从星历初始数据计算得到：\(H_{20}=H_2-(k_\omega C_\omega+k_mC_m)|_{t=0}\)。
总得来说这样修改之后形变理论就很实用了。即使对小行星和妊神星这样形状怪异的天体，我们都可以将其二阶球谐势中和原有形变理论不吻合的部分直接塞进\(H_{20}\)，非常方便~~~而高阶球谐势我们简单地假设不会随形变发生变化，直接按实测数值计算即可，反正影响不大。
最后一个问题是行星的形状，也就是行星海平面，或者叫参考椭球的形状。这里直接给出结论：行星的形状即为半径为\(R\)的球面上的点，经变换矩阵\[T=I+H_2+(C_\omega+C_m)\]作用后形成的椭球面。可以验证，在这一椭球面上，行星自身的质点势因距质心距离不同造成的差异，恰好可以抵消自转影响\(C_\omega\)加上潮汐影响\(C_m\)加上行星自身的球谐势场\(H_2\)的影响，使得总势能为一常值。当然，还是忽略高阶球谐势的影响。
~~~暂且完结~~~
目前我正在搞太阳系的星历积分，也就是写一个积分器，然后拟合一份初值，让积分器积出来的结果在尽可能长的过去未来中与实际星历的偏差尽可能小。这个工作量非常大，主要困难在于提高数值积分器的速度和精度。目前还没做完，不过已经有了一些初步结果：比如在1980~2020年范围内的月球轨道及自转偏差均可降低到2m以下；而在1800~2200年范围内，地轴偏差不超过10m，月球自转偏差不超过4m，轨道偏差不超过150m。我预想是搞出来一个包括太阳、八大行星和冥王星以及它们的73颗较大的卫星、343颗小行星带天体以及30颗柯依伯带天体的综合星历，覆盖上下两万年（BC 18000 ~ AD 22000）范围。
