太阳的位置和经度

∟太阳的位置和经度

本节介绍了计算太阳位置的公式和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