Position of Sun - Longitude Angle

Chinese Calendar Algorithm - Year 1901 to 2100

∟Chinese Calendar Background Information

∟Position of Sun - Longitude Angle

This section provides a tutorial example on how revise and run Keith's QBasic program to calculate the position of the Sun, including longitude angle, using a formula for years between 1950 and 2050.

In the previous tutorial, I said there is no easy ways to calculate Equinoxes and Solstices. But I actually found a simple QBasic program written by Keith Burnett "Position of the Sun" at http://www.stargazing.net/kepler/sun.html.

Keith's program is based on a formula published in "1996 Astronomical Almanac" that can calculates Sun's longitude angle and position to an accuracy of 0.01 degree between the years 1950 and 2050.

Here is Keith's program source code:

 
'*********************************************************
'   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
'*********************************************************

Since the above source code was written for a very old version of QBasic, I have to revise and run it with QB64 from https://www.qb64.org on my macOS computer.

'*********************************************************
'   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
'*********************************************************

Here is the first test session of Position-Of-Sun-QB64.bas:

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

The output matches exactly the example provided by Keith.

Now we can validate Keith's program with some examples from Fred's Equinoxes and Solstices data.

Example 1 - Spring Equinox in 2001, the first date/time of "2001= 20 13:31, 21 07:38, 22 23:05, 21 19:22". The output show that the Sun is at longitude of 0.01 degrees, within the expected error range.

  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

Example 2 - Summer Solstice in 2021, the second date/time of "2021= 20 09:37, 21 03:32, 22 19:21, 21 15:59". The output show that the Sun is at longitude of 90.01 degrees, within the expected error range.

  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

Example 3 - Fall Equinox in 2060, the third date/time of "2060= 19 20:37, 20 13:44, 22 05:47, 21 03:00". The output show that the Sun is perfectly at longitude of 180.00 degrees.

  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

Example 4 - Winter Solstice in 2100, the fourth date/time of "2100= 20 13:04, 21 05:32, 22 22:00, 21 19:51". The output show that the Sun is at longitude of 270.85 degrees, with an error of 0.15 degrees. This is expected, since Keith's program only is valid only between the years 1950 and 2050.

  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

If you want to see the formula provided in "1996 Astronomical Almanac", you can visit its PDF image at https://babel.hathitrust.org/cgi/pt?id=mdp.39015036953076.