Chinese Calendar Algorithm - Year 1901 to 2100 - v4.15, by Herong Yang
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.
Here is a picture of the formula:
Table of Contents
►Chinese Calendar Background Information
Astronomical Bases of Calendars
Equinoxes, Solstices and Seasons
►Position of Sun - Longitude Angle
Chinese Calendar Algorithm and Program
Chinese Calendars: Year 1901 to 1910
Chinese Calendars: Year 1911 to 1920
Chinese Calendars: Year 1921 to 1930
Chinese Calendars: Year 1931 to 1940
Chinese Calendars: Year 1941 to 1950
Chinese Calendars: Year 1951 to 1960
Chinese Calendars: Year 1961 to 1970
Chinese Calendars: Year 1971 to 1980
Chinese Calendars: Year 1981 to 1990
Chinese Calendars: Year 1991 to 2000
Chinese Calendars: Year 2001 to 2010
Chinese Calendars: Year 2011 to 2020
Chinese Calendars: Year 2021 to 2030
Chinese Calendars: Year 2031 to 2040
Chinese Calendars: Year 2041 to 2050
Chinese Calendars: Year 2051 to 2060
Chinese Calendars: Year 2061 to 2070
Chinese Calendars: Year 2071 to 2080
Chinese Calendars: Year 2081 to 2090