The Museum of HP CalculatorsHP Articles Forum
[Return to the Index ]
[ Previous | Next ]
Sunrise and Sunset
Posted by Thomas Klemm on 28 Oct 2010, 12:05 p.m.
At the horticultural exhibition Grün 80 I discovered this precision sundial punctual to the minute and still remember exactly how much I was impressed watching the time go by. It is still there if you ever happen to be in the vicinity.
Years later, I attended a lecture by Christian Blatter: "Von den Keplerschen Gesetzen zu einer minutengenauen Sonnenuhr" (From Kepler's laws to an exact sundial). I'm sorry that the paper is only available in German but here's a similar paper in English: "Another design for sundials" by Norbert Hungerbühler. Though I knew about the equation of time and Kepler's law I wasn't aware of the interconnection until then.
Now that we can calculate the position of the sun for any date of the year I thought we could go one step further and calculate sunrise and sunset. However I was quiet disappointed by the result because it was very inaccurate compared to what I got for instance from http://www.sunrisesunset.com. It took me some time to understand why: It all depends on how you define sunrise and sunset. The usual value of the altitude of the center of the solar disc is about -0.83° while I assumed this to be zero. This may lead to an aberration of several minutes. On the other hand you're always late for the sunset if the upper edge just disappears.
The program is probably not the most precise. There are others around that do much better. My goal was to understand what's going on.
Best regards
Thomas
From the paper I mentioned above the formula on page 158 was used in combination with the ansatz: psi(t) = t + d(t).
5
psi(t) = t + 2 sin t K + - sin(2t) K2
4
To make the calculation a bit easier I used the double-angle formula:
sin(2t) = 2 sin(t) cos(t)
This leads to:
5
psi(t) = t + K sin t (2 + - K cos t)
2
Both values K cos t and K sin t can be calculated in one step.
The expression in formula (6) on page 161 for mju(t) is stupendous. However it's far easier to use a suitable combination of polar/rectangular coordinate transformation. Don't tell me now you miss something on your HP-35s.
auxiliary quantity:
1 + K
L2 = -------
1 - K
............................................ : Main Program : : (dd.mm -> sunrise sunset) ............................................ 001 XEQ A 002 309 003 - 004 365 005 / 006 2 007 * 008 PI 009 * 010 XEQ B 011 XEQ C 012 GTO D............................................ : Difference days : : (dd.mm -> days) : : calculates difference of days to March 1st ............................................ A01 LBL A A02 IP A03 12 A04 LASTx A05 FP A06 100 A07 * A08 3 A09 - A10 x<0? A11 + A12 30.6 A13 * A14 0.5 A15 + A16 IP A17 R^ A18 + A19 RTN
............................................ : Position of sun (approximation) : : (t -> psi) : : psi(t) = t + 2 K sin(t) + 5/4 K^2 sin(2t) + ... ............................................ B01 LBL B B02 ENTER B03 ENTER B04 RAD B05 RCL K BO6 P->R B07 5 B06 * B09 2 B10 / B11 2 B12 + B13 * B14 + B15 RTN
............................................ : Transformation to equatorial coordinates : : (psi -> mju delta) : : phi = psi - alpha : : x = cos(phi) = cos(delta) cos(AR) : y = sin(phi) cos(eps) = cos(delta) sin(AR) : z = sin(phi) sin(eps) = sin(delta) : : mju = AR - (t - alpha) ............................................ C01 LBL C C02 RCL- P C03 RCL E C04 x<>y C05 1 C06 P->R C07 RDN C08 P->R C09 R^ C10 R->P C11 x<>y C12 RDN C13 R->P C14 x<>y C15 RDN C16 RDN C17 RCL- P C18 - C19 RTN
............................................ : Sunrise and Sunset : : (mju delta -> Sunrise Sunset) ............................................ D01 LBL D D02 RCL+ U D03 ENTER D04 ENTER D05 R^ D06 TAN D07 RCL V D08 TAN D09 * D10 ACOS D11 + D12 x<>y D13 LASTx D14 - D15 XEQ E D16 x<>y
............................................ : Angle to Hours : : (RAD -> HMS) ............................................ E01 LBL E E02 2 E03 PI E04 * E05 / E06 FP E07 1 E08 + E09 FP E10 24 E11 * E12 ->HMS E13 RTN
............................................ : Test approximative solution of Kepler Equation : : (psi -> t) ............................................ F01 LBL F F02 2 F03 / F04 TAN F05 RCL/ L F06 ATAN F07 2 F08 * F09 ENTER F10 SIN F11 RCL* K F12 - F13 RTN
Let me add a few things:
.___________________.___________________.___________________.___________________.___________________.
| X | Y | Z | T | L |
| dd.mm | | | | |
001 IP | d | | | | dd.mm |
002 12 | 12 | d | | | dd.mm |
003 LASTx | dd.mm | 12 | d | | dd.mm |
004 FP | 0.mm | 12 | d | | |
005 100 | 100 | 0.mm | 12 | d | |
006 * | m | 12 | d | d | |
007 3 | 3 | m | 12 | d | |
008 - | m - 3 | 12 | d | d | |
009 x<0? | m - 3 | 12 | d | d | |
010 + | n | ? | d | d | |
011 30.6 | 30.6 | n | ? | d | |
012 * | 30.6 n | ? | d | d | |
013 0.5 | 0.5 | 30.6 n | ? | d | |
014 + | 30.6 n + 0.5 | ? | d | d | |
015 IP | D | ? | d | d | |
016 R^ | d | D | | | |
017 + | days | | | | |
018 309 | 309 | days | | | |
019 - | days - 309 | | | | |
020 365 | 365 | days - 309 | | | |
021 / | T | | | | |
022 2 | 2 | T | | | |
023 * | 2 T | | | | |
024 PI | pi | 2 T | | | |
025 * | t | | | | |
026 ENTER | t | t | | | |
027 ENTER | t | t | t | | |
028 RAD | t | t | t | | |
029 RCL K | K | t | t | t | |
030 P->R | K cos t | K sin t | t | t | |
031 5 | 5 | K cos t | K sin t | t | |
032 * | 5 K cos t | K sin t | t | t | |
033 2 | 2 | 5 K cos t | K sin t | t | |
034 / | 5/2 K cos t | K sin t | t | t | |
035 2 | 2 | 5/2 K cos t | K sin t | t | |
036 + | 2 + 5/2 K cos t | K sin t | t | t | |
037 * | d(t) | t | t | t | |
038 + | psi | t | t | t | |
039 RCL- P | phi | t | t | t | |
040 RCL E | eps | phi | t | t | |
041 x<>y | phi | eps | t | t | |
042 1 | 1 | phi | eps | t | |
043 P->R | x | r | eps | t | |
044 RDN | r | eps | t | x | |
045 P->R | y | z | t | x | |
046 R^ | x | y | z | t | |
047 R->P | r | AR | z | t | |
048 x<>y | AR | r | z | t | |
049 RDN | r | z | t | AR | |
050 R->P | 1 | delta | t | AR | |
051 x<>y | delta | 1 | t | AR | |
052 RDN | 1 | t | AR | delta | |
053 RDN | t | AR | delta | | |
054 RCL- P | t - alpha | AR | delta | | |
055 - | mju | delta | | | |
056 RCL+ U | tau | delta | | | |
057 ENTER | tau | tau | delta | | |
058 ENTER | tau | tau | tau | delta | |
059 R^ | delta | tau | tau | tau | |
060 TAN | tan delta | tau | tau | tau | |
061 RCL V | theta | tan delta | tau | tau | |
062 TAN | tan theta | tan delta | tau | tau | |
063 * | cos omega | tau | tau | tau | |
064 ACOS | omega | tau | tau | tau | |
065 + | rho | tau | tau | tau | omega |
066 x<>y | tau | rho | tau | tau | omega |
067 LASTx | omega | tau | rho | | omega |
068 - | sigma | rho | | | |
069 XEQ E | set | rho | | | |
070 x<>y | rho | set | | | |
E01 LBL E | rho | set | | | |
E02 2 | 2 | rho | set | | |
E03 PI | pi | 2 | rho | set | |
E04 * | 2 pi | rho | set | | |
E05 / | r | set | | | |
E06 FP | [r] | set | | | |
E07 1 | 1 | [r] | set | | |
E08 + | r' | set | | | |
E09 FP | [r'] | set | | | |
E10 24 | 24 | [r'] | set | | |
E11 * | r" | set | | | |
E12 ->HMS | mm.hhss | set | | | |
E13 RTN | rise | set | | | |
Quote:
it's straightforward to transform that program to the HP-41
01 LBL "SUN" 31 P-R 61 ENTER 02 INT 32 5 62 ENTER 03 12 33 * 63 R^ 04 LASTX 34 2 64 TAN 05 FRC 35 / 65 RCL 04 06 100 36 2 66 TAN 07 * 37 + 67 * 08 3 38 * 68 ACOS 09 - 39 + 69 + 10 X<0? 40 RCL 02 70 X<>Y 11 + 41 - 71 LASTX 12 30.6 42 RCL 01 72 - 13 * 43 X<>Y 73 XEQ 00 14 0.5 44 1 74 X<>Y 15 + 45 P-R 75 LBL 00 16 INT 46 RDN 76 2 17 R^ 47 P-R 77 PI 18 + 48 R^ 78 * 19 309 49 R-P 79 / 20 - 50 X<>Y 80 FRC 21 365 51 RDN 81 1 22 / 52 R-P 82 + 23 2 53 X<>Y 83 FRC 24 * 54 RDN 84 24 25 PI 55 RDN 85 * 26 * 56 RCL 02 86 HMS 27 ENTER 57 - 87 END 28 ENTER 58 - 29 RAD 59 RCL 03 30 RCL 00 60 +
Earth's perihelion and aphelion
Year Perihelion Aphelion
Date Hour Date Hour
2007 January 3 20:00 July 7 00:00
2008 January 3 00:00 July 4 08:00
2009 January 4 15:00 July 4 02:00
2010 January 3 00:00 July 6 11:00
2011 January 3 19:00 July 4 15:00
2012 January 5 00:00 July 5 03:00
2013 January 2 05:00 July 5 15:00
2014 January 4 12:00 July 4 00:00
2015 January 4 07:00 July 6 19:00
2016 January 2 23:00 July 4 16:00
2017 January 4 14:00 July 3 20:00
2018 January 3 06:00 July 6 17:00
2019 January 3 05:00 July 4 22:00
2020 January 5 08:00 July 4 12:00
Orbital Paramameters for the Earth
Long. of
Year Eccentri Obliquity Perihel.
(A.D.) city (degrees) (degrees)
------ -------- --------- --------
1990 .016708 23.4411 282.724
1991 .016707 23.4409 282.741
1992 .016707 23.4408 282.758
1993 .016707 23.4407 282.776
1994 .016706 23.4405 282.793
1995 .016706 23.4404 282.810
1996 .016705 23.4403 282.827
1997 .016705 23.4402 282.844
1998 .016704 23.4400 282.861
1999 .016704 23.4399 282.878
2000 .016704 23.4398 282.895
2001 .016703 23.4396 282.913
2002 .016703 23.4395 282.930
2003 .016702 23.4394 282.947
2004 .016702 23.4392 282.964
2005 .016702 23.4391 282.981
2006 .016701 23.4390 282.998
2007 .016701 23.4389 283.015
2008 .016700 23.4387 283.033
2009 .016700 23.4386 283.050
2010 .016700 23.4385 283.067
2011 .016699 23.4383 283.084
2012 .016699 23.4382 283.101
2013 .016698 23.4381 283.118
2014 .016698 23.4379 283.135
2015 .016698 23.4378 283.152
2016 .016697 23.4377 283.170
2017 .016697 23.4376 283.187
2018 .016696 23.4374 283.204
2019 .016696 23.4373 283.221
2020 .016696 23.4372 283.238
Edited: 1 Nov 2010, 10:20 a.m.
Go back to the main exhibit hall