- 浏览: 40842 次
- 性别:
- 来自: 北京
最新评论
中国农历的朔望月是农历历法的基础,而朔望月又是严格以日月合朔发生的那一天作为月首,因此日月合朔时间的计算是制定农历历法的关键。本文将介绍ELP-2000/82月球运行理论,以及如何用ELP-2000/82月球运行理论计算日月合朔时间。
要计算日月合朔时间,首先要对日月合朔这一天文现象进行数学定义。朔望月是在地球上观察到的月相周期,平均长度约等于29.53059日,而恒星月(天文月)是月亮绕地球公转一周的时间,长度约27.32166日。月相周期长度比恒星月长大约两天,这是因为在月球绕地球旋转一周的同时,地球还带着它绕太阳旋转了一定的角度的缘故,所以月相周期不仅与月球运行有关,还和太阳运行有关。日月合朔的时候,太阳、月亮和地球三者接近一条直线,月亮未被照亮的一面对着地球,因此地球上看不到月亮,此时又被称为新月。图(1)就是日月合朔天文现象的示意图:
图(1)日月天文现象示意图
月亮绕太阳公转的白道面和地球绕太阳公转的黄道面存在一个最大约5°的夹角,因此大多数情况下,日月合朔时都不是严格在同一条直线上,不过也会发生在同一直线的情况,此时就会发生日食。图(1-b)显示了日月合朔时侧切面上月亮的三种可能的位置情况,当月亮处在位置2时就会发生日食。由图(1)可知,日月合朔的数学定义就是太阳和月亮的地心视黄经差为0的时刻。
要计算日月合朔,需要知道太阳地心视黄经和月亮地心视黄经的计算方法。“日历生成算法”系列的第三篇《用天文方法计算二十四节气》一文已经介绍了如何用VSOP82/87行星理论计算太阳的地心视黄经,本文将继续介绍如何用ELP-2000/82月球理论计算月亮的地心视黄经。ELP-2000/82月球理论是M. Chapront-Touze和J. Chapront在1983年提出的一个月球位置的半解析理论,和其它半解析理论一样,ELP-2000/82理论也包含一套计算方法和相应的迭代周期项。这套理论共包含37862个周期项,其中20560个用于计算月球经度,7684个用于计算月球纬度,9618个用于计算地月距离。但是这些周期项中有很多都是非常小的值,例如一些计算经纬度的项对结果的增益只有0.00001角秒,还有一些地月距离周期项对距离结果的增益只有0.02米,对于精度不高的历法计算,完全可以忽略。
有很多基于ELP-2000/82月球理论的改进或简化理论,《天文算法》一书的第四十五章就介绍了一种改进算法,其周期项参数都是从ELP-2000/82理论的周期项参数转换来的,忽略了小的周期项。使用该方法计算的月球黄经精度只有10”,月亮黄纬精度只有4”,但是只用计算60个周期项,速度很快,本文就采用这种修改过的ELP-2000/82理论计算月亮的地心视黄经。这种计算方法的周期项分三部分,分别用来计算月球黄经,月球黄纬和地月距离,三部分的周期项的内容一样,由四个计算辐角的系数和一个正弦(或余弦)振幅组成。计算月球黄经和地月距离使用正弦表达式求和:A * sin(θ),计算月球黄纬用余弦表达式求和:A * cos(θ),其中辐角θ的计算公式是:
θ = a * D + b * M + c * M’ + d * F (4.1式)
4.1式中的四个辐角系数a、b、c和d由每个迭代周期项给出,日月距角D、太阳平近地角M、月亮平近地角M’以及月球生交点平角距F则分别有4.2式-4.5式进行计算:
D = 297.8502042 + 445267.1115168 * T - 0.0016300 * T2
+ T3 / 545868 - T4 / 113065000 (4.2式)
M = 357.5291092 + 35999.0502909 * T - 0.0001536 * T2
+ T3 / 24490000 (4.3式)
M' = 134.9634114 + 477198.8676313 * T + 0.0089970 * T2
+ T3 / 69699 - T4 / 14712000 (4.4式)
F = 93.2720993 + 483202.0175273 * T - 0.0034029 * T2
- T3 / 3526000 + T4 / 863310000 (4.5式)
以上各式计算结果的单位是度,其中T是儒略世纪数,T计算由4.6式计算:
T = (JDE - 2451545.0) / 36525.0 (4.6式)
以计算月球黄经的周期项第二项的计算为例,第二项数据如下,辐角系数a = 2,b = 0,c = -1,d = 0,振幅A = 1274027,黄经计算用正弦表达式,则I2的计算如下所示:
I2 = 1274027 * sin(2D – M’) (4.7式)
在套用4.7式计算出60个月球黄经周期项值的时候,需要注意对包含了太阳平近地角M的项进行修正,因为M的值与地球公转轨道的离心率有关,因为离心率是个与时间有关的变量,导致振幅A实际上是个变量,需要根据时间进行修正。月球黄经周期项的修正方法是:如果辐角中包含了M或-M时,需要乘以系数E修正;如果辐角中包含了2M或-2M,则需要乘以系数E的平方进行修正。系数E的计算表达式如下:
E = 1 - 0.002516 * T - 0.0000074 * T2 (4.8式)
其中T值由4.6式计算。上面的计算月球黄经的第二个周期项中M对应的系数是0,因此I2不需要修正,但是第五个周期项中M对应的系数是1,因此I5需要乘以E进行修正。套用4.7式计算出60个月球黄经周期项值I1-I60之和ΣI:
ΣI = I1 + I2 + … + I60 (4.9式)
月球黄纬的周期项和Σb的计算方法与月球黄经周期项和ΣI的计算方法一样,每个月球黄纬周期项也包含振幅A和四个辐角系数a、b、c和d,对于太阳平近地角M的系数b不是0的情况也需要乘以E或E2进行修正,唯一的区别是计算月球黄纬周期项用余弦表达式。地月距离的周期项和Σr也可以按照上面的方法计算,计算地月距离的目的是为了计算月亮光行差,因为地月距离较小,从地球观察月亮产生的光行差也很小,相对于本文算法的精度(月球黄经精度10”,月亮黄纬精度4”)来说,可以忽略光行差修正,因此就不用计算地月距离。
由于金星和木星对月球的摄动影响,需要对计算出的月球黄经周期项和ΣI和月球黄纬周期项和Σb金星摄动修正,修正的方法如下:
ΣI += +3958 * sin( A1 ) + 1962 * sin( L' - F ) + 318 * sin( A2 ) (4.10式)
Σb += -2235 * sin( L' ) + 382 * sin( A3) + 175 * sin( A1 - F ) + 175 * sin( A1 + F )
+ 127 * sin( L' - M') - 115 * sin( L' + M') (4.11式)
其中M’和F分别由4.4式和4.5式计算得到,L’是月球平黄经,计算方法是:
L'=218.3164591 + 481267.88134236 * T - 0.0013268 * T2
+ T3 / 538841 - T4 / 65194000 (4.12式)
A1、A2和A4是摄动角修正量,计算方法如下:
A1 = 119.75 + 131.849 * T (4.13式)
A2 = 53.09 + 479264.290 * T (4.14式)
A3 = 313.45 + 481266.484 * T (4.15式)
完成所有修正后,就可以用4.16式和4.17式最终得到月亮的地心视黄经和地心视黄纬:
λ = L'+ ΣI / 1000000.0 (4.16式)
β = Σb / 1000000.0 (4.17式)
ΣI和Σb最后要除以1000000.0是因为周期项系数中振幅A的单位是0.000001度,因此λ和β的单位是度。下面给出计算月球地心视黄经的代码:
123 double GetMoonEclipticLongitudeEC(double dbJD)
124 {
125 double Lp,D,M,Mp,F,E;
126 double dt = (dbJD - JD2000) / 36525.0; /*儒略世纪数*/
127
128 GetMoonEclipticParameter(dt, &Lp, &D, &M, &Mp, &F, &E);
129
130 /*计算月球地心黄经周期项*/
131 double EI = CalcMoonECLongitudePeriodicTbl(D, M, Mp, F, E);
132
133 /*修正金星,木星以及地球扁率摄动*/
134 EI += CalcMoonLongitudePerturbation(dt, Lp, F);
135
136 /*计算月球地心黄经*/
137 double longitude = Lp + EI / 1000000.0;
138
139 /*计算天体章动干扰*/
140 longitude += CalcEarthLongitudeNutation(dt / 10.0);
141
142 longitude = Mod360Degree(longitude);
143 return longitude;
144 }
函数参数dbJD是力学时儒略日时间,返回以度为单位的月球视黄经。其中GetMoonEclipticParameter()函数分别根据4.2式、4.3式、4.4式、4.5式、4.8式和4.12式计算日月距角D、太阳平近地角M、月亮平近地角M’、月球生交点平角距F、修正系数E和月球平黄经L’,不需多说明,只要根据以上各式直接计算即可。CalcMoonECLongitudePeriodicTbl()函数计算60个月球黄经周期项和,并根据M值系数的情况进行修正,算法实现如下:
42 double CalcMoonECLongitudePeriodicTbl(double D, double M, double Mp, double F, double E)
43 {
44 double EI = 0.0 ;
45
46 for(int i = 0; i < COUNT_ITEM(Moon_longitude); i++)
47 {
48 double sita = Moon_longitude[i].D * D + Moon_longitude[i].M * M + Moon_longitude[i].Mp * Mp + Moon_longitude[i].F * F;
49 sita = DegreeToRadian(sita);
50 EI += (Moon_longitude[i].eiA * sin(sita) * pow(E, fabs(Moon_longitude[i].M)));
51 }
52
53 return EI;
54 }
CalcMoonLongitudePerturbation()函数计算月球黄经摄动修正量,使用了4.13式和4.14式给出的计算方法:
87 double CalcMoonLongitudePerturbation(double dt, double Lp, double F)
88 {
89 double T = dt; /*T是'ca?从'b4?J2000起'c6?算'cb?的'b5?儒'c8?略'c2?世'ca?纪'bc?数'ca?*/
90 double A1 = 119.75 + 131.849 * T;
91 double A2 = 53.09 + 479264.290 * T;
92
93 A1 = Mod360Degree(A1);
94 A2 = Mod360Degree(A2);
95
96 double result = 3958.0 * sin(DegreeToRadian(A1));
97 result += (1962.0 * sin(DegreeToRadian(Lp - F)));
98 result += (318.0 * sin(DegreeToRadian(A2)));
99
100 return result;
101 }
至此,本文已经介绍了使用ELP-2000/82月球理论计算任意时刻月亮地心视黄经的方法,结合“日历生成算法”系列的第三篇《用天文方法计算二十四节气》一文介绍的计算太阳地心视黄经的方法,就可以计算日月合朔的准确时间了。由于ELP-2000/82月球理论也没有根据月球黄经反算时间的方法,因此本文也采用和《用天文方法计算二十四节气》一文中一样的牛顿迭代法计算日月合朔时间。
关于牛顿迭代法可以参考相关的数学资料,“日历生成算法”系列的第三篇《用天文方法计算二十四节气》一文对如何使用牛顿迭代法有简单的介绍,可以参考一下。总的来说,就是要先定义需要求解的方程f(x),根据上文的介绍,我们需要求解的是太阳的地心黄经和月亮的地心黄经差值是0的时候的时间,《用天文方法计算二十四节气》一文已经介绍了求太阳地心黄经的函数GetSunEclipticLongitudeECDegree(),本文也给出了求月亮地心黄经的函数GetMoonEclipticLongitudeECDegree(),因此可以定义方程为:
f(x) = GetMoonEclipticLongitudeECDegree(x) – GetSunEclipticLongitudeECDegree(x) = 0
其中x是儒略日单位的,我们要用牛顿迭代法求方程f(x)=0时的解x,也就是时间值。牛顿迭代法求解的迭代式是:
Xn+1 = Xn – f(Xn)/f’(Xn)
这里也不多解释了。导函数仍然使用近似公式,也不解释了,直接上迭代求解的代码了:
102 double CalculateMoonShuoJD(double tdJD)
103 {
104 double JD0, JD1,stDegree,stDegreep;
105
106 JD1 = tdJD;
107 do
108 {
109 JD0 = JD1;
110 double moonLongitude = GetMoonEclipticLongitudeECDegree(JD0);
111 double sunLongitude = GetSunEclipticLongitudeECDegree(JD0);
112
113 stDegree = moonLongitude - sunLongitude;
114
115
116 stDegreep = (GetMoonEclipticLongitudeECDegree(JD0 + 0.000005) - GetSunEclipticLongitudeECDegree(JD0 + 0.000005) - GetMoonEclipticLongitudeECDegree(JD0 - 0.000005) + GetSunEclipticLongitudeECDegree(JD0 - 0.000005)) / 0.00001;
117 JD1 = JD0 - stDegree / stDegreep;
118 }while((fabs(JD1 - JD0) > 0.00000001));
119
120 return JD1;
121 }
至本文结束,我们已经能够使用半解析算法计算太阳的黄经和月亮的黄经,并且能够通过牛顿迭代法或者24节气的准确时间和日月合朔的准确时间,在这基础上就可以进行中国农历的推算了,“日历生成算法”系列的下一篇将介绍中国农历的历法规则和推算方法。
再次说明一下,以上算法中讨论的时间都是力学时时间(TD),与国际协调时(UTC)以及各个时区的本地时间都有不同,以上计算出来的时间都需要调整成本地时间,比如中国的中原地区就是东八区标准时(UTC + 8)。应用本文的算法计算出2012年前后15个日月合朔时间如下(已经转换为东八区标准时):
2011-11-25, 14:09:41.25
2011-12-25, 02:06:27.25
2012-01-23, 15:39:24.16
2012-02-22, 06:34:40.84
2012-03-22, 22:37:08.91
2012-04-21, 15:18:22.12
2012-05-21, 07:46:59.97
2012-06-19, 23:02:06.39
2012-07-19, 12:24:02.83
2012-08-17, 23:54:28.03
2012-09-16, 10:10:36.99
2012-10-15, 20:02:30.98
2012-11-14, 06:08:05.90
2012-12-13, 16:41:37.60
2013-01-12, 03:43:31.34
要计算日月合朔时间,首先要对日月合朔这一天文现象进行数学定义。朔望月是在地球上观察到的月相周期,平均长度约等于29.53059日,而恒星月(天文月)是月亮绕地球公转一周的时间,长度约27.32166日。月相周期长度比恒星月长大约两天,这是因为在月球绕地球旋转一周的同时,地球还带着它绕太阳旋转了一定的角度的缘故,所以月相周期不仅与月球运行有关,还和太阳运行有关。日月合朔的时候,太阳、月亮和地球三者接近一条直线,月亮未被照亮的一面对着地球,因此地球上看不到月亮,此时又被称为新月。图(1)就是日月合朔天文现象的示意图:
图(1)日月天文现象示意图
月亮绕太阳公转的白道面和地球绕太阳公转的黄道面存在一个最大约5°的夹角,因此大多数情况下,日月合朔时都不是严格在同一条直线上,不过也会发生在同一直线的情况,此时就会发生日食。图(1-b)显示了日月合朔时侧切面上月亮的三种可能的位置情况,当月亮处在位置2时就会发生日食。由图(1)可知,日月合朔的数学定义就是太阳和月亮的地心视黄经差为0的时刻。
要计算日月合朔,需要知道太阳地心视黄经和月亮地心视黄经的计算方法。“日历生成算法”系列的第三篇《用天文方法计算二十四节气》一文已经介绍了如何用VSOP82/87行星理论计算太阳的地心视黄经,本文将继续介绍如何用ELP-2000/82月球理论计算月亮的地心视黄经。ELP-2000/82月球理论是M. Chapront-Touze和J. Chapront在1983年提出的一个月球位置的半解析理论,和其它半解析理论一样,ELP-2000/82理论也包含一套计算方法和相应的迭代周期项。这套理论共包含37862个周期项,其中20560个用于计算月球经度,7684个用于计算月球纬度,9618个用于计算地月距离。但是这些周期项中有很多都是非常小的值,例如一些计算经纬度的项对结果的增益只有0.00001角秒,还有一些地月距离周期项对距离结果的增益只有0.02米,对于精度不高的历法计算,完全可以忽略。
有很多基于ELP-2000/82月球理论的改进或简化理论,《天文算法》一书的第四十五章就介绍了一种改进算法,其周期项参数都是从ELP-2000/82理论的周期项参数转换来的,忽略了小的周期项。使用该方法计算的月球黄经精度只有10”,月亮黄纬精度只有4”,但是只用计算60个周期项,速度很快,本文就采用这种修改过的ELP-2000/82理论计算月亮的地心视黄经。这种计算方法的周期项分三部分,分别用来计算月球黄经,月球黄纬和地月距离,三部分的周期项的内容一样,由四个计算辐角的系数和一个正弦(或余弦)振幅组成。计算月球黄经和地月距离使用正弦表达式求和:A * sin(θ),计算月球黄纬用余弦表达式求和:A * cos(θ),其中辐角θ的计算公式是:
θ = a * D + b * M + c * M’ + d * F (4.1式)
4.1式中的四个辐角系数a、b、c和d由每个迭代周期项给出,日月距角D、太阳平近地角M、月亮平近地角M’以及月球生交点平角距F则分别有4.2式-4.5式进行计算:
D = 297.8502042 + 445267.1115168 * T - 0.0016300 * T2
+ T3 / 545868 - T4 / 113065000 (4.2式)
M = 357.5291092 + 35999.0502909 * T - 0.0001536 * T2
+ T3 / 24490000 (4.3式)
M' = 134.9634114 + 477198.8676313 * T + 0.0089970 * T2
+ T3 / 69699 - T4 / 14712000 (4.4式)
F = 93.2720993 + 483202.0175273 * T - 0.0034029 * T2
- T3 / 3526000 + T4 / 863310000 (4.5式)
以上各式计算结果的单位是度,其中T是儒略世纪数,T计算由4.6式计算:
T = (JDE - 2451545.0) / 36525.0 (4.6式)
以计算月球黄经的周期项第二项的计算为例,第二项数据如下,辐角系数a = 2,b = 0,c = -1,d = 0,振幅A = 1274027,黄经计算用正弦表达式,则I2的计算如下所示:
I2 = 1274027 * sin(2D – M’) (4.7式)
在套用4.7式计算出60个月球黄经周期项值的时候,需要注意对包含了太阳平近地角M的项进行修正,因为M的值与地球公转轨道的离心率有关,因为离心率是个与时间有关的变量,导致振幅A实际上是个变量,需要根据时间进行修正。月球黄经周期项的修正方法是:如果辐角中包含了M或-M时,需要乘以系数E修正;如果辐角中包含了2M或-2M,则需要乘以系数E的平方进行修正。系数E的计算表达式如下:
E = 1 - 0.002516 * T - 0.0000074 * T2 (4.8式)
其中T值由4.6式计算。上面的计算月球黄经的第二个周期项中M对应的系数是0,因此I2不需要修正,但是第五个周期项中M对应的系数是1,因此I5需要乘以E进行修正。套用4.7式计算出60个月球黄经周期项值I1-I60之和ΣI:
ΣI = I1 + I2 + … + I60 (4.9式)
月球黄纬的周期项和Σb的计算方法与月球黄经周期项和ΣI的计算方法一样,每个月球黄纬周期项也包含振幅A和四个辐角系数a、b、c和d,对于太阳平近地角M的系数b不是0的情况也需要乘以E或E2进行修正,唯一的区别是计算月球黄纬周期项用余弦表达式。地月距离的周期项和Σr也可以按照上面的方法计算,计算地月距离的目的是为了计算月亮光行差,因为地月距离较小,从地球观察月亮产生的光行差也很小,相对于本文算法的精度(月球黄经精度10”,月亮黄纬精度4”)来说,可以忽略光行差修正,因此就不用计算地月距离。
由于金星和木星对月球的摄动影响,需要对计算出的月球黄经周期项和ΣI和月球黄纬周期项和Σb金星摄动修正,修正的方法如下:
ΣI += +3958 * sin( A1 ) + 1962 * sin( L' - F ) + 318 * sin( A2 ) (4.10式)
Σb += -2235 * sin( L' ) + 382 * sin( A3) + 175 * sin( A1 - F ) + 175 * sin( A1 + F )
+ 127 * sin( L' - M') - 115 * sin( L' + M') (4.11式)
其中M’和F分别由4.4式和4.5式计算得到,L’是月球平黄经,计算方法是:
L'=218.3164591 + 481267.88134236 * T - 0.0013268 * T2
+ T3 / 538841 - T4 / 65194000 (4.12式)
A1、A2和A4是摄动角修正量,计算方法如下:
A1 = 119.75 + 131.849 * T (4.13式)
A2 = 53.09 + 479264.290 * T (4.14式)
A3 = 313.45 + 481266.484 * T (4.15式)
完成所有修正后,就可以用4.16式和4.17式最终得到月亮的地心视黄经和地心视黄纬:
λ = L'+ ΣI / 1000000.0 (4.16式)
β = Σb / 1000000.0 (4.17式)
ΣI和Σb最后要除以1000000.0是因为周期项系数中振幅A的单位是0.000001度,因此λ和β的单位是度。下面给出计算月球地心视黄经的代码:
123 double GetMoonEclipticLongitudeEC(double dbJD)
124 {
125 double Lp,D,M,Mp,F,E;
126 double dt = (dbJD - JD2000) / 36525.0; /*儒略世纪数*/
127
128 GetMoonEclipticParameter(dt, &Lp, &D, &M, &Mp, &F, &E);
129
130 /*计算月球地心黄经周期项*/
131 double EI = CalcMoonECLongitudePeriodicTbl(D, M, Mp, F, E);
132
133 /*修正金星,木星以及地球扁率摄动*/
134 EI += CalcMoonLongitudePerturbation(dt, Lp, F);
135
136 /*计算月球地心黄经*/
137 double longitude = Lp + EI / 1000000.0;
138
139 /*计算天体章动干扰*/
140 longitude += CalcEarthLongitudeNutation(dt / 10.0);
141
142 longitude = Mod360Degree(longitude);
143 return longitude;
144 }
函数参数dbJD是力学时儒略日时间,返回以度为单位的月球视黄经。其中GetMoonEclipticParameter()函数分别根据4.2式、4.3式、4.4式、4.5式、4.8式和4.12式计算日月距角D、太阳平近地角M、月亮平近地角M’、月球生交点平角距F、修正系数E和月球平黄经L’,不需多说明,只要根据以上各式直接计算即可。CalcMoonECLongitudePeriodicTbl()函数计算60个月球黄经周期项和,并根据M值系数的情况进行修正,算法实现如下:
42 double CalcMoonECLongitudePeriodicTbl(double D, double M, double Mp, double F, double E)
43 {
44 double EI = 0.0 ;
45
46 for(int i = 0; i < COUNT_ITEM(Moon_longitude); i++)
47 {
48 double sita = Moon_longitude[i].D * D + Moon_longitude[i].M * M + Moon_longitude[i].Mp * Mp + Moon_longitude[i].F * F;
49 sita = DegreeToRadian(sita);
50 EI += (Moon_longitude[i].eiA * sin(sita) * pow(E, fabs(Moon_longitude[i].M)));
51 }
52
53 return EI;
54 }
CalcMoonLongitudePerturbation()函数计算月球黄经摄动修正量,使用了4.13式和4.14式给出的计算方法:
87 double CalcMoonLongitudePerturbation(double dt, double Lp, double F)
88 {
89 double T = dt; /*T是'ca?从'b4?J2000起'c6?算'cb?的'b5?儒'c8?略'c2?世'ca?纪'bc?数'ca?*/
90 double A1 = 119.75 + 131.849 * T;
91 double A2 = 53.09 + 479264.290 * T;
92
93 A1 = Mod360Degree(A1);
94 A2 = Mod360Degree(A2);
95
96 double result = 3958.0 * sin(DegreeToRadian(A1));
97 result += (1962.0 * sin(DegreeToRadian(Lp - F)));
98 result += (318.0 * sin(DegreeToRadian(A2)));
99
100 return result;
101 }
至此,本文已经介绍了使用ELP-2000/82月球理论计算任意时刻月亮地心视黄经的方法,结合“日历生成算法”系列的第三篇《用天文方法计算二十四节气》一文介绍的计算太阳地心视黄经的方法,就可以计算日月合朔的准确时间了。由于ELP-2000/82月球理论也没有根据月球黄经反算时间的方法,因此本文也采用和《用天文方法计算二十四节气》一文中一样的牛顿迭代法计算日月合朔时间。
关于牛顿迭代法可以参考相关的数学资料,“日历生成算法”系列的第三篇《用天文方法计算二十四节气》一文对如何使用牛顿迭代法有简单的介绍,可以参考一下。总的来说,就是要先定义需要求解的方程f(x),根据上文的介绍,我们需要求解的是太阳的地心黄经和月亮的地心黄经差值是0的时候的时间,《用天文方法计算二十四节气》一文已经介绍了求太阳地心黄经的函数GetSunEclipticLongitudeECDegree(),本文也给出了求月亮地心黄经的函数GetMoonEclipticLongitudeECDegree(),因此可以定义方程为:
f(x) = GetMoonEclipticLongitudeECDegree(x) – GetSunEclipticLongitudeECDegree(x) = 0
其中x是儒略日单位的,我们要用牛顿迭代法求方程f(x)=0时的解x,也就是时间值。牛顿迭代法求解的迭代式是:
Xn+1 = Xn – f(Xn)/f’(Xn)
这里也不多解释了。导函数仍然使用近似公式,也不解释了,直接上迭代求解的代码了:
102 double CalculateMoonShuoJD(double tdJD)
103 {
104 double JD0, JD1,stDegree,stDegreep;
105
106 JD1 = tdJD;
107 do
108 {
109 JD0 = JD1;
110 double moonLongitude = GetMoonEclipticLongitudeECDegree(JD0);
111 double sunLongitude = GetSunEclipticLongitudeECDegree(JD0);
112
113 stDegree = moonLongitude - sunLongitude;
114
115
116 stDegreep = (GetMoonEclipticLongitudeECDegree(JD0 + 0.000005) - GetSunEclipticLongitudeECDegree(JD0 + 0.000005) - GetMoonEclipticLongitudeECDegree(JD0 - 0.000005) + GetSunEclipticLongitudeECDegree(JD0 - 0.000005)) / 0.00001;
117 JD1 = JD0 - stDegree / stDegreep;
118 }while((fabs(JD1 - JD0) > 0.00000001));
119
120 return JD1;
121 }
至本文结束,我们已经能够使用半解析算法计算太阳的黄经和月亮的黄经,并且能够通过牛顿迭代法或者24节气的准确时间和日月合朔的准确时间,在这基础上就可以进行中国农历的推算了,“日历生成算法”系列的下一篇将介绍中国农历的历法规则和推算方法。
再次说明一下,以上算法中讨论的时间都是力学时时间(TD),与国际协调时(UTC)以及各个时区的本地时间都有不同,以上计算出来的时间都需要调整成本地时间,比如中国的中原地区就是东八区标准时(UTC + 8)。应用本文的算法计算出2012年前后15个日月合朔时间如下(已经转换为东八区标准时):
2011-11-25, 14:09:41.25
2011-12-25, 02:06:27.25
2012-01-23, 15:39:24.16
2012-02-22, 06:34:40.84
2012-03-22, 22:37:08.91
2012-04-21, 15:18:22.12
2012-05-21, 07:46:59.97
2012-06-19, 23:02:06.39
2012-07-19, 12:24:02.83
2012-08-17, 23:54:28.03
2012-09-16, 10:10:36.99
2012-10-15, 20:02:30.98
2012-11-14, 06:08:05.90
2012-12-13, 16:41:37.60
2013-01-12, 03:43:31.34
发表评论
-
用天文方法计算二十四节气
2012-12-15 09:08 2472二十四节气在中国古代历法中扮演着非常重要的角色 ... -
Positional Astronomy
2012-12-10 05:49 886http://www.jgiesen.de/elevaz/ba ... -
Ephemeris trail
2012-07-13 14:46 855http://www.astrosurf.com/jephem ... -
de406行星历表的结构
2012-07-12 17:14 1697以下 header.405 KSIZE= 2036 NC ... -
DE405/406星历表算法
2012-07-12 17:16 7175#pragma hdrstop #pragma argsuse ... -
农历24节气算法
2012-07-12 17:16 2596许剑伟 [摘要] 古老而 ... -
万年历计算之节气
2012-06-29 14:35 2011一、基本知识 二 ... -
万年历计算之干支
2012-06-29 14:20 18911、基本知识 ... -
行星的位置计算
2010-11-20 17:25 3827许剑伟 莆田十中 [摘 ...
相关推荐
算法系列之十八:用天文方法计算二十四节气(上) 本文将介绍二十四节气的基本知识,以及如何使用 VSOP82/87 行星运行理论计算二十四节气发生的准确时间。二十四节气在中国古代历法中扮演着非常重要的角色,它是...
天文算法万年历 寿星万年历的易语言源码!
算法系列之十八:用天文方法计算二十四节气(下) 本文主要介绍了使用天文方法计算二十四节气的算法,具体来说是计算二十四节气的天文坐标值。计算二十四节气的坐标需要进行章动修正,因为地球的章动会影响天球上...
《天文算法》是由许剑伟翻译的一本详细探讨天文计算的著作,这本书为软件开发者提供了构建精确天文模型的理论基础和实用方法。书中涵盖了多个关键的天文计算领域,旨在帮助读者理解和实现各种天文现象的计算。 第1...
"position_by_sun.rar_vc 计算_天文_天文航海计算"这个压缩包文件,显然包含了一个使用Visual C++(简称VC)编写的程序,用于进行天文航海中的位置估算。这个程序的核心功能是基于对太阳高度的观测数据,结合船舶的...
标题“prgram test.zip_天文_计算天文相角”表明这是一个与天文相角计算相关的程序,可能是用Fortran语言编写(因为压缩包中的文件名prgram test.f90是Fortran程序的典型命名格式)。让我们深入探讨一下这个主题。 ...
本文深入探讨了《太阳能辐射计算_第一讲太阳能中天文参数的计算》这一主题,聚焦于日地距离、太阳赤纬角及时差等关键天文参数的计算方法,这些参数对太阳能系统的效率评估与设计具有重大意义。 #### 日地距离的精确...
计算新月时间可调用oyyq.calendar.util.NewMoonCalculator.getJulianDayInYearAndMonthForNewMoon(int, int)计算得到。 以上得到的时间是以地球时计算的儒略日,如需转换成实际时间,还需要转成UTC时间再转成日历...
其中儒略日是一种从公元前4713年1月1日开始算起的连续计时系统,至今仍被广泛使用于天文学中,对于精确计算天文事件发生的具体时间至关重要。 在《天文算法英文版——Jean Meeus》一书中,读者可以找到许多与天文学...
在天文学中,计算方位角偏差对于精确定位天体、导航和天文观测至关重要。MATLAB作为一种强大的数值计算和数据可视化工具,被广泛用于天文学领域的软件开发。本项目旨在利用MATLAB开发一个程序,能够计算天体的方位角...
5. **程序实现**:ASTRO.FOR 文件很可能是用FORTRAN编程语言编写的天文辐射计算程序,这是一种早期广泛用于科学计算的语言。此类程序通常包含复杂的算法来处理上述的各种因素,并可能提供用户界面供输入特定的气象和...
随着天文观测技术的发展,我们获取的数据量和复杂度以指数级增长,这使得传统的数据处理和分析方法不再适用。计算机技术的介入,不仅提升了我们处理这些数据的能力,而且为天文学研究带来了全新的视角和工具,让我们...
该软件在天文学的数据处理与展示上提供了精准的计算和丰富的功能,满足了用户对于天文学知识的追求与日常的天文现象查询需求。 许剑伟老师以其深厚的专业知识和不懈的研究精神,对《寿星天文历》进行了持续的优化与...
航天器自主天文导航原理与方法,本人花钱买的,现在共享出来给大家
寿星万年历是一款采用现代天文算法制作的农历历算程序,含有公历与回历信息,可以很方便的...在过去几百年中,寿星万年历的误差是非常小的,节气时刻计算及日月合朔时刻的平均误差小于1秒,太阳坐标的最大可能误差为
天文导航是利用天体(如星星、太阳、月亮)的位置来确定位置的方法。这种导航方式历史悠久,但在现代科技的推动下,天文导航已经发展到利用高精度的星敏感器来快速识别和追踪天体,从而计算出载体的位置。 惯性导航...
本文档《图解球体表面积和体积正确计算方法及计算公式.doc》详细解析了球体表面积和体积的计算方法,并通过图解的方式加深理解。 首先,我们来探讨球体表面积的正确计算方法。球体表面积的计算方法较为直观,其原理...
天文测量是地球科学中的一种基础方法,主要通过观测宇宙中的天体来确定地面点的位置和时间。其关键任务包括测定天文坐标(天文经度与天文纬度)、天文方位角和时刻。天文测量的结果对于大地测量至关重要,它们提供了...
- **章节内容**:本章讨论了在进行天文计算时精度的重要性以及如何提高计算精度的方法。 - **关键知识点**: - 不同类型的数值误差及其来源。 - 提高计算精度的技术,如使用高精度的数据类型等。 #### 3. 插值 - ...