万年历计算 之 节气（续）
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
本文给出了一种通过计算得到某一年的节气发生日期，相比之下比通过查表的方法的优点不言而喻，经本人比对（与日梭万年历）在范围公元３００－３０００年间误差很小，如果对历法进一步熟悉，你会发现历法的计算非常困难，计算的误差也随时间跨度的增加增长的很快，即便是所有查阅到的最权威的书籍数据，也只能保证在前后２００年间的准确。
最后欢迎所有朋友提出你的意见。（全文完）
导入论坛 引用链接 收藏 分享给好友 推荐到圈子 管理 举报
TAG:
标题搜索
日历
|
|||||||||
| 日 | 一 | 二 | 三 | 四 | 五 | 六 | |||
| 1 | |||||||||
| 2 | 3 | 4 | 5 | 6 | 7 | 8 | |||
| 9 | 10 | 11 | 12 | 13 | 14 | 15 | |||
| 16 | 17 | 18 | 19 | 20 | 21 | 22 | |||
| 23 | 24 | 25 | 26 | 27 | 28 | 29 | |||
| 30 | |||||||||
我的存档
数据统计
- 访问量: 1155
- 日志数: 3
- 建立时间: 2007-01-29
- 更新时间: 2008-05-07
