The Museum of HP Calculators

HP Articles Forum

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

## 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)
```

```                           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
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
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
:
............................................
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:

• Should you be more interested in the equation of time just remove the call to subroutine D in the main program. The result is in RAD and may be transformed to hours using subroutine E.
• The subroutine F isn't really needed. It may be used to check the approximation used in subroutine B. These two programs are (not exactly) inverse to each other.
• I organized the program into smaller modules that can be used and tested independently. However the program can be unfolded to save some space.
• In the listing below I've added the status of the stack for each step to ease reading the program.
• Initially the program was written for the HP-11C but I preferred using letters for the registers. That's why I have selected the HP-32SII. However it's straightforward to transform that program to the HP-41 or HP-15C.

## Annotated Listing

```                    .___________________.___________________.___________________.___________________.___________________.
|         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

## Registers:

```
00   K = 0.016722  numerical orbital eccentricity of the earth
01   E = 23.45     obliquity of the ecliptic
02   P = 78.5      perihelion (relating to vernal equinox)
03   U = 8:30'30"  longitude (Zurich)
04   V = 47:24'36" latitude (Zurich)

```

## Program for HP-41C:

``` 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 -
30 RCL 00              60 +
```

## Appendix:

```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
```
```                                    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.