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) K^{2}

^{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

L^{2}= -------

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*