上篇：日月引潮力探究 讲述了引潮力原理，此篇再次简要说明并着重于计算。

引潮力原理

任意双星稳定系统（稳定是指运行轨道相对稳定，不明显变化，如地月系统），各星球绕两者的质量中心运行。根据万有引力公式，两者之间的引力为：

$$ F = G\frac{M_AM_B}{r^2} $$

其中G为引力常数，为$6.67259\times10^{-11}\mbox{N}\cdot\mbox{m^2}\cdot\mbox{kg}^{-2}$，$M_A, M_B$分别为两个星球的质量，r为AB两星质心距离。在A星质心处，向心力由B星提供，B星的万有引力与A星质点处的惯性离心力相互抵消，且A星上任意一点处受到的惯性离心力大小方向都和质心处相同。因此，A星上任意一处单位质量（下文均使用单位质量）受到的惯性离心力$\overrightarrow{f}$与质心处受到的引力相等。

$$ \overrightarrow{f} = G\frac{M_B}{r^3}\overrightarrow{r} $$

设A星上任一点到B星质心距离为$d$，那么该点受到的B星对它的万有引力$\overrightarrow{g}$为

$$ \overrightarrow{g} = G\frac{M_B}{d^3}\overrightarrow{d} $$

则该点受到的引潮力，为惯性离心力和万有引力的合力：

$$ \overrightarrow{F} = \overrightarrow{f}+\overrightarrow{g} $$

地月系下引潮力计算

以地月系为例，讨论引潮力的计算：

A为地球，C为月球，r为地月质心距离，B为地月质心连线在地表的交点，称为月下点，R为地球直径。$\alpha$ 为点P、地球质心与月下点的夹角。因引潮力沿地月质心连线轴对称分布，因此可用$\alpha$表示点P位置，$\alpha \in [0,\pi]$。

引潮力大小

首先惯性离心力为月球对地球质心的单位质量的万有引力，方向为从月质心指向地质心，大小为。

$$ f = G\frac{M_m}{r^2} $$

其垂直分量与水平分量分别为：

$$ f_h = -G\frac{M_m}{r^2}\sin\alpha\\ f_v = -G\frac{M_m}{r^2}\cos\alpha $$

符号为定义力的方向，垂直分量以向上为正，水平分量以沿球面指向月下点方向为正。

点P到月质心的距离d，通过余弦定理得$d^2=r^2+R^2-2Rr\cos\alpha$，则月球对点P处引力大小为：

$$ g = G\frac{M_m}{d^2}=G\frac{M_m}{r^2+R^2+2Rr\cos\alpha} $$

根据正弦定理和余弦定律，得

$$ \begin{align} &\sin\beta=\frac{r\sin\alpha}{d}\\ &\cos\beta=\frac{d^2+R^2-r^2}{2dR} \end{align} $$

因此，引力的垂直分量和水平分量为：

$$ \begin{align} &g_h = G\frac{M_m}{d^2}\sin\beta =G\frac{M_mr\sin\alpha}{d^3} \\ &g_v = -G\frac{M_m}{d^2}\cos\beta =-G\frac{M_m(d^2+R^2-r^2)}{2d^3R} \end{align} $$

所以，引潮力的垂直分量和水平分量为：

$$ \begin{align} &F_h = f_h+g_h=GM_m\sin\alpha\left(\frac{r}{d^3}-\frac{1}{r^2}\right)\\ &F_v = f_v+g_v = -GM\left(\frac{d^2+R^2-r^2}{2d^3R}+\frac{\cos\alpha}{r^2}\right) \end{align} $$

至此，计算引潮力大小的公式已给出，只需要知道$\alpha, R, r, M_m$ 即可计算引潮力。下面给出$\alpha$和引潮力方向的计算。

引潮力方向

首先需要计算点P与月下点的夹角$\alpha$，使用两点间球心角公式：

$$ \alpha = \cos\varphi\cos\varphi_m\cos(\lambda-\lambda_m)+\sin\varphi\sin\varphi_m $$

其中$\varphi$为点P纬度，$\lambda$为点P经度；$\varphi_m,\lambda_m$分别为月下点的纬度和经度。

下一步是判断引潮力方向，垂直引潮力正值为上，负值为下；水平引潮力则正值指向月下点，负值指向地月质心连线与地表相交的远点，即与月下点地心对称的点。水平引潮力的方向，使用方位角表示，即0°为北、90°为东、180°为南、270°为西，因此需要计算月下点（或其对称点）对于点P的方位角。
方位角计算公式引自Calculate distance, bearing and more between Latitude/Longitude points，方位角$\theta$为：

$$ \theta = \mbox{atan2}(\sin(\lambda-\lambda_m)\cdot\cos\varphi_m , \cos\varphi\cdot\sin\varphi_m − \sin\varphi\cdot\cos\varphi_m\cdot\cos(\lambda-\lambda_m)) $$

其中atan2为计算方位角的函数。

地月、地日双系统下的引潮力

地日系的引潮力与地月系相同，计算得到月球与太阳的引潮力后，两者相加即得到日月共同形成影响的引潮力。
垂直方向上，直击相加即可；水平方向上，因为日月两个水平引潮力方向不同，需要进行力的合成计算。
以下是Python程序示例，输入日下点、月下点、地表点的经纬度和日月距离（天文单位），返回该情况下的瞬时垂直、水平引潮力和方向，需注意角度和弧度的转换：

import numpy as np
def tidal_force(sun_lon, sun_lat, moon_lon, moon_lat, lon, lat, sun_r, moon_r):
    R = 6378137  # m
    moon_m = 7.349e22  # kg
    sun_m = 1.9891e30
    G = 6.67259e-11
    au = 149597870700

    sun_lon = sun_lon * np.pi / 180
    sun_lat = sun_lat * np.pi / 180
    moon_lon = moon_lon * np.pi / 180
    moon_lat = moon_lat * np.pi / 180
    lon = lon * np.pi / 180
    lat = lat * np.pi / 180
    sun_r = sun_r * au
    moon_r = moon_r * au

    # 月引潮力
    alpha = np.arccos(np.cos(moon_lat) * np.cos(lat) * np.cos(moon_lon - lon) + np.sin(moon_lat) * np.sin(lat))
    dd = moon_r ** 2 + R ** 2 - 2 * R * moon_r * np.cos(alpha)
    d = np.sqrt(dd)
    fh_moon = G * moon_m * np.sin(alpha) * (
                1 / (moon_r + R ** 2 / moon_r - 2 * R * np.cos(alpha)) * (1 / d) - (1 / moon_r ** 2))
    fv_moon = -G * moon_m * (((dd + R ** 2 - moon_r ** 2) / (2 * d * R)) * (
                1 / (moon_r ** 2 + R ** 2 - 2 * R * moon_r * np.cos(alpha))) +
                             np.cos(alpha) / moon_r ** 2)
    # 月水平引潮力方位角
    if fh_moon < 0:
        moon_lon = (moon_lon + np.pi) % (2 * np.pi)
        moon_lat *= -1
    theta_moon = np.arctan2(np.sin(moon_lon - lon) * np.cos(moon_lat), np.cos(lat) * np.sin(moon_lat) - np.sin(lat) *
                            np.cos(moon_lat) * np.cos(moon_lon - lon))
    theta_moon = (theta_moon * 180 / np.pi + 360) % 360

    # 日引潮力
    alpha = np.arccos(np.cos(sun_lat) * np.cos(lat) * np.cos(sun_lon - lon) + np.sin(sun_lat) * np.sin(lat))
    dd = sun_r ** 2 + R ** 2 - 2 * R * sun_r * np.cos(alpha)
    d = np.sqrt(dd)
    fh_sun = G * sun_m * np.sin(alpha) * (
            1 / (sun_r + R ** 2 / sun_r - 2 * R * np.cos(alpha)) * (1 / d) - (1 / sun_r ** 2))
    fv_sun = -G * sun_m * (((dd + R ** 2 - sun_r ** 2) / (2 * d * R)) * (
            1 / (sun_r ** 2 + R ** 2 - 2 * R * sun_r * np.cos(alpha))) +
                           np.cos(alpha) / sun_r ** 2)
    # 日水平引潮力方位角
    if fh_sun < 0:
        sun_lon = (sun_lon + np.pi) % (2 * np.pi)
        sun_lat *= -1
    theta_sun = np.arctan2(np.sin(sun_lon - lon) * np.cos(sun_lat), np.cos(lat) * np.sin(sun_lat) - np.sin(lat) *
                           np.cos(sun_lat) * np.cos(sun_lon - lon))
    theta_sun = (theta_sun * 180 / np.pi + 360) % 360
    
    # 力合成
    fv = fv_sun + fv_moon
    x = fh_sun * np.sin(theta_sun * np.pi / 180) + fh_moon * np.sin(theta_moon * np.pi / 180)
    y = fh_sun * np.cos(theta_sun * np.pi / 180) + fh_moon * np.cos(theta_moon * np.pi / 180)
    fh = np.sqrt(x ** 2 + y ** 2)
    theta = np.angle(np.complex64(x) + np.complex64(y) * np.complex64(1j), deg=True)
    if 0 <= theta <= 90:
        theta = 90 - theta
    elif theta < 0:
        theta = 90 - theta
    else:
        theta = 360 - (theta - 90)
    return fv, fh, theta

模型检验

根据上面给出的程序，画出大潮（日下点0E, 0N、月下点0E, 0N）和小潮（日下点0E, 0N、月下点90E, 0N）时，赤道上不同经度引潮力情况：

图中蓝色线为大潮、红色线为小潮，数值变化符合预期，数量级也符合预期。