算法系列之十八：用天文方法计算二十四节气（下）

【接上篇】

 

        经过上述计算转换得到坐标值是理论值，或者说是天体的几何位置，但是FK5系统是一个目视系统，也就是说体现的是人眼睛观察效果（光学位置），这就需要根据地球的物理环境、大气环境等信息做进一步的修正，使其和人类从地球上观察星体的观测结果一致。

        首先需要进行章动修正。章动是指地球沿自转轴的指向绕黄道极缓慢旋转过程中，由于地球上物质分布不均匀性和月球及其它行星的摄动力造成的轻微抖动。英国天文学家詹姆斯·布拉德利(1693—1762)最早发现了章动，章动可以沿着黄道分解为水平分量和垂直分量，黄道上的水平分量记为Δψ，称为黄经章动，它影响了天球上所有天体的经度。黄道上的垂直分量记为Δε，称为交角章动，它影响了黄赤交角。目前编制天文年历所依据的章动理论是伍拉德在1953年建立的，它是以刚体地球模型为基础的。1977年，国际天文联合会的一个专家小组建议采用非刚体地球模型――莫洛坚斯基II模型代替刚体地球模型计算章动，1979年的国际天文学联合会第十七届大会正式通过了这一建议，并决定于1984年正式实施。

        地球章动主要是月球运动引起的，也具有一定的周期性，可以描述为一些周期项的和，主要项的周期是6798.4日(18.6年)，但其它项是一些短周期项(小于10天)。本文采用的计算方法取自国际天文联合会的IAU1980章动理论，周期项系数数据来源于《天文算法》一书第21章的表21－A，该表忽略了IAU1980章动理论中系数小于0.0003"的周期项，因此只有63项。每个周期项包括计算黄经章动（Δψ）的正弦系数（相位内项系数）、计算交角章动的（Δε）余弦系数（相位外项系数）以及计算辐角的5个基本角距（M、M'、D、F、Ω）的线性组合系数。5个基本角距的计算公式是：

 

平距角(日月对地心的角距离)：
D = 297.85036 + 455267.111480 * T - 0.0019142 * T² + T³ / 189474        （3.10式）
太阳（地球）平近点角：
M = 357.52772 + 35999.050340 * T - 0.0001603 * T² - T³ / 300000         （3.11式）
月球平近点角
M'= 134.96298 + 477198.867398 * T + 0.0086972 * T² + T³ / 56250        （3.12式）

月球纬度参数：
F = 93.27191 + 483202.017538 * T - 0.0036825 * T² + T³ / 327270          （3.13式）
黄道与月球平轨道升交点黄经：
Ω= 125.04452 - 1934.136261 * T + 0.0020708 * T² + T³ / 450000            （3.14式）

 

以上各式中的T是儒略世纪数，计算出来的5个基本角距的单位都是度，在计算正弦或余弦时要转换为弧度单位。计算每一个周期项的黄经章动过程是这样的，首先将3.10－3.14式计算出来的值与对应的5个基本角距系数组合，计算出辐角。以本文使用的章动周期项系数表中的第七项为例，5个基本角距对应的系数分别是1、0、-2、2和2，辐角θ的值就是：-2D + M + 2F + 2Ω。计算出辐角后就可以计算周期项的值：

 

S = (S1+ S2 * T) * sin(θ)                          （3.15式）

 

仍以第七项为例，S的值就是（-517 + 1.2 * T）* sin(θ)。对每一项的值S累加就可得到黄经章动，单位是0.0001"。交角章动的计算方法与黄经章动的计算类似，辐角θ的值是一样的，只是计算章动使用的是余弦系数：

C = (C1 + C2 * T) * cos(θ)                          （3.16式）

 

CalcEarthLongitudeNutation()函数就是计算黄经章动的实现代码：

double CalcEarthLongitudeNutation(double dt)  {      double T = dt * 10;      double D,M,Mp,F,Omega;        GetEarthNutationParameter(dt, &D, &M, &Mp, &F, &Omega);        double resulte = 0.0 ;      for(int i = 0; i < sizeof(nutation) / sizeof(nutation[0]); i++)      {          double sita = nutation[i].D * D + nutation[i].M * M + nutation[i].Mp * Mp + nutation[i].F * F + nutation[i].omega * Omega;                    resulte += (nutation[i].sine1 + nutation[i].sine2 * T ) * sin(sita);      }        /*先乘以章动表的系数 0.0001，然后换算成度的单位*/      return resulte * 0.0001 / 3600.0;  }

 GetEarthNutationParameter()辅助函数用于计算5个基本角距：

void GetEarthNutationParameter(double dt, double *D, double *M, double *Mp, double *F, double *Omega)