4. 交点因子与天文幅角初值的逐日计算
为什么先算“中点时刻”的 f 与 (u+V)
在潮汐调和分析中，分潮的实际贡献可写成 h(t)=H·f·cos(V+u−g)。其中 f（交点因子）与 u（交点订正角）随月球升交点 N 的缓慢运动而变化，周期约 18.6 年；V 是天文幅角，随时间线性增长。工程上通常取分析时段中点时刻 t₀ 计算一组 f 与 (u+V)，作为整个时段的代表值。本节任务就是：给定 t₀，逐分潮算出 f 与 (u+V)，并检查其量级与相位是否合理。
计算所需的天文基本变量
以 2024-01-01 00:00（UTC）为起算历元，取分析时段中点 t₀ = 2024-07-01 00:00 UTC，即距历元 D = 182.0 天。按常用线性天文公式（单位：度）计算五个基本幅角：
s = 218.316 + 13.176396·D # 月球平黄经
h = 280.466 + 0.985647·D # 太阳平黄经
p = 83.353 + 0.111404·D # 月球近地点平黄经
N = 125.045 − 0.0529539·D # 月球升交点平黄经
p1= 282.938 + 0.0000471·D # 太阳近地点平黄经
D = 182.0 代入并取模 360：
s = 218.316 + 2398.104 = 2616.420 → 96.420°
h = 280.466 + 179.388 = 459.854 → 99.854°
p = 83.353 + 20.276 = 103.629°
N = 125.045 − 9.638 = 115.407°
p1 = 282.938 + 0.009 = 282.947°注意 N 是负向缓慢减小的量，2024 年 N 约在 115° 附近，属于“N 接近 90°~120°”的区间，此时各分潮的 f 会明显偏离 1，正是检验计算是否正确的良好样本。
交点因子 f 与交点订正角 u 的公式
对主要分潮，采用含 N 的经典近似式（角度单位度，结果取模 360）：
| 分潮 | f | u（度） |
|---|---|---|
| M2 | 1.000 − 0.037·cos N + 0.0002·cos 2N | −2.14·sin N |
| S2 | 1.000（无交点效应） | 0 |
| K1 | 1.006 + 0.115·cos N − 0.009·cos 2N + 0.002·cos 3N | −8.86·sin N + 0.68·sin 2N − 0.07·sin 3N |
| O1 | 1.009 + 0.187·cos N − 0.015·cos 2N | 10.80·sin N − 1.34·sin 2N + 0.19·sin 3N |
天文幅角 V 由 Doodson 数线性组合基本幅角得到，例如：
V(M2) = 2s − 2h
V(S2) = 0
V(K1) = s + h + 90°
V(O1) = −2s + h + 270°逐分潮计算结果
先算三角函数中间量：N=115.407°，sin N=0.9032，cos N=−0.4293，sin 2N=sin 230.814°=−0.7754，cos 2N=−0.6315，sin 3N=sin 346.221°=−0.2384，cos 3N=0.9712。
| 分潮 | f 计算 | f | u（度） | V（度） | (u+V) mod 360 |
|---|---|---|---|---|---|
| M2 | 1.000−0.037×(−0.4293)+0.0002×(−0.6315) | 1.0158 | −1.933 | 2×96.420−2×99.854=−6.868→353.132 | 351.199 |
| S2 | 恒为 1 | 1.0000 | 0 | 0 | 0.000 |
| K1 | 1.006+0.115×(−0.4293)−0.009×(−0.6315)+0.002×0.9712 | 0.9632 | −8.86×0.9032+0.68×(−0.7754)−0.07×(−0.2384)=−8.503 | 96.420+99.854+90=286.274 | 277.771 |
| O1 | 1.009+0.187×(−0.4293)−0.015×(−0.6315) | 1.0382 | 10.80×0.9032−1.34×(−0.7754)+0.19×(−0.2384)=10.756 | −2×96.420+99.854+270=177.014 | 187.770 |
合理性检查：量级与相位
算完不能直接使用，需做三项检查：
- f 的量级：M2 的 f 应在 0.96~1.04 之间，K1 在 0.88~1.12 之间，O1 在 0.80~1.20 之间。本例 M2=1.0158、K1=0.9632、O1=1.0382，均落在合理区间；若某分潮 f 出现 0.5 或 1.5 这类值，多半是 N 未取模或三角函数用了弧度。
- u 的量级：M2 的 |u| 不超过约 2.2°，K1 不超过约 9.6°，O1 不超过约 12.3°。本例 M2=−1.933°、K1=−8.503°、O1=10.756°，符合预期。
- (u+V) 的相位：S2 的 V 恒为 0，是天然的“零点基准”；M2 的 V=2s−2h 应等于 2×(s−h)，本例 s−h=−3.434°，2 倍为 −6.868°，取模后 353.132°，与表中一致。若 M2 的 (u+V) 与 S2 的差不是约 350° 量级，说明 s、h 的历元或系数写错。
可复用的计算脚本
把上述过程写成 Python 函数，输入中点时刻的儒略日偏移 D，输出各分潮的 f 与 (u+V)，并内置量级检查：
import math
def norm360(x):
return x % 360.0
def astro_args(D):
s = norm360(218.316 + 13.176396 * D)
h = norm360(280.466 + 0.985647 * D)
p = norm360( 83.353 + 0.111404 * D)
N = norm360(125.045 - 0.0529539 * D)
p1 = norm360(282.938 + 0.0000471 * D)
return s, h, p, N, p1
def f_u_V(D):
s, h, p, N, p1 = astro_args(D)
Nr = math.radians(N)
sinN, cosN = math.sin(Nr), math.cos(Nr)
sin2N, cos2N = math.sin(2*Nr), math.cos(2*Nr)
sin3N, cos3N = math.sin(3*Nr), math.cos(3*Nr)
res = {}
# M2
f = 1.000 - 0.037*cosN + 0.0002*cos2N
u = -2.14*sinN
V = norm360(2*s - 2*h)
res['M2'] = (f, norm360(u + V))
# S2
res['S2'] = (1.0, 0.0)
# K1
f = 1.006 + 0.115*cosN - 0.009*cos2N + 0.002*cos3N
u = -8.86*sinN + 0.68*sin2N - 0.07*sin3N
V = norm360(s + h + 90.0)
res['K1'] = (f, norm360(u + V))
# O1
f = 1.009 + 0.187*cosN - 0.015*cos2N
u = 10.80*sinN - 1.34*sin2N + 0.19*sin3N
V = norm360(-2*s + h + 270.0)
res['O1'] = (f, norm360(u + V))
return res
def check(res):
limits = {'M2': (0.96, 1.04), 'S2': (0.999, 1.001),
'K1': (0.88, 1.12), 'O1': (0.80, 1.20)}
for name, (f, uv) in res.items():
lo, hi = limits[name]
ok = lo <= f <= hi
print(f"{name}: f={f:.4f} (u+V)={uv:.3f}° f量级{'正常' if ok else '异常'}")
if __name__ == '__main__':
r = f_u_V(182.0)
check(r)运行后应输出与上表一致的数值：M2 f=1.0158、(u+V)=351.199°；S2 f=1.0000、(u+V)=0.000°；K1 f=0.9632、(u+V)=277.771°；O1 f=1.0382、(u+V)=187.770°，且四项均判为“正常”。
结果的使用方式
得到的 f 与 (u+V) 将作为常数进入设计矩阵：对每个时刻 t，第 j 个分潮的列向量为 f_j·cos(ω_j·t + (u+V)_j) 与 f_j·sin(ω_j·t + (u+V)_j)，其中 ω_j 为分潮角速度。这样最小二乘解出的系数再经反算即可得到振幅 H 与迟角 g。若时段跨越数年，应改用逐日或逐年变化的 f、u，而非单一中点值。