Algorithm to calculate sun/moon rise/set with example

For the last few weeks I’ve been trying to write some code to calculate the rise and set times for the sun and moon for a given day. I’ve been following an explanation of the algorithm on these pages:

http://www.stjarnhimlen.se/comp/ppcomp.html
http://www.stjarnhimlen.se/comp/riset.html

and I just can’t seem to get the moon part to work right. The sun part of it seems to be working okay and I’ve verified most of the moon code against the tutorial on the first page but the moon example doesn’t go all the way to the end and that is where I’m doing something wrong. I just can’t seem to figure out what it is.

I was hoping the author of the page had more details but he hasn’t been very helpful so far. If I just had a tutorial that went all the way through I could figure it out. This is driving me completely crazy - the problem seems to be somewhere in about 10 lines of code but I haven’t puzzled it out yet.

Does anyone know where I can find a good reference that explains the algorithm from beginning to end? I’ve done some Google searching but haven’t turned up anything that quite fits the bill.

I have run across code that uses some other algorithm but I’ve gotten almost all the way there with this algorithm and I don’t want to totally give up yet. =)

Thanks!

Well, this is not something we can do in a few minutes but I will point out what is probably obvious to you already. To get the rise and set times you need to calculate the declination and Hour Angle coordinates first and once you have them the rest is trivial. Did you check your results against the published Almanac? Am I correct in assuming your problem is in computing these magnitudes? I am not sure why you are doing this but there are programs out there which will calculate this for you.

Thanks for the reply - some extra details:

I am doing this as part of a weather monitoring application I’m working on as a hobby. I have a bunch of sensors around the house (temperature, humidity, atmospheric pressure, wind speed and direction, etc) and my application reads the sensor values for display. As part of this application I want to list some astronomical statistics for the day like sun rise/set, moon rise/set, and length of day.

I did some Google searching for algorithms and I came across the pages I mentioned in the OP. Using the web pages and a sample application written by the author as a guide I was able to write my own code to calculate sun rise/set times and the length of the day.

Once I got the sun part working I moved on to the moon calculations which are an extension of the sun calculations. Unfortuantely the author didn’t have any sample code for the moon and the “walkthrough” only goes as far as calculating the right ascension and the declination. I’ve compared my code to the walkthrough and I am getting the correct right ascension and declination for the date and location the author used as an example.

The problem now is the gap between having the position information and getting the rise/set times. The author says that it works the same as for the sun except that you need to iterate using the last computed time as the starting point for the next calculation until the times do not change significantly. Somewhere in here I am doing something wrong because the values I get are nowhere near the values given by an almanac. I’ve been staring at it for quite some time I can’t seem to get anywhere.

I’m hoping that somebody knows of another web page that has a walkthrough of the rest of the algorithm and/or some sample I could compare against. I mailed the author of the web pages but unfortunately he hasn’t been very helpful so far.

I know that this is a nasty problem and I don’t expect anyone to write it for me - I’m hoping somebody knows of a guide I can use to finish this up.

If all else fails I found some javascript code that implements a different algorithm that I could translate but I’ve come so far with this one I don’t want to give up. :slight_smile:

Computing dec and RA for the moon is the difficult part. If you have that right then calculating sunrise and sunset is trivial.

You can check your results against http://aa.usno.navy.mil/
In fact if you just want the end data you can get it from http://aa.usno.navy.mil/data/docs/RS_OneYear.html
You can get astronomical data from http://aa.usno.navy.mil/data/docs/WebMICA_2.html
http://aa.usno.navy.mil/data/docs/celnavtable.html

Anyway, if you want to continue doing all the raw computing then confirm if you are computing your moon coordinates right and then we can go from there. If that is the case it should be easy to compute rise & set times.

BTW, for such an application you really do not need any degree of accuracy and rough data with interpolation would suffice.

Just getting the table from the USNO seems the easiest way as you can get a year at a time.

Sailor has the site I use, but I know there is a book that gives nice, simple calculations for these. I’ve used it at work, 6 or 7 years ago. Darned if I can remember what it is though. :rolleyes: Have you searched Amazon? They can’t have too many books with algorithms for calculating astrodynamic quantities.

Well, there is no simple way to compute the Moon’s coordinates and any degree of precision which would take into account planetary perturbations I think is out of the reach of amateurs (although I think for the purposes of the OP no degree of precision is needed).

What I have used in the past for navigation was the USNO computer almanac which listed coeficients of chebyschev polynomials which allowed the calculation of the variables. I also have an old DOS program DAVIS INSTRUMENTS - PC ASTRO-NAVIGATOR which does the Almanac computations with pretty good accuracy, enough for navigation For the purposes of the OP probably a precision of a minute or two of time is enough.

First of all, are you also factoring in the perturbations and topocentric position of the Moon?

If so, then you’re having the same problem I am, using information from the same source. In my case, it’s for an astronomy program which seems to be having trouble nailing down eclipses to within less than a day or so. Reversing the signs of all the perturbations seems to help, but still doesn’t get it exactly right. Does your program also do eclipses? If we’re getting the same problems then perhaps there is a typo on the page someplace.

Yup - I’m factoring in the perturbations and calculating the topocentric position. If I compare against the test case in the tutorial I get the same topRA and topDecl values so it seems like I have it right - at least for that one date. :slight_smile:

I did a quick check of the values against the USNO values and they seem to be accurate enough for my needs. I need to look more closely just to be safe.

So it seems like I’m not converting from position to time correctly. The document says it should be exactly the same as the sun conversion but it isn’t working.

If you are obtaining dec and RA close enough then obtaining rise and set times are easy. To what degree of precsision do you need?

With your screen name I hope you do realise with the moon you need to take parallax into account?

The USNO calculates moonrise and moonset as:

Horizontal parallax ranges roughly between 54 and 61 minutes of arc but is of opposite sign as semidiameter and augmentation. All things considered the total correction for these concepts ranges from 40.0’ to 44.5’ Both semidiameter (SD) and horizontal parallax (HP) are functions of the same variable (distance of the moon to the earth) and can easily be computed as a single formula.

The formula I use for semidiameter plus augmentation is 0.2733 * HP

HP can easily be calculated as a function of the moon’s distance.

I hope this helps

This page, for the purpose of computing moonrise and moonset, rather than calculate the actual semidiameter and parallax proposes using an average oif 57’:

As I say, if you have dec and RA correctly computed, caculating rise and set times is quite trivial by just using the generic formula for H:

Sin(Hc) = Sin(lat) * Sin(dec) + Cos(lat) * Cos(dec) Cos(LHA)

Here is what I’ve been doing - I’ll try your latest suggestion and report back.

Date 3/27/2004
Time 0:00 UT
Lat 42.8 = 42 deg 48 ’
Long -71.1 = -71 deg 6 ’

My application:

RA 4.9342308084187958
LST 7.577966726399997
Decl 25.026935194582265

USNO:

RA 4 51 39.20 = 4.860888
LST 7 34 39.4250 = 7.577618
Decl 25 06 45.07 = 25.1125194

Based on this I do the following:

Altitude = -0.833

dTimeMoonInSouth = CorrectHours((RA - LST - Long) / 15.0); // CorrectHours brings the result into the range of [0,24)

dLocalHourAngle = (SinDegrees(Altitude) - SinDegrees(Lat) * SinDegrees(Decl)) / (CosDegrees(Lat) * CosDegrees(Decl));

dLocalHour = AcosDegrees(dLocalHourAngle) / 15;

dRiseTime = CorrectHours(dTimeMoonInSouth - dLocalHour);
dSetTime = CorrectHours(dTimeMoonInSouth + dLocalHour);

dRiseTime = 20.762840816711275
dSetTime = 12.36466106089123

USNO_Rise = 13:57 = 13.95
USNO_Set = 22:03 = 22.05

Something ain’t right. =)

I notice you are using the dec and RA at 00:00 UT for your calculations. I do not know what kind of accuracy you need but the coordinates of the moon change fast enough that you would have to use the coordinates at the time of moonrise or moonset. The moon moves very fast in the sky.

I don’t have time to go over the whole thing right now but you may want to look at that. If you post an example with the intermediate calculated values I can go over it and see where we diverge.

My quick calculations give:

moonrise at 42º 48’ N, 71º 6’W , 27 march 13:57 UT
Moon GHA 313º 39.2’, dec 26º 33.4’ N
Works out correctly for me (H=0)

I have never had any problem with this type of thing and I often used to plot a line of position by observing with naked eye when the sun first touched the horizon and I was always surprised that the result was pretty good in spite of the tricky refraction at such low altitude.

The moon moves fast enough on the celestial grid that you need to recalculate its RA and dec. continually and cannot just use data which is several hours out of time. From 00:00 UT to 13:57, the actual time of moonrise, the moon has moved several degrees in RA and more than one degree in declination. It depends on the precision you need but such an error would make the calculations completely useless for navigation. Since you do not tell us what kind of precision you are looking for it is difficult to give useful ideas.

I do not think there is any way to calculate the moon rise and set times directly with a single calculation and you need to iterate. What I would implement is an algorithm which would determine moonrise/set by succesive approximations. You can start the iteration by seeding it with the time of last similar event and add the average time between events. You can get a more precise number but using the rule of thumb that moonrise happes 50 minutes later the following day, if yesterday’s moon rise was at 13:57 start today with 13:57 + 24:50 = 14:47. Feed 14:47 into the loop and let it converge which it should do rather quickly if all you need is a precision of one minute of time. In fact the USNO gives 14:43 so you see a couple of iterations would do it.

if you had to do it by hand then you could just use the coordinates (dec, RA) of the moon for the first approximation because it would only be off by some minutes of time. But you cannot use data which is hours of time off.

So, with all this, the simplest way (and least precise) to calculate moonrise would be:

  1. Take yesterday’s moon rise time and add 24:50 (or a better approximation), call it t.
  2. Compute RA & dec. for moon at that time which will not be recalculated. The same values will be used in succesive iterations as they are close enough
  3. Compute LHA as a function of t
  4. Compute H as a function of (lat, dec, LHA)
  5. If ABS(H)< acceptable error then t is the answer = END
    . . . Else IF (H<0) then make t=t+1 and GOTO (3)
    . . . but IF (H>0) then make t=t-1 and GOTO (3)

Again, it all depends on the precision you need. For navigation you need to take into account pretty much everything. Time has to be with a precision greater than one second and you need to account for refraction, parallax, augmentation etc. But if you just want to determine moonrise for general weather and astronomical information one minute of time is more than sufficient precision.