Sunrise and Sunset



#19

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


#20

Nice! very interesting, will dig into the papers tonight, thanks for sharing!

Cheers

Peter

#21

Hallo Thomas,

beeindruckend! Damit deine Arbeit auf Dauer leichter gefunden wird, schlage ich vor, dass du sie in den Bereich "Articles" kopierst.

(Impressing work! I'd suggest you copy your work to the articles section for easier access in the future.)

Grüße,

Walter


#22

Quote:
I'd suggest you copy your work to the articles section

cf. 1014: Sunrise and Sunset

Best regards

Thomas


#23

Thomas --

Good work. Maybe a more-appropriate choice than the Articles Archive would be the Software Library and HP-41 Software Library.

-- KS

Edited: 28 Oct 2010, 11:26 p.m.


#24

Hi Karl

Maybe you're right. But I wouldn't recommend this program to actually calculate sunrise and sunset.
Why would someone do that at all? It is indeed written in the newspaper every day.

I probably would not have posted this at all, unless Geir Isene recently would have drawn my attention
to the module ASTRO 2010 by Jean-Marc Baillard. So I suppose there are a few of you who are interested in the topic as well.

There is an artist, I know. She noticed a spot of sunlight on the floor of the barn of her grandmother.
Every day at the same time, she drew its outline on a large sheet of paper. With different colors. What led to an image over the time.
But it's not about whether that's art. It's about perception.
What happens when we do something every day at the same time? Why does this give a figure eight, symbol of infinity?

Or the days after the winter solstice. Why are the days getting longer again but still the sunrise happens later and later?

21.12                X:       8.2414   Y:      16.3906
22.12 X: 8.2447 Y: 16.3933
23.12 X: 8.2517 Y: 16.4002
24.12 X: 8.2545 Y: 16.4034
25.12 X: 8.2610 Y: 16.4109
26.12 X: 8.2632 Y: 16.4147
27.12 X: 8.2652 Y: 16.4227
28.12 X: 8.2708 Y: 16.4310
29.12 X: 8.2722 Y: 16.4355
30.12 X: 8.2733 Y: 16.4443
31.12 X: 8.2741 Y: 16.4534
1.01 X: 8.2746 Y: 16.4627
2.01 X: 8.2748 Y: 16.4722
3.01 X: 8.2747 Y: 16.4819
4.01 X: 8.2743 Y: 16.4919
5.01 X: 8.2736 Y: 16.5021
6.01 X: 8.2727 Y: 16.5125
7.01 X: 8.2714 Y: 16.5231
8.01 X: 8.2659 Y: 16.5339
9.01 X: 8.2640 Y: 16.5449
10.01 X: 8.2619 Y: 16.5600

And then the question of why I wrote a program for a thirty year old calculator and didn't simply key in the formula into Mathematica?
Most probably because for me it is still the best way to understand these things.

Kind regards

Thomas


#25

Quote:
But I wouldn't recommend this program to actually calculate sunrise and sunset.
Why would someone do that at all? It is indeed written in the newspaper every day.

I don't have a need to use a calculator to compute sunrise and sunset. But a newspaper is simply not adequate so my computer does that calculation for me every day. My computer needs to know when to expect the sunrise and sunset so it can perform tasks relative to those events. It is much easier for it to calculate sunrise and sunset than to fetch a newspaper and read it.

For example, at 0645 my computer starts to gradually bring up a light in my bedroom, but only if the sun is not yet up. Then it turns off the light about 30minutes after sunrise.


#26

You might be interested in Egan's article: Control the World with HP-IL/RS-232/41CX + X10

#27

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.

Cheers

Thomas

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

Edited: 27 Oct 2010, 3:57 p.m.

#28

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 -
29 RAD 59 RCL 03
30 RCL 00 60 +

Edited: 1 Nov 2010, 10:23 a.m. after one or more responses were posted


#29

For those who may be interested, here you can find my "Sun Tracker" program, for the 40G and the 49G. The 49G version also runs on the 50G.

http://www.physique.lu/hpcalculators/sunposition/sunposition.htm

I also hava an 48G version which is not online. If people are interested I will check whether I have more actual versions on my calculators (I think so). I haven't updated the code for some time.


#30

Quote:
I also hava an 48G version

Yes, please post the code here. I'd be interested.

Many thanks in advance

Thomas

#31

Patrick --

That looks like a good opportunity for me to try I/O on the HP-49G and actually use it for something...

And also my HP-48G, if and when the code is available.

-- KS

Edited: 28 Oct 2010, 11:26 p.m.


#32

Hello Thomas, Hello Karl.

Give me some time. I will update the code within the weekend and will let you know here.

#33

I did it immediately, here you can find the actual versions:

http://www.physique.lu/hpcalculators/sunposition/sunposition.htm

In this version the errors during midnight sun and polar night are gone. If you use Urania then you will see that the program uses the same variables for latitude, longitude and time zone. There is one "minor bug" left. The times for sunrise and sunset are only calculated once when you start the program. So if you let the program run for more than a day, they will not be accurate because they will not be updated. But who will do that?

Pay attention that you now have to start the program with the "SUNT" command (for sun tracker).

The 49G versions also runs on the 50G. The 2 last lines of the screen will remain blank. I also added the version for the HP48G series, but did not have the time to check it on my 48G (fuzzy file transfer ...), but it worked fine on the emulator. Just let me know if there is a problem.

The programs are written in system RPL.


#34

Thanks a lot

Thomas

#35

Thomas, if you agree this seems a good one to add to the -ASTRO U/I module, which still has plenty of available space. I can add the alpha prompts to make it more intuitive - or you'd like to?

Any other contribution? Geir?


#36

Hi Ángel

Feel free to do whatever you like with the program. I agree that the UI is rather sparse.

Just be warned: the program uses currently the following formula to calculate the hour angle:

This correction to that equation is used by Patrick's "Sun Tracker":

with the altitude (a) of the center of the solar disc set to -0.833°.

However I doubt that this alone explains the deviation.

Best regards

Thomas


cf.
Sunrise equation


Possibly Related Threads...
Thread Author Replies Views Last Post
  Bearing to sunrise & sunset. r. d. bärtschiger. 2 134 06-27-2002, 12:37 PM
Last Post: andy d

Forum Jump: