万年历计算之节气

一、基本知识
二十四节气起源于黄河流域。远在春秋时代，就定出仲春、仲夏、仲秋和仲冬等四个节气。以后不断地改进与完善，到秦汉年间，二十四节气已完全确立。公元前104年，由邓平等制定的《太初历》，正式把二十四节气订于历法，明确了二十四节气的天文位置。
　　太阳从黄经零度起，沿黄经每运行15度所经历的时日称为“一个节气”。每年运行360度，共经历24个节气，每月2个。其中，每月第一个节气为“节气”，即：立春、惊蛰、清明、立夏、芒种、小暑、立秋、白露、寒露、立冬、大雪和小寒等12个节气；每月的第二个节气为“中气”，即：雨水、春分、谷雨、小满、夏至、大暑、处暑、秋分、霜降、小雪、冬至和大寒等12个节气。“节气” 和“中气”交替出现，各历时15天，现在人们已经把“节气”和“中气”统称为“节气”。
　　二十四节气反映了太阳的周年视运动，所以节气在现行的公历中日期基本固定，上半年在6日、21日，下半年在8日、23日，前后不差1～2天。
二、基本概念
儒略日：Julian Day number 或 Julian Day 或 JD。是指从公元前 4712 年开始连续计算日数得出的天数及不满一日的尾数。传统上儒略日的计数是从格林尼治平午，即世界时间 12 点开始的。
格里历：即公历或格列高利历，是现行国际通行的历法，属于阳历的一种，通称阳历，其前身是奥古斯都历，而奥古斯都历的前身是儒略历。其历年为一个回归年(365.2425日)，划分为12个历月。是教皇格里高利13世(也译格雷果里)在公元1582年改革儒略历制定的历法，儒略历的回归年为365.25，与实际的回归365.2422相差甚多，当时儒略历和地球实际位置的误差已达14天，格里历将误差纠正，确定所有整数世纪年除了可被400整除的外一律不设闰年，同时规定1582年10月4日之后的那天为1582年10月15日，但原有星期不变。新颁布的历法理论上可以达到两万年内误差不超过一天，但由于地球自转的变化，实际到公元4909年误差就可达一天。与儒略历相比，公历是新制定的历法，所以有时候，包括中华民国《中国国家标准》CNS 7648《数据元及交换格式–信息交换–日期及时间的表示法》，又称新历。
黄经：地球绕太阳运转一周约365天5小时多，运转94，000万公里。地球的公转在地球上人看来就表现为太阳周年视运动，其运行线路（即地球公转轨道在天球上的反映）称为黄道。黄经就是黄道上的度量坐标（经度），按天文学惯例，以春分点为起点自西向东度量，分360度。我国古人把太阳黄经的360度划分成24等份，每份15度，为一个节气。两个节气间相隔日数为15天左右，全年即有二十四节气。二十四节气的太阳黄经度分别为：立春315度，雨水330度，惊蛰345度，春分360度，清明15度，谷雨30度，立夏45度，小满60度，芒种75度，夏至90度，小暑105度，大暑120度，立秋135度，处暑150度，白露165度，秋分180度，寒露195度，霜降210度，立冬225度，小雪240度，大雪255度，冬至270度，小寒285度，大寒300度。
黄纬：黄道上的纬度。
三、基本算法：
3.1：儒略日的计算（算法描述请参考相关书籍，为了节省空间和时间，这里请罗列出 C++ 实现代码）

// ------------------------------------------------------------------------
//
// 儒略日转换为格历日
//
// dbGD[IN]：儒略日
// year[OUT]：年
// month[OUT]：月
// day[OUT]：日
// hour[OUT]：时
// minute[OUT]：分
// second[OUT]：秒
// ------------------------------------------------------------------------
void w_JDToGD(const double &dbJD, int &year, int &month, int &day, int &hour, int &minute, int &second)
{
double dJDM = dbJD + 0.5;
unsigned long ulZ = static_cast<unsigned long>(floor(dJDM));
double dF = dJDM - floor(dJDM);
unsigned long ulA, ulB, ulC, ulD;
int nE, nQ;
if (dbJD < 2299161)
ulA = ulZ;
else
{
nQ = static_cast<int>((ulZ - 1867216.25) / 36524.25);
ulA = ulZ + 1 + nQ - static_cast<int>(nQ >> 2);
}
ulB = ulA + 1524 ;
ulC = static_cast<int>(( ulB - 122.1) / 365.25);
ulD = static_cast<int>(365.25 * ulC);
nE = static_cast<int>((ulB - ulD) / 30.6001);
// 计算日
day = ulB - ulD - int ( 30.6001 * nE ) ;
// 计算时
hour = static_cast<int>(floor(dF * 24.0));
// 计算分
minute = static_cast<int>(((dF * 24.0) - floor(dF * 24.0)) * 60.0);
// 计算秒
second = static_cast<int>((((dF * 24.0) * 60.0) - floor(( dF * 24.0) * 60.0)) * 60.0);
// 计算月
if(nE < 14)
month = nE - 1;
if(nE == 14 || nE == 15)
month = nE - 13;
// 计算年
if(month > 2)
year = ulC - 4716;
if(month == 1 || month == 2)
year = ulC - 4715;
}

// ------------------------------------------------------------------------
//
// 格历日转换为儒略日
//
// year[IN]：年
// month[IN]：月
// day[IN]：日
// hour[IN]：时
// minute[IN]：分
// second[IN]：秒
//
// 返回对应的儒略日
// ------------------------------------------------------------------------
double w_GDToJD(const int &year, const int &month, const int &day, const int &hour, const int &minute, const int &second)
{
int nY = year, nM = month;
double dD = static_cast<double>(day) + hour / 24.0 + (minute / 60.0) / 24.0 + ((second / 60.0) / 60.0) / 24.0;
if(month == 1 || month == 2)
{
nY = year - 1;
nM = month + 12;
}
int nA = nY / 100;
int nB = 2 - nA + (nA >> 2);
return static_cast<double>(static_cast<int>(365.25 * (nY + 4716)) + static_cast<int>(30.6001 * (nM + 1)) + dD + nB - 1524.5);
}
3.2、角度调整
// ------------------------------------------------------------------------
//
// 调整角度到 0-360 之间
//
// dbDegrees[IN]：角度
//
// 返回调整后的角度
// ------------------------------------------------------------------------
double w_MapTo0To360Range(const double &dbDegrees)
{
double dbValue = dbDegrees;
// map it to the range 0 - 360
while(dbValue < 0.0)
dbValue += 360.0;
while(dbValue > 360.0)
dbValue -= 360.0;
return dbValue;
}
四、步骤简述
4.1、按照 5.2 和 5.3 节介绍方法计算出某一时刻太阳黄经和黄纬。
4.2、按照 5.4 、5.5、5.6 节介绍的三次校正方法对黄经进行校正，即可求得某时刻太阳黄经度数。
4.3、综合步骤演示在第 六 节。
4.4、第 七 节演示如何使用二分法求指定节气的时间。
4.5、第 八 节介绍了 本地时间（LST）和格林尼治时间（UTC）之间互换的函数。
4.6、第 九 节罗列出了本代码计算出的 2008年24节气时间（北京时间）
五、太阳黄经、黄纬、向量半径以及校正量的计算
5.1、为了代码可读性，定义类型变量。
//
// 圆周率
const double PI = 3.1415926535897932384626433832795;
//
// 一度代表的弧度
const double dbUnitRadian = PI / 180.0;
//
// 计算太阳黄经赤纬所需类型变量
typedef struct
{
double dA;
double dB;
double dC;
} VSOP87COEFFICIENT, *PVSOP87COEFFICIENT
// 二次修正黄经赤纬所需的天体章动系数类型变量
typedef struct
{
int nD;
int nM;
int nMprime;
int nF;
int nOmega;
int nSincoeff1;
double dSincoeff2;
int nCoscoeff1;
double dCoscoeff2;
} NUTATIONCOEFFICIENT, *PNUTATIONCOEFFICIENT;
//
// 节气枚举
enum SOLARTERMS
{
ST_VERNAL_EQUINOX = 0, // 春分
ST_CLEAR_AND_BRIGHT = 15, // 清明
ST_GRAIN_RAIN = 30, // 谷雨
ST_SUMMER_BEGINS = 45, // 立夏
ST_GRAIN_BUDS = 60, // 小满
ST_GRAIN_IN_EAR = 75, // 芒种
ST_SUMMER_SOLSTICE = 90, // 夏至
ST_SLIGHT_HEAT = 105, // 小暑
ST_GREAT_HEAT = 120, // 大暑
ST_AUTUMN_BEGINS = 135, // 立秋
ST_STOPPING_THE_HEAT = 150, // 处暑
ST_WHITE_DEWS = 165, // 白露
ST_AUTUMN_EQUINOX = 180, // 秋分
ST_COLD_DEWS = 195, // 寒露
ST_HOAR_FROST_FALLS = 210, // 霜降
ST_WINTER_BEGINS = 225, // 立冬
ST_LIGHT_SNOW = 240, // 小雪
ST_HEAVY_SNOW = 255, // 大雪
ST_WINTER_SOLSTICE = 270, // 冬至
ST_SLIGHT_COLD = 285, // 小寒
ST_GREAT_COLD = 300, // 大寒
ST_SPRING_BEGINS = 315, // 立春
ST_THE_RAINS = 330, // 雨水
ST_INSECTS_AWAKEN = 345 // 惊蛰
};
5.2、黄经：（算法描述请参考相关书籍，为了节省空间和时间，这里请罗列出 C++ 实现代码）
//
// 计算太阳黄经用参数
const VSOP87COEFFICIENT Earth_SLG0[64] =
{
{ 175347046.0 , 0.0000000 , 000000.0000000 } ,
{ 3341656.0 , 4.6692568 , 6283.0758500 } ,
{ 34894.0 , 4.6261000 , 12566.1517000 } ,
{ 3497.0 , 2.7441000 , 5753.3849000 } ,
{ 3418.0 , 2.8289000 , 3.5231000 } ,
{ 3136.0 , 3.6277000 , 77713.7715000 } ,
{ 2676.0 , 4.4181000 , 7860.4194000 } ,
{ 2343.0 , 6.1352000 , 3930.2097000 } ,
{ 1324.0 , 0.7425000 , 11506.7698000 } ,
{ 1273.0 , 2.0371000 , 529.6910000 } ,
{ 1199.0 , 1.1096000 , 1577.3435000 } ,
{ 990.0 , 5.2330000 , 5884.9270000 } ,
{ 902.0 , 2.0450000 , 26.2980000 } ,
{ 857.0 , 3.5080000 , 398.1490000 } ,
{ 780.0 , 1.1790000 , 5223.6940000 } ,
{ 753.0 , 2.5330000 , 5507.5530000 } ,
{ 505.0 , 4.5830000 , 18849.2280000 } ,
{ 492.0 , 4.2050000 , 775.5230000 } ,
{ 357.0 , 2.9200000 , 000000.0670000 } ,
{ 317.0 , 5.8490000 , 11790.6290000 } ,
{ 284.0 , 1.8990000 , 796.2880000 } ,
{ 271.0 , 0.3150000 , 10977.0790000 } ,
{ 243.0 , 0.3450000 , 5486.7780000 } ,
{ 206.0 , 4.8060000 , 2544.3140000 } ,
{ 205.0 , 1.8690000 , 5573.1430000 } ,
{ 202.0 , 2.4580000 , 6069.7770000 } ,
{ 156.0 , 0.8330000 , 213.2990000 } ,
{ 132.0 , 3.4110000 , 2942.4630000 } ,
{ 126.0 , 1.0830000 , 20.7750000 } ,
{ 115.0 , 0.6450000 , 000000.9800000 } ,
{ 103.0 , 0.6360000 , 4694.0030000 } ,
{ 102.0 , 0.9760000 , 15720.8390000 } ,
{ 102.0 , 4.2670000 , 7.1140000 } ,
{ 99.0 , 6.2100000 , 2146.1700000 } ,
{ 98.0 , 0.6800000 , 155.4200000 } ,
{ 86.0 , 5.9800000 , 161000.6900000 } ,
{ 85.0 , 1.3000000 , 6275.9600000 } ,
{ 85.0 , 3.6700000 , 71430.7000000 } ,
{ 80.0 , 1.8100000 , 17260.1500000 } ,
{ 79.0 , 3.0400000 , 12036.4600000 } ,
{ 75.0 , 1.7600000 , 5088.6300000 } ,
{ 74.0 , 3.5000000 , 3154.6900000 } ,
{ 74.0 , 4.6800000 , 801.8200000 } ,
{ 70.0 , 0.8300000 , 9437.7600000 } ,
{ 62.0 , 3.9800000 , 8827.3900000 } ,
{ 61.0 , 1.8200000 , 7084.9000000 } ,
{ 57.0 , 2.7800000 , 6286.6000000 } ,
{ 56.0 , 4.3900000 , 14143.5000000 } ,
{ 56.0 , 3.4700000 , 6279.5500000 } ,
{ 52.0 , 0.1900000 , 12139.5500000 } ,
{ 52.0 , 1.3300000 , 1748.0200000 } ,
{ 51.0 , 0.2800000 , 5856.4800000 } ,
{ 49.0 , 0.4900000 , 1194.4500000 } ,
{ 41.0 , 5.3700000 , 8429.2400000 } ,
{ 41.0 , 2.4000000 , 19651.0500000 } ,
{ 39.0 , 6.1700000 , 10447.3900000 } ,
{ 37.0 , 6.0400000 , 10213.2900000 } ,
{ 37.0 , 2.5700000 , 1059.3800000 } ,
{ 36.0 , 1.7100000 , 2352.8700000 } ,
{ 36.0 , 1.7800000 , 6812.7700000 } ,
{ 33.0 , 0.5900000 , 17789.8500000 } ,
{ 30.0 , 0.4400000 , 83996.8500000 } ,
{ 30.0 , 2.7400000 , 1349.8700000 } ,
{ 25.0 , 3.1600000 , 4690.4800000 }
};
const VSOP87COEFFICIENT Earth_SLG1[34] =
{
{ 628331966747.0 , 0.000000 , 00000.0000000 } ,
{ 206059.0 , 2.678235 , 6283.0758500 } ,
{ 4303.0 , 2.635100 , 12566.1517000 } ,
{ 425.0 , 1.590000 , 3.5230000 } ,
{ 119.0 , 5.796000 , 26.2980000 } ,
{ 109.0 , 2.966000 , 1577.3440000 } ,
{ 93.0 , 2.590000 , 18849.2300000 } ,
{ 72.0 , 1.140000 , 529.6900000 } ,
{ 68.0 , 1.870000 , 398.1500000 } ,
{ 67.0 , 4.410000 , 5507.5500000 } ,
{ 59.0 , 2.890000 , 5223.6900000 } ,
{ 56.0 , 2.170000 , 155.4200000 } ,
{ 45.0 , 0.400000 , 796.3000000 } ,
{ 36.0 , 0.470000 , 775.5200000 } ,
{ 29.0 , 2.650000 , 7.1100000 } ,
{ 21.0 , 5.430000 , 00000.9800000 } ,
{ 19.0 , 1.850000 , 5486.7800000 } ,
{ 19.0 , 4.970000 , 213.3000000 } ,
{ 17.0 , 2.990000 , 6275.9600000 } ,
{ 16.0 , 0.030000 , 2544.3100000 } ,
{ 16.0 , 1.430000 , 2146.1700000 } ,
{ 15.0 , 1.210000 , 10977.0800000 } ,
{ 12.0 , 2.830000 , 1748.0200000 } ,
{ 12.0 , 3.260000 , 5088.6300000 } ,
{ 12.0 , 5.270000 , 1194.4500000 } ,
{ 12.0 , 2.080000 , 4694.0000000 } ,
{ 11.0 , 0.770000 , 553.5700000 } ,
{ 10.0 , 1.300000 , 6286.6000000 } ,
{ 10.0 , 4.240000 , 1349.8700000 } ,
{ 9.0 , 2.700000 , 242.7300000 } ,
{ 9.0 , 5.640000 , 951.7200000 } ,
{ 8.0 , 5.300000 , 2352.8700000 } ,
{ 6.0 , 2.650000 , 9437.7600000 } ,
{ 6.0 , 4.670000 , 4690.4800000 }
};
const VSOP87COEFFICIENT Earth_SLG2[20] =
{
{ 52919.0 , 0.0000 , 00000.0000 } ,
{ 8720.0 , 1.0721 , 6283.0758 } ,
{ 309.0 , 0.8670 , 12566.1520 } ,
{ 27.0 , 0.0500 , 3.5200 } ,
{ 16.0 , 5.1900 , 26.3000 } ,
{ 16.0 , 3.6800 , 155.4200 } ,
{ 10.0 , 0.7600 , 18849.2300 } ,
{ 9.0 , 2.0600 , 77713.7700 } ,
{ 7.0 , 0.8300 , 775.5200 } ,
{ 5.0 , 4.6600 , 1577.3400 } ,
{ 4.0 , 1.0300 , 7.1100 } ,
{ 4.0 , 3.4400 , 5573.1400 } ,
{ 3.0 , 5.1400 , 796.3000 } ,
{ 3.0 , 6.0500 , 5507.5500 } ,
{ 3.0 , 1.1900 , 242.7300 } ,
{ 3.0 , 6.1200 , 529.6900 } ,
{ 3.0 , 0.3100 , 398.1500 } ,
{ 3.0 , 2.2800 , 553.5700 } ,
{ 2.0 , 4.3800 , 5223.6900 } ,
{ 2.0 , 3.7500 , 00000.9800 }
};
const VSOP87COEFFICIENT Earth_SLG3[ 7] =
{
{ 289.0 , 5.844 , 6283.076 } ,
{ 35.0 , 0.000 , 00000.000 } ,
{ 17.0 , 5.490 , 12566.150 } ,
{ 3.0 , 5.200 , 155.420 } ,
{ 1.0 , 4.720 , 3.520 } ,
{ 1.0 , 5.300 , 18849.230 } ,
{ 1.0 , 5.970 , 242.730 }
};
const VSOP87COEFFICIENT Earth_SLG4[ 3] =
{
{ 114.0 , 3.142 , 00000.00 } ,
{ 8.0 , 4.130 , 6283.08 } ,
{ 1.0 , 3.840 , 12566.15 }
};
const VSOP87COEFFICIENT Earth_SLG5[ 1] =
{
{ 1.0 , 3.14 , 0.0 }
};
// ------------------------------------------------------------------------
//
// 计算太阳在黄道面上的经度（单位：度）
//
// dbJD[IN]：儒略日（计算该时刻太阳在黄道面上的经度）
//
// 返回太阳黄经度数
// ------------------------------------------------------------------------
double w_GetSunLongitude(const double & dbJD)
{
// 计算τ
double dt = (dbJD - 2451545.0) / 365250.0;
double dL = 0.0, dL0 = 0.0, dL1 = 0.0, dL2 = 0.0, dL3 = 0.0, dL4 = 0.0, dL5 = 0.0;
// L0 38x3
for(int i = 0; i < sizeof(Earth_SLG0) / sizeof(VSOP87COEFFICIENT); i++)
dL0 += (Earth_SLG0[i].dA * cos((Earth_SLG0[i].dB + Earth_SLG0[i].dC * dt)));
// L1 16x3
for(int i = 0; i < sizeof(Earth_SLG1) / sizeof(VSOP87COEFFICIENT); i++)
dL1 += (Earth_SLG1[i].dA * cos((Earth_SLG1[i].dB + Earth_SLG1[i].dC * dt)));
// L2 10x3
for(int i = 0; i < sizeof(Earth_SLG2) / sizeof(VSOP87COEFFICIENT); i++)
dL2 += (Earth_SLG2[i].dA * cos((Earth_SLG2[i].dB + Earth_SLG2[i].dC * dt)));
// L3 8x3
for(int i = 0; i < sizeof(Earth_SLG3) / sizeof(VSOP87COEFFICIENT); i++)
dL3 += (Earth_SLG3[i].dA * cos((Earth_SLG3[i].dB + Earth_SLG3[i].dC * dt)));
// L4 6x3
for(int i = 0; i < sizeof(Earth_SLG4) / sizeof(VSOP87COEFFICIENT); i++)
dL4 += (Earth_SLG4[i].dA * cos((Earth_SLG4[i].dB + Earth_SLG4[i].dC * dt)));
// L5 1x3
for(int i = 0; i < sizeof(Earth_SLG5) / sizeof(VSOP87COEFFICIENT); i++)
dL5 += (Earth_SLG5[i].dA * cos((Earth_SLG5[i].dB + Earth_SLG5[i].dC * dt)));
// 计算 L = ( L0 + L1 * τ^1 + L2 * τ^2 + L3 * τ^3 + L4 * τ^4 + L5 * τ^5 ) / 10^8 ;（单位弧度）
dL = (dL0 + (dL1 * dt) + (dL2 * (dt * dt)) + (dL3 * (dt * dt * dt)) * (dL4 * (dt * dt * dt * dt)) + (dL5 * (dt * dt * dt * dt * dt))) / 100000000.0;
// 转化为度θ = L + 180 （单位度）;
return (w_MapTo0To360Range(w_MapTo0To360Range(dL / dbUnitRadian) + 180.0));
}
5.3、黄纬：（算法描述请参考相关书籍，为了节省空间和时间，这里请罗列出 C++ 实现代码）
//
// 计算太阳黄纬用参数
const VSOP87COEFFICIENT Earth_SLT0[ 5] =
{
{ 280.0 , 3.199 , 84334.662 } ,
{ 102.0 , 5.422 , 5507.553 } ,
{ 80.0 , 3.880 , 5223.690 } ,
{ 44.0 , 3.700 , 2352.870 } ,
{ 32.0 , 4.000 , 1577.340 }
};
const VSOP87COEFFICIENT Earth_SLT1[ 2] =
{
{ 9.0 , 3.90 , 5507.55 } ,
{ 6.0 , 1.73 , 5223.69 }
};
const VSOP87COEFFICIENT Earth_SLT2[ 4] =
{
{ 22378.0 , 3.38509 , 10213.28555 } ,
{ 282.0 , 0.00000 , 00000.00000 } ,
{ 173.0 , 5.25600 , 20426.57100 } ,
{ 27.0 , 3.87000 , 30639.86000 }
};
const VSOP87COEFFICIENT Earth_SLT3[ 4] =
{
{ 647.0 , 4.992 , 10213.286 } ,
{ 20.0 , 3.140 , 00000.000 } ,
{ 6.0 , 0.770 , 20426.570 } ,
{ 3.0 , 5.440 , 30639.860 }
};
const VSOP87COEFFICIENT CWDCalendarEngine::Earth_SLT4[ 1] =
{
{ 14.0 , 0.32 , 10213.29 }
};
// ------------------------------------------------------------------------
//
// 计算太阳在黄道面上的纬度（单位：度）
//
// dbJD[IN]：儒略日（计算该时刻太阳在黄道面上的纬度）
//
// 返回太阳黄纬度数
// ------------------------------------------------------------------------
double w_GetSunLatitude(const double & dbJD)
{
// 计算τ
double dbt = (dbJD - 2451545.0) / 365250.0;
double dbB = 0.0, dbB0 = 0.0, dbB1 = 0.0, dbB2 = 0.0, dbB3 = 0.0, dbB4 = 0.0;
// B0 5x3
for(int i = 0; i < sizeof(Earth_SLT0) / sizeof(VSOP87COEFFICIENT); i++)
dbB0 += (Earth_SLT0[i].dA * cos((Earth_SLT0[i].dB + Earth_SLT0[i].dC * dbt)));
// B1 2x3
for(int i = 0; i < sizeof(Earth_SLT1) / sizeof(VSOP87COEFFICIENT); i++)
dbB1 += (Earth_SLT1[i].dA * cos((Earth_SLT1[i].dB + Earth_SLT1[i].dC * dbt)));
// B2 4x3
for(int i = 0; i < sizeof(Earth_SLT2) / sizeof(VSOP87COEFFICIENT); i++)
dbB2 += (Earth_SLT2[i].dA * cos((Earth_SLT2[i].dB + Earth_SLT2[i].dC * dbt)));

// B3 4x3
for(int i = 0; i < sizeof(Earth_SLT3) / sizeof(VSOP87COEFFICIENT); i++)
dbB3 += (Earth_SLT3[i].dA * cos((Earth_SLT3[i].dB + Earth_SLT3[i].dC * dbt)));

// B4 1x3
for(int i = 0; i < sizeof(Earth_SLT4) / sizeof(VSOP87COEFFICIENT); i++)
dbB4 += (Earth_SLT4[i].dA * cos((Earth_SLT4[i].dB + Earth_SLT4[i].dC * dbt)));

// 计算 B = ( B0 + B1 * τ^1 + B2 * τ^2 + B3 * τ^3 + B4 * τ^4 ) / 10^8 ;（单位弧度）
dbB = (dbB0 + (dbB1 * dbt) + (dbB2 * (dbt * dbt)) + (dbB3 * (dbt * dbt * dbt)) * (dbB4 * (dbt * dbt * dbt * dbt))) / 100000000.0;
// 计算 θ（单位度）;
return -(dbB / dbUnitRadian);
}


5.4、第一次校正黄经：（算法描述请参考相关书籍，为了节省空间和时间，这里请罗列出 C++ 实现代码）
// ------------------------------------------------------------------------
//
// 修正某时刻太阳在黄道上的经度（黄经）
//
// dbSrcLongitude[IN]：黄经
// dbSrcLatitude[IN]：黄纬
// dbJD[IN]：儒略日（修正在该时刻太阳的黄经）
//
// 返回太阳黄经度数
// ------------------------------------------------------------------------
double w_CorrectionCalcSunLongitude(const double &dbSrcLongitude, const double &dbSrcLatitude, const double &dbJD)
{
double dbT = (dbJD - 2451545.0) / 36525.0;
double dbLdash = dbSrcLongitude - 1.397 * dbT - 0.00031 * dbT * dbT;
// 转换为弧度
dbLdash *= dbUnitRadian;
return (-0.09033 + 0.03916 * (cos(dbLdash) + sin(dbLdash)) * tan(dbSrcLatitude * dbUnitRadian)) / 3600.0;
}
5.5、第二次校正黄经（章动）：（算法描述请参考相关书籍，为了节省空间和时间，这里请罗列出 C++ 实现代码）
//
// 二次修正黄经黄纬所需的天体章动系数
const NUTATIONCOEFFICIENT Nutation_Gene[63] =
{
{ 0, 0, 0, 0, 1, -171996, -174.2, 92025, 8.9 },
{ -2, 0, 0, 2, 2, -13187, -1.6, 5736, -3.1 },
{ 0, 0, 0, 2, 2, -2274, -0.2, 977, -0.5 },
{ 0, 0, 0, 0, 2, 2062, 0.2, -895, 0.5 },
{ 0, 1, 0, 0, 0, 1426, -3.4, 54, -0.1 },
{ 0, 0, 1, 0, 0, 712, 0.1, -7, 0 },
{ -2, 1, 0, 2, 2, -517, 1.2, 224, -0.6 },
{ 0, 0, 0, 2, 1, -386, -0.4, 200, 0 },
{ 0, 0, 1, 2, 2, -301, 0, 129, -0.1 },
{ -2, -1, 0, 2, 2, 217, -0.5, -95, 0.3 },
{ -2, 0, 1, 0, 0, -158, 0, 0, 0 },
{ -2, 0, 0, 2, 1, 129, 0.1, -70, 0 },
{ 0, 0, -1, 2, 2, 123, 0, -53, 0 },
{ 2, 0, 0, 0, 0, 63, 0, 0, 0 },
{ 0, 0, 1, 0, 1, 63, 0.1, -33, 0 },
{ 2, 0, -1, 2, 2, -59, 0, 26, 0 },
{ 0, 0, -1, 0, 1, -58, -0.1, 32, 0 },
{ 0, 0, 1, 2, 1, -51, 0, 27, 0 },
{ -2, 0, 2, 0, 0, 48, 0, 0, 0 },
{ 0, 0, -2, 2, 1, 46, 0, -24, 0 },
{ 2, 0, 0, 2, 2, -38, 0, 16, 0 },
{ 0, 0, 2, 2, 2, -31, 0, 13, 0 },
{ 0, 0, 2, 0, 0, 29, 0, 0, 0 },
{ -2, 0, 1, 2, 2, 29, 0, -12, 0 },
{ 0, 0, 0, 2, 0, 26, 0, 0, 0 },
{ -2, 0, 0, 2, 0, -22, 0, 0, 0 },
{ 0, 0, -1, 2, 1, 21, 0, -10, 0 },
{ 0, 2, 0, 0, 0, 17, -0.1, 0, 0 },
{ 2, 0, -1, 0, 1, 16, 0, -8, 0 },
{ -2, 2, 0, 2, 2, -16, 0.1, 7, 0 },
{ 0, 1, 0, 0, 1, -15, 0, 9, 0 },
{ -2, 0, 1, 0, 1, -13, 0, 7, 0 },
{ 0, -1, 0, 0, 1, -12, 0, 6, 0 },
{ 0, 0, 2,-2, 0, 11, 0, 0, 0 },
{ 2, 0, -1, 2, 1, -10, 0, 5, 0 },
{ 2, 0, 1, 2, 2, -8, 0, 3, 0 },
{ 0, 1, 0, 2, 2, 7, 0, -3, 0 },
{ -2, 1, 1, 0, 0, -7, 0, 0, 0 },
{ 0, -1, 0, 2, 2, -7, 0, 3, 0 },
{ 2, 0, 0, 2, 1, -7, 0, 3, 0 },
{ 2, 0, 1, 0, 0, 6, 0, 0, 0 },
{ -2, 0, 2, 2, 2, 6, 0, -3, 0 },
{ -2, 0, 1, 2, 1, 6, 0, -3, 0 },
{ 2, 0, -2, 0, 1, -6, 0, 3, 0 },
{ 2, 0, 0, 0, 1, -6, 0, 3, 0 },
{ 0, -1, 1, 0, 0, 5, 0, 0, 0 },
{ -2, -1, 0, 2, 1, -5, 0, 3, 0 },
{ -2, 0, 0, 0, 1, -5, 0, 3, 0 },
{ 0, 0, 2, 2, 1, -5, 0, 3, 0 },
{ -2, 0, 2, 0, 1, 4, 0, 0, 0 },
{ -2, 1, 0, 2, 1, 4, 0, 0, 0 },
{ 0, 0, 1,-2, 0, 4, 0, 0, 0 },
{ -1, 0, 1, 0, 0, -4, 0, 0, 0 },
{ -2, 1, 0, 0, 0, -4, 0, 0, 0 },
{ 1, 0, 0, 0, 0, -4, 0, 0, 0 },
{ 0, 0, 1, 2, 0, 3, 0, 0, 0 },
{ 0, 0, -2, 2, 2, -3, 0, 0, 0 },
{ -1, -1, 1, 0, 0, -3, 0, 0, 0 },
{ 0, 1, 1, 0, 0, -3, 0, 0, 0 },
{ 0, -1, 1, 2, 2, -3, 0, 0, 0 },
{ 2, -1, -1, 2, 2, -3, 0, 0, 0 },
{ 0, 0, 3, 2, 2, -3, 0, 0, 0 },
{ 2, -1, 0, 2, 2, -3, 0, 0, 0 }
};
// ------------------------------------------------------------------------
//
// 计算天体章动
// 二次修正某时刻太阳在黄道上的纬度（单位：度）
// 使用天体章动系数修正，消除扰动影响
//
// dbJD[IN]：儒略日（计算该时刻天体章动补偿量）
//
// 返回天体扰动干扰量
// ------------------------------------------------------------------------
double w_GetNutationJamScalar(const double &dbJD)
{
double dbT = (dbJD - 2451545.0) / 36525.0;
double dbTsquared = dbT * dbT;
double dbTcubed = dbTsquared * dbT;
double dbD = 297.85036 + 445267.111480 * dbT - 0.0019142 * dbTsquared + dbTcubed / 189474.0;
dbD = w_MapTo0To360Range(dbD);
double dbM = 357.52772 + 35999.050340 * dbT - 0.0001603 * dbTsquared - dbTcubed / 300000.0;
dbM = w_MapTo0To360Range(dbM);
double dbMprime = 134.96298 + 477198.867398 * dbT + 0.0086972 * dbTsquared + dbTcubed / 56250.0;
dbMprime = w_MapTo0To360Range(dbMprime);
double dbF = 93.27191 + 483202.017538 * dbT - 0.0036825 * dbTsquared + dbTcubed / 327270.0;
dbF = w_MapTo0To360Range(dbF);
double dbOmega = 125.04452 - 1934.136261 * dbT + 0.0020708 * dbTsquared + dbTcubed / 450000.0;
dbOmega = w_MapTo0To360Range(dbOmega);
double dbResulte = 0.0 ;
for(int i = 0; i < sizeof(Nutation_Gene) / sizeof(NUTATIONCOEFFICIENT); i++)
{
double dbRadargument = (Nutation_Gene[i].nD * dbD + Nutation_Gene[i].nM * dbM + Nutation_Gene[i].nMprime * dbMprime + Nutation_Gene[i].nF * dbF + Nutation_Gene[i].nOmega * dbOmega) * dbUnitRadian;
dbResulte += (Nutation_Gene[i].nSincoeff1 + Nutation_Gene[i].dSincoeff2 * dbT ) * sin(dbRadargument) * 0.0001;
}
return dbResulte;
}

5.6、第三次校正黄经（太阳半径向量）：（算法描述请参考相关书籍，为了节省空间和时间，这里请罗列出 C++ 实现代码）
//
// 计算太阳向量半径用参数
const VSOP87COEFFICIENT Earth_SRV0[40] =
{
{ 100013989 , 0 , 0 },
{ 1670700 , 3.0984635 , 6283.0758500 },
{ 13956 , 3.05525 , 12566.15170 },
{ 3084 , 5.1985 , 77713.7715 },
{ 1628 , 1.1739 , 5753.3849 },
{ 1576 , 2.8469 , 7860.4194 },
{ 925 , 5.453 , 11506.770 },
{ 542 , 4.564 , 3930.210 },
{ 472 , 3.661 , 5884.927 },
{ 346 , 0.964 , 5507.553 },
{ 329 , 5.900 , 5223.694 },
{ 307 , 0.299 , 5573.143 },
{ 243 , 4.273 , 11790.629 },
{ 212 , 5.847 , 1577.344 },
{ 186 , 5.022 , 10977.079 },
{ 175 , 3.012 , 18849.228 },
{ 110 , 5.055 , 5486.778 },
{ 98 , 0.89 , 6069.78 },
{ 86 , 5.69 , 15720.84 },
{ 86 , 1.27 , 161000.69 },
{ 65 , 0.27 , 17260.15 },
{ 63 , 0.92 , 529.69 },
{ 57 , 2.01 , 83996.85 },
{ 56 , 5.24 , 71430.70 },
{ 49 , 3.25 , 2544.31 },
{ 47 , 2.58 , 775.52 },
{ 45 , 5.54 , 9437.76 },
{ 43 , 6.01 , 6275.96 },
{ 39 , 5.36 , 4694.00 },
{ 38 , 2.39 , 8827.39 },
{ 37 , 0.83 , 19651.05 },
{ 37 , 4.90 , 12139.55 },
{ 36 , 1.67 , 12036.46 },
{ 35 , 1.84 , 2942.46 },
{ 33 , 0.24 , 7084.90 },
{ 32 , 0.18 , 5088.63 },
{ 32 , 1.78 , 398.15 },
{ 28 , 1.21 , 6286.60 },
{ 28 , 1.90 , 6279.55 },
{ 26 , 4.59 , 10447.39 }
};
const VSOP87COEFFICIENT Earth_SRV1[10] =
{
{ 103019 , 1.107490 , 6283.075850 },
{ 1721 , 1.0644 , 12566.1517 },
{ 702 , 3.142 , 0 },
{ 32 , 1.02 , 18849.23 },
{ 31 , 2.84 , 5507.55 },
{ 25 , 1.32 , 5223.69 },
{ 18 , 1.42 , 1577.34 },
{ 10 , 5.91 , 10977.08 },
{ 9 , 1.42 , 6275.96 },
{ 9 , 0.27 , 5486.78 }
};
const VSOP87COEFFICIENT Earth_SRV2[ 6] =
{
{ 4359 , 5.7846 , 6283.0758 },
{ 124 , 5.579 , 12566.152 },
{ 12 , 3.14 , 0 },
{ 9 , 3.63 , 77713.77 },
{ 6 , 1.87 , 5573.14 },
{ 3 , 5.47 , 18849.23 }
};
const VSOP87COEFFICIENT Earth_SRV3[ 2] =
{
{ 145 , 4.273 , 6283.076 },
{ 7 , 3.92 , 12566.15 }
};
const VSOP87COEFFICIENT Earth_SRV4[ 1] =
{
{ 4 , 2.56 , 6283.08 }
};
// ------------------------------------------------------------------------
//
// 计算某时刻太阳半径向量
// 三次修正某时刻太阳在黄道上的经度（单位：弧度）
//
// dbJD[IN]：儒略日（计算该时刻太阳半径向量）
//
// 返回太阳半径向量
// ------------------------------------------------------------------------
double w_GetSunRadiusVector(const double &dbJD)
{
// 计算τ
double dbt = (dbJD - 2451545.0) / 365250.0;
double dbR = 0.0, dbR0 = 0.0, dbR1 = 0.0, dbR2 = 0.0, dbR3 = 0.0, dbR4 = 0.0;
// R0 40x3
for(int i = 0; i < sizeof(Earth_SRV0) / sizeof(VSOP87COEFFICIENT); i++)
dbR0 += (Earth_SRV0[i].dA * cos((Earth_SRV0[i].dB + Earth_SRV0[i].dC * dbt)));
// R1 10x3
for(int i = 0; i < sizeof(Earth_SRV1) / sizeof(VSOP87COEFFICIENT); i++)
dbR1 += (Earth_SRV1[i].dA * cos((Earth_SRV1[i].dB + Earth_SRV1[i].dC * dbt)));
// R2 6x3
for(int i = 0; i < sizeof(Earth_SRV2) / sizeof(VSOP87COEFFICIENT); i++)
dbR2 += (Earth_SRV2[i].dA * cos((Earth_SRV2[i].dB + Earth_SRV2[i].dC * dbt)));
// R3 2x3
for(int i = 0; i < sizeof(Earth_SRV3) / sizeof(VSOP87COEFFICIENT); i++)
dbR3 += (Earth_SRV3[i].dA * cos((Earth_SRV3[i].dB + Earth_SRV3[i].dC * dbt)));
// R4 1x3
for(int i = 0; i < sizeof(Earth_SRV4) / sizeof(VSOP87COEFFICIENT); i++)
dbR4 += (Earth_SRV4[i].dA * cos((Earth_SRV4[i].dB + Earth_SRV4[i].dC * dbt)));
// 计算 R = ( R0 + R1 * τ^1 + R2 * τ^2 + R3 * τ^3 + R4 * τ^4 ) / 10^8 ;（单位弧度）
return ((dbR0 + (dbR1 * dbt) + (dbR2 * (dbt * dbt)) + (dbR3 * (dbt * dbt * dbt)) * (dbR4 * (dbt * dbt * dbt * dbt))) / 100000000.0);
}
六、计算指定时刻太阳黄道经度度数（算法描述请参考相关书籍，为了节省空间和时间，这里请罗列出 C++ 实现代码）
利用以上介绍算法，如下：
// ------------------------------------------------------------------------
//
// 计算某时刻太阳黄经黄纬
//
// dbJD[IN]：儒略日（计算该时刻太阳在黄道面上的经度和纬度）
// dbLongitude[OUT]：黄经
// dbLatitude[OUT]：黄纬
// ------------------------------------------------------------------------
void w_CalcEclipticLongLat(const double & dbJD, double &dbLongitude, double &dbLatitude)
{
// 计算太阳黄经
dbLongitude = w_GetSunLongitude(dbJD);
// 计算太阳黄纬
dbLatitude = w_GetSunLatitude(dbJD);
// 一次校正经度
dbLongitude += w_CorrectionCalcSunLongitude(dbLongitude, dbLatitude, dbJD);
// 二次校正天体章动
dbLongitude += w_GetNutationJamScalar(dbJD) / 3600.0;
// 三次校正太阳半径向量
dbLongitude -= (20.4898 / w_GetSunRadiusVector(dbJD)) /3600.0;
// 校正太阳黄纬
dbLatitude += w_CorrectionCalcSunLatitude(dbLongitude, dbJD);
}
七、使用二分法计算指定节气的时间
// ------------------------------------------------------------------------
//
// 计算节气
//
// year[IN]：公历年份（计算该年份，指定节气的时间）
// ST_SolarTerms[IN]：节气指定的节气
//
// 返回指定节气的儒略日时间
// ------------------------------------------------------------------------
double w_CalcSolarTerms(const int &year, const SOLARTERMS &ST_SolarTerms)
{
// 节气月份
int SolarTermsMonth = static_cast<int>(ceil(static_cast<double>((ST_SolarTerms + 90.0) / 30.0)));
SolarTermsMonth = SolarTermsMonth > 12 ? SolarTermsMonth - 12 : SolarTermsMonth;
// 节令的发生日期基本都在每月 4 - 9 号间
int LowerLimitSolarTermsDay = ST_SolarTerms % 15 == 0 && ST_SolarTerms % 30 != 0 ? 4 : 16;
// 节气的发生日期基本都在每月 16 - 24 号间
int UpperLimitSolarTermsDay = ST_SolarTerms % 15 == 0 && ST_SolarTerms % 30 != 0 ? 9 : 24;
// 采用二分法逼近计算
double dbLowerLinit = w_GDToJD(year, SolarTermsMonth, LowerLimitSolarTermsDay, 0, 0, 0);
double dbUpperLinit = w_GDToJD(year, SolarTermsMonth, UpperLimitSolarTermsDay, 23, 59, 59);
// 二分法分界点日期
double dbDichotomyDivisionDayJD = 0;
// 太阳黄经角度
double dbLongitude = 0;
// 对比二分法精度是否达到要求
for(; fabs(dbLongitude - static_cast<double>(ST_SolarTerms)) >= 0.00001;)
{
dbDichotomyDivisionDayJD = ((dbUpperLinit - dbLowerLinit) / 2.0) + dbLowerLinit;
// 计算太阳黄经
dbLongitude = w_GetSunLongitude(dbDichotomyDivisionDayJD);
// 一次校正经度
dbLongitude += w_CorrectionCalcSunLongitude(dbLongitude, w_GetSunLatitude(dbDichotomyDivisionDayJD), dbDichotomyDivisionDayJD);
// 二次校正天体章动
dbLongitude += w_GetNutationJamScalar(dbDichotomyDivisionDayJD) / 3600.0;
// 三次校正太阳半径向量
dbLongitude -= (20.4898 / w_GetSunRadiusVector(dbDichotomyDivisionDayJD)) /3600.0;
// 由于春分这天黄经为 0 度，比较特殊，因此专门判断（如不加以特殊对待则会导致计算范围覆盖整个 360 度角）
dbLongitude = ((ST_SolarTerms == ST_VERNAL_EQUINOX) && (dbLongitude > 345.0)) ? -dbLongitude : dbLongitude;
// 调整二分法上下限
(dbLongitude > static_cast<double>(ST_SolarTerms)) ? dbUpperLinit = dbDichotomyDivisionDayJD : dbLowerLinit = dbDichotomyDivisionDayJD;
}
return dbDichotomyDivisionDayJD;
}
八、本地时间（LST）和格林尼治时间（UTC）的互换（仅限于 VC 代码）
// ------------------------------------------------------------------------
//
// 格林威治时间转本地时间（以格里历表示本地时间）
//
// ------------------------------------------------------------------------
void w_UTCToLST(int &year, int &month, int &day, int &hour, int &minute, int &second)
{
_tzset () ;
// 计算本地时间和标准时间的时差（单位：秒）
int nDifference_hour = static_cast<int>(_timezone / 3600);
int nDifference_minute = static_cast<int>((_timezone - nDifference_hour * 3600) / 60);
int nDifference_second = static_cast<int>((_timezone - nDifference_hour * 3600) - nDifference_minute * 60);
// 格林威治时间 + 时差 = 本地时间
// 秒
second = second - nDifference_second;
if(second >= 60 || second < 0)
{
minute = second > 0 ? minute + 1 : minute - 1 ;
second = abs(abs(second) - 60);
}
// 分
minute = minute - nDifference_minute;
if(minute >= 60 || minute < 0)
{
hour = minute > 0 ? hour + 1 : hour - 1;
minute = abs(abs(minute) - 60);
}
// 时
hour = hour - nDifference_hour;
if(hour >= 24 || hour < 0)
{
day = (hour >= 24 || hour == 0) ? day + 1 : day - 1;
hour = abs(abs(hour) - 24);
}
// 日
int nDaysOfMonth = w_GetDaysOfMonth(year, month);
if(day > nDaysOfMonth || day <= 0)
{
if(day > nDaysOfMonth)
month++;
if(day < nDaysOfMonth || day <= 0)
month--;
day = abs(abs(day) - nDaysOfMonth);
}
// 月
if(month > 12 || month <= 0)
{
year = month > 0 ? year + 1 : year - 1;
month = month > 0 ? abs(month - 12) : abs(12 + month);
}
}
// ------------------------------------------------------------------------
//
// 本地时间转格林威治时间（以格里历表示）
//
// ------------------------------------------------------------------------
void w_LSTToUTC(int &year, int &month, int &day, int &hour, int &minute, int &second)
{
_tzset () ;
// 计算本地时间和标准时间的时差（单位：秒）
int nDifference_hour = static_cast<int>(_timezone / 3600);
int nDifference_minute = static_cast<int>((_timezone - nDifference_hour * 3600) / 60);
int nDifference_second = static_cast<int>((_timezone - nDifference_hour * 3600) - nDifference_minute * 60);
// 本地时间 - 时差 = 格林威治时间
// 秒
second = second + nDifference_second;
if(second >= 60 || second < 0)
{
minute = second > 0 ? minute + 1 : minute - 1 ;
second = abs(abs(second) - 60);
}
// 分
minute = minute + nDifference_minute;
if(minute >= 60 || minute < 0)
{
hour = minute > 0 ? hour + 1 : hour - 1;
minute = abs(abs(minute) - 60);
}
// 时
hour = hour + nDifference_hour;
if(hour >= 24 || hour < 0)
{
day = (hour >= 24 || hour == 0) ? day + 1 : day - 1;
hour = abs(abs(hour) - 24);
}
// 日
int nDaysOfMonth = w_GetDaysOfMonth(year, month);
if(day > nDaysOfMonth || day <= 0)
{
if(day > nDaysOfMonth)
month++;
if(day < nDaysOfMonth || day <= 0)
month--;
day = abs(abs(day) - nDaysOfMonth);
}
// 月
if(month > 12 || month <= 0)
{
year = month > 0 ? year + 1 : year - 1;
month = month > 0 ? abs(month - 12) : abs(12 + month);
}
}


九、2008年全年 24 节气发生时间（北京时间）
春分：03月20日 13:49:24
清明：04月04日 17:46:58
谷雨：04月20日 00:52:17
立夏：05月05日 11:04:35
小满：05月21日 00:02:03
芒种：06月05日 15:12:53
夏至：06月21日 08:00:30
小暑：07月07日 01:27:59
大暑：07月22日 18:55:58
立秋：08月07日 11:17:19
处暑：08月23日 02:03:22
白露：09月07日 14:15:15
秋分：09月22日 23:45:36
寒露：10月08日 05:57:46
霜降：10月23日 09:09:49
立冬：11月07日 09:11:45
小雪：11月22日 06:45:32
大雪：12月07日 02:03:28
冬至：12月21日 20:04:54
小寒：01月06日 07:25:55
大寒：01月21日 00:44:38
立春：02月04日 19:01:30
雨水：02月19日 14:50:40
惊蛰：03月05日 12:59:54
本文给出了一种通过计算得到某一年的节气发生日期，相比之下比通过查表的方法的优点不言而喻，经本人比对（与日梭万年历）在范围公元３００－３０００年间误差很小，如果对历法进一步熟悉，你会发现历法的计算非常困难，计算的误差也随时间跨度的增加增长的很快，即便是所有查阅到的最权威的书籍数据，也只能保证在前后２００年间的准确。