太阳的位置和经度

∟太阳的位置和经度

本節介紹了計算太陽位置的公式和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