首页 / 软件开发 / 数据结构与算法 / 算法系列（十八） 用天文方法计算二十四节气（下）

算法系列（十八） 用天文方法计算二十四节气（下）2014-05-24 csdn博客 吹泡泡的小猫【接上篇】

经过上述计算转换得到坐标值是理论值，或者说是天体的几何位置，但是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 * T2 + T3 / 189474        （ 3.10式）
太阳（地球）平近点角：
M = 357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300000         （3.11式）
月球平近点角
M"= 134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 56250         （3.12式）

月球纬度参数：
F = 93.27191 + 483202.017538 * T - 0.0036825 * T2 + T3 / 327270          （3.13式）
黄道与月球平轨道升交点黄经：
Ω= 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 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）{double T = dt * 10; /*T是从J2000起算的儒略世纪数*/ double T2 = T * T;double T3 = T2 * T;/*平距角（如月对地心的角距离）*/ *D = 297.85036 + 445267.111480 * T - 0.0019142 * T2 + T3 / 189474.0;/*太阳（地球）平近点角*/ *M = 357.52772 + 35999.050340 * T - 0.0001603 * T2 - T3 / 300000.0;/*月亮平近点角*/ *Mp = 134.96298 + 477198.867398 * T + 0.0086972 * T2 + T3 / 56250.0;/*月亮纬度参数*/ *F = 93.27191 + 483202.017538 * T - 0.0036825 * T2 + T3 / 327270.0;/*黄道与月亮平轨道升交点黄经*/ *Omega = 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 450000.0;}