中国农历二百年算法及年历 - 和荣笔记 - v4.15, by 杨和荣
太阳的位置和经度
本节介绍了计算太阳位置的公式和Keith的QBasic程序,修改后可以在QB64上运行。该程序可以计算验证二分点和二至点的日期和时间,在1950年至2050年之间的误差为0.01度。
前面我讲到没有简单的方法计算二分点和二至点。 但后来我发现了 Keith Burnett 编写的 QBasic 简单程序,发表于 "太阳的位置"( http://www.stargazing.net/kepler/sun.html)。
Keith的程序采用了《1996年天文年鉴》书中的 计算太阳的经度和位置的公式, 精度在1950年至2050年之间的0.01度。
下面是Keith的程序的源代码:
'********************************************************* ' This program will calculate the position of the Sun ' using a low precision method found on page C24 of the ' 1996 Astronomical Almanac. ' ' The method is good to 0.01 degrees in the sky over the ' period 1950 to 2050. ' ' QBASIC program by Keith Burnett (kburnett@geocity.com) ' ' ' Work in double precision and define some constants ' DEFDBL A-Z pr1$ = "\ \#####.##" pr2$ = "\ \#####.#####" pr3$ = "\ \#####.###" pi = 4 * ATN(1) tpi = 2 * pi twopi = tpi degs = 180 / pi rads = pi / 180 ' ' Get the days to J2000 ' h is UT in decimal hours ' FNday only works between 1901 to 2099 - see Meeus chapter 7 ' DEF FNday (y, m, d, h) = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 _ + 275 * m \ 9 + d - 730531.5 + h / 24 ' ' define some arc cos and arc sin functions and a modified inverse ' tangent function ' DEF FNacos (x) s = SQR(1 - x * x) FNacos = ATN(s / x) END DEF DEF FNasin (x) c = SQR(1 - x * x) FNasin = ATN(x / c) END DEF ' ' the atn2 function below returns an angle in the range 0 to two pi ' depending on the signs of x and y. ' DEF FNatn2 (y, x) a = ATN(y / x) IF x < 0 THEN a = a + pi IF (y < 0) AND (x > 0) THEN a = a + tpi FNatn2 = a END DEF ' ' the function below returns the true integer part, ' even for negative numbers ' DEF FNipart (x) = SGN(x) * INT(ABS(x)) ' ' the function below returns an angle in the range ' 0 to two pi ' DEF FNrange (x) b = x / tpi a = tpi * (b - FNipart(b)) IF a < 0 THEN a = tpi + a FNrange = a END DEF ' ' Find the ecliptic longitude of the Sun ' DEF FNsun (d) ' ' mean longitude of the Sun ' L = FNrange(280.461 * rads + .9856474# * rads * d) ' ' mean anomaly of the Sun ' g = FNrange(357.528 * rads + .9856003# * rads * d) ' ' Ecliptic longitude of the Sun ' FNsun = FNrange(L + 1.915 * rads * SIN(g) + .02 * rads * SIN(2 * g)) ' ' Ecliptic latitude is assumed to be zero by definition ' END DEF ' ' ' CLS ' ' get the date and time from the user ' INPUT " year : ", y INPUT " month : ", m INPUT " day : ", day INPUT "hour UT : ", h INPUT " minute : ", mins h = h + mins / 60 d = FNday(y, m, day, h) ' ' Use FNsun to find the ecliptic longitude of the ' Sun ' lambda = FNsun(d) ' ' Obliquity of the ecliptic ' obliq = 23.439 * rads - .0000004# * rads * d ' ' Find the RA and DEC of the Sun ' alpha = FNatn2(COS(obliq) * SIN(lambda), COS(lambda)) delta = FNasin(SIN(obliq) * SIN(lambda)) ' ' Find the Earth - Sun distance ' r = 1.00014 - .01671 * COS(g) - .00014 * COS(2 * g) ' ' Find the Equation of Time ' equation = (L - alpha) * degs * 4 ' ' print results in decimal form ' PRINT PRINT "Position of Sun" PRINT "===============" PRINT PRINT USING pr2$; " days : "; d PRINT USING pr1$; "longitude : "; lambda * degs PRINT USING pr3$; " RA : "; alpha * degs / 15 PRINT USING pr1$; " DEC : "; delta * degs PRINT USING pr2$; " distance : "; r PRINT USING pr1$; "eq time : "; equation END '*********************************************************
由于这个源代码是用非常老的Qbasic语言编写的,我做了一些修改后, 才能在QB64(https://www.qb64.org) 上运行:
'********************************************************* ' Position-Of-Sun-QB64.bas ' ' This program will calculate the position of the Sun ' using a low precision method found on page C24 of the ' 1996 Astronomical Almanac. ' ' The method is good to 0.01 degrees in the sky over the ' period 1950 to 2050. ' ' QBASIC program by Keith Burnett (kburnett@geocity.com) ' HY: Revised for QB64 by Herong Yang (herong_yang@yahoo.com) ' ' Work in double precision and define some constants ' DEFDBL A-Z ' HY: Need to declare global variables to be SHARED ' DIM SHARED pr1$ DIM SHARED pr2$ DIM SHARED pr3$ DIM SHARED pi DIM SHARED tpi DIM SHARED twopi DIM SHARED degs DIM SHARED rads pr1$ = "\ \#####.##" pr2$ = "\ \#####.#####" pr3$ = "\ \#####.###" pi = 4 * ATN(1) tpi = 2 * pi twopi = tpi degs = 180 / pi rads = pi / 180 ' HY: Need to generate some junk lines first, since the QB64 only ' HY: show the lower half of the screen ' SCREEN _NEWIMAGE(800, 660, 256) CLS FOR x = 1 TO 21 PRINT "..." NEXT x ' ' get the date and time from the user ' INPUT " year : ", y INPUT " month : ", m INPUT " day : ", day INPUT "hour UT : ", h INPUT " minute : ", mins h = h + mins / 60 d = FNday(y, m, day, h) ' ' Use FNsun to find the ecliptic longitude of the ' Sun ' lambda = FNsun(d) ' ' Obliquity of the ecliptic ' obliq = 23.439 * rads - .0000004# * rads * d ' ' Find the RA and DEC of the Sun ' alpha = FNatn2(COS(obliq) * SIN(lambda), COS(lambda)) delta = FNasin(SIN(obliq) * SIN(lambda)) ' ' Find the Earth - Sun distance ' HY: Need to recalculate g instead of borrow it from FNsun() ' g = FNrange(357.528 * rads + .9856003# * rads * d) r = 1.00014 - .01671 * COS(g) - .00014 * COS(2 * g) ' ' Find the Equation of Time ' HY: Need to recalculate L instead of borrow it from FNsun() ' L = FNrange(280.461 * rads + .9856474# * rads * d) equation = (L - alpha) * degs * 4 ' ' print results in decimal form ' PRINT PRINT "Position of Sun" PRINT "===============" PRINT PRINT USING pr2$; " days : "; d PRINT USING pr1$; "longitude : "; lambda * degs PRINT USING pr3$; " RA : "; alpha * degs / 15 PRINT USING pr1$; " DEC : "; delta * degs PRINT USING pr2$; " distance : "; r PRINT USING pr1$; "eq time : "; equation END ' HY: Need to move all custom functions to the end ' HY: Need to change DEF statements in to FUNCTION statements ' ' Get the days to J2000 ' h is UT in decimal hours ' FNday only works between 1901 to 2099 - see Meeus chapter 7 ' 'DEF FNday (y, m, d, h) = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 _ ' + 275 * m \ 9 + d - 730531.5 + h / 24 ' FUNCTION FNday (y, m, d, h) FNday = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 _ + 275 * m \ 9 + d - 730531.5 + h / 24 END FUNCTION ' ' define some arc cos and arc sin functions and a modified inverse ' tangent function ' FUNCTION FNacos (x) s = SQR(1 - x * x) FNacos = ATN(s / x) END FUNCTION FUNCTION FNasin (x) c = SQR(1 - x * x) FNasin = ATN(x / c) END FUNCTION ' ' the atn2 function below returns an angle in the range 0 to two pi ' depending on the signs of x and y. ' FUNCTION FNatn2 (y, x) a = ATN(y / x) IF x < 0 THEN a = a + pi IF (y < 0) AND (x > 0) THEN a = a + tpi FNatn2 = a END FUNCTION ' ' the function below returns the true integer part, ' even for negative numbers ' FUNCTION FNipart (x) FNipart = SGN(x) * INT(ABS(x)) END FUNCTION ' ' the function below returns an angle in the range ' 0 to two pi ' FUNCTION FNrange (x) b = x / tpi a = tpi * (b - FNipart(b)) IF a < 0 THEN a = tpi + a FNrange = a END FUNCTION ' ' Find the ecliptic longitude of the Sun ' FUNCTION FNsun (d) ' ' mean longitude of the Sun ' L = FNrange(280.461 * rads + .9856474# * rads * d) ' ' mean anomaly of the Sun ' g = FNrange(357.528 * rads + .9856003# * rads * d) ' ' Ecliptic longitude of the Sun ' FNsun = FNrange(L + 1.915 * rads * SIN(g) + .02 * rads * SIN(2 * g)) ' ' Ecliptic latitude is assumed to be zero by definition ' END FUNCTION '*********************************************************
下面是首次运行测试的结果:
herong$ cd qb64 herong$ ./qb64 -x Position-Of-Sun-QB64.bas QB64 Compiler V1.4 Beginning C++ output from QB64 code... first pass finished. Translating code... [..................................................] 100% Compiling C++ code into executable... Output: Position-Of-Sun-QB64 herong$ ./Position-Of-Sun-QB64 Display Type: Built-In Retina LCD Resolution: 2880 x 1800 Retina (You see a new QB64 screen) (Enter input data as prompted on the new screen) year : 1997 month : 8 day : 7 hour UT : 11 minute : 0 Position of Sun =============== days : -877.04167 longitude : 134.98 RA : 9.163 DEC : 16.34 distance : 1.01408 eq time : -5.75
测试结果跟Keith的运行的例子完全吻合。
现在我们可以用Fred的二分点和二至点的数据来验证Keith的程序:
例子一 - 2001年的春分,"2001= 20 13:31, 21 07:38, 22 23:05, 21 19:22" 的第一个数据。结果显示太阳的经角为0.01度,符合误差范围:
year : 2001 month : 3 day : 20 hour UT : 13 minute : 31 Position of Sun =============== days : 444.06319 longitude : 0.01 RA : 0.001 DEC : 0.00 distance : 0.99599 eq time : 1432.56
例子二 - 2021年的夏至,"2021= 20 09:37, 21 03:32, 22 19:21, 21 15:59" 的第二个数据。结果显示太阳的经角为90.01度,符合误差范围:
year : 2021 month : 6 day : 21 hour UT : 3 minute : 32 Position of Sun =============== days : 7841.64722 longitude : 90.01 RA : 6.00 DEC : 23.44 distance : 1.01625 eq time : -1.78
例子三 - 2060年的秋分,"2060= 19 20:37, 20 13:44, 22 05:47, 21 03:00" 的第三个数据。结果显示太阳的经角为180.00度,完全吻合:
year : 2060 month : 9 day : 22 hour UT : 5 minute : 47 Position of Sun =============== days :22179.74097 longitude : 180.00 RA : 12.00 DEC : -0.00 distance : 1.00377 eq time : 7.46
例子四 - 2100年的冬至,"2100= 20 13:04, 21 05:32, 22 22:00, 21 19:51" 的第四个数据。结果显示太阳的经角为270.85度,有0.15度的误差。 这在预料之中,因为Keith的程序只能保证在1950到2050年有0.01度的误差。
year : 2100 month : 12 day : 21 hour UT : 15 minute : 59 Position of Sun =============== days :36880.16597 longitude : 270.85 RA : 18.062 DEC : -23.42 distance : 0.98376 eq time : 1.50
如果你对《1996 Astronomical Almanac》书中的公式感兴趣, 可以在 https://babel.hathitrust.org/cgi/pt?id=mdp.39015036953076 提取PDF图像。
下面是公式的截屏:
Table of Contents