Sunrise and Sunset
Message #1 Posted by Thomas Klemm on 26 Oct 2010, 6:52 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", 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
Formulae
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.
Constants:
- K = 0.016722 numerical orbital eccentricity of the earth
- E = 23.45 obliquity of the ecliptic
- P = 78.5 perihelion (relating to vernal equinox)
- U = 8:30'30" longitude (Zurich)
- V = 47:24'36" latitude (Zurich)
auxiliary quantity:
1 + K
L2 = -------
1 - K
Initialization: (all angles must be stored in RAD)
- 0.016722 STO K
- 23.45 ->RAD STO E
- 78.5 ->RAD STO P
- 8.3030 ->HR ->RAD STO U
- 47.2436 ->HR ->RAD STO V
- 1 RCL+ K 1 RCL- K / SQRT STO L
Program for HP-32SII:
............................................
: 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
Appendix:
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: 27 Oct 2010, 4:15 p.m. after one or more responses were posted
| 
Go back to the main exhibit hall