Inverse of the sunrise equation - finding locations with a given sunrise time on a given day

I'm working on a project for fun where I represent some sleep data geographically. For a given day, I have a date, a time for falling asleep that night, and a time for waking up the next day. The idea is something like this xkcd comic. My attempts at taking the inverse of the sunrise/sunset equations I've found online have not been fruitful.

Basically, I would like to find an equation or algorithm to give a rough location on Earth with a specific sunset time on some specific day and a specific sunrise time the next day. I am not too concerned with accuracy since this is just for fun, so I'm hoping elevation information isn't necessary.

I'm a programmer, so the goal is really to be able to type something like

>>>> locate(date = '2014-08-27', sunset = '10:00 PM', sunrise = '7:30 AM')


and get back

Latitude: 49.887, Longitude: 96.141


I could program it myself if I just understood enough astronomy. Any help is greatly appreciated.

• There wouldn't be a just one location, there would usually be a line across the globe where the Sun sets at a specific time. The length and shape of the line will vary with the date and time. – James Screech Jul 7 '16 at 7:18
• Hmmm... at a given date & time, we should be able to find an equation that tells us where the Sun is (perhaps the Lat/Long co-ordinates of a point in the Tropics, where the Sun is currently directly overhead.) Then the sunrise/sunset points form a great circle centred on that point. (By the way this technique has been used for following bird migrations.) – Andy Jul 7 '16 at 7:56
• @JamesScreech I realize that there is a line across the globe where the sun sets at a specific time, but at each point on that line, the sun will rise at a different time the next day. The goal is to find a location where the sun sets at one specific time and rises at another specific time. – xanxerus Jul 8 '16 at 7:47
• I'm assuming you've seen the wiki page for the Sunrise equation. Note that what you want to determine is $\phi$ and $\lambda$, given $J_{set}$ and $J_{rise}$. I'd suggest putting all the euqations together into one big one that has only those four variables, then, since your equation is transcendental, solve it numerically. – zephyr Aug 23 '16 at 12:52
• OK I'll bite - do you mean local time or UTC? If you specify local time there can often be something like 24 equivalent spots on the earth - one for each time zone. UTC would get rid of the degeneracy. The rise/set pair will work nicely but not in all cases. Also, did you get the two numbers swapped here? Latitude: 96.141, Longitude: 49.887 (latitude greater than 90°) – uhoh Aug 23 '16 at 17:44

To avoid the headaches associated with finding roughly two dozen similar solutions in the zig-zagging time zone boundaries , you might want to stick with UTC. If you want to deal with local time - it's going to get a little crazy if you also try to include daylight savings toggling on and off.

The points on the earth where the center of the sun coincides with the horizon (ignoring topography, oblateness, atmospheric refraction, the finite speed of light and other small effects) is just the circle on the earth where a cone drawn from the Sun intersects a spherical Earth.

above: A sphere inside a cone, from http://mathcentral.uregina.ca.

But the next step is hard - getting from that circle to the lat/lon coordinates on the surface, because the axis is tilted, and because the motion around the sun speeds up and slows down as the Earth revolves closer and farther from the Sun. These can be approximated by periodic functions with some coeficients, and that's what you will find in the math behind the first answer. For the second answer, I'll list two Python solutions - PyEphem and Skyfield. Both are easy to use, but you are separated from the actual math (in one case it's an ephemeris/table. The third answer is really a collection of NASA/JPL routines that are highly regarded but may involve more time for you to get up to speed compared to the Python packages.

This is something you'll have to dig into a bit, but if you like to program, it may be exactly what you're looking for. The Gaisma website is one of my favorites on the internet - easy to use and presents a bunch of information in easy-to-understand graphics. Click around!

I believe that this site uses algorithms from the collection found at this NOAA site. Click around there as well. They provide Excel spreadsheets which contain the algorithms and other resources. The "main" resource is a collection of algorithms published in the book Astronomical Algorithms - Jean Meeus. From that Amazon page you can see that there are many similarly titled books. I'd recommend going to a library if possible, because it's (in my opinion) always good to go to libraries. However parts of these can be found on-line. For example, a few pages shown from the book Astronomical Formulae for Calculators (1988) include an interesting table of contents.

Answer 2: Python packages PyEphem and Skyfield

I'll copy a bit of the text from this answer:

1. The Python package PyEphem has been around and well supported, and is the pythonic reincarnation of XEphem. I haven't used it, but I believe it keeps enough information about orbital parameters at certain epochs to generate an ephemeris, including some gravitational perturbations. In other words, it's much more than planets moving on fixed elliptical orbits around a fixed sun. So I believe it runs without internet connection.

2. I never used it because I was recommended to look at Skyfield and it's exactly what I needed. It downloads a standard JPL ephemeris that you choose, and then just uses it from your hard drive after that. However, in order to deal with leap seconds and other time related effects, it occasionally needs to check the internet for leap second information updates, since these are arbitrary.

I don't know if Skyfield has a mode to avoid that. Actually that's a good question. If you work with a timescale that doesn't have leap seconds, I am not sure if it will run in its current version.

Both Skyfield and PyEphem Python packages have been written and are maintained by @BrandonRhodes.

I've included a simple Python script just to illustrate how Skyfield can be used. If you are comfortable with coding without curly braces, I would strongly recommend you give this a try. It's incredibly powerful and Pythonic.

This is just a starter - you need to add some better housekeeping to detect sunrises vs sunsets, and maybe a more global-type search in some cases. Actually there is some slightly tedious housekeeping necessary to make this work robustly.

note: you can turn on atmospheric refraction using the arguments in the apparent() method. See the Skyfield API documentation for more info, and for a discussion about iterating using Skyfield methods - especially solving for times, see this helpful answer.

def alt_lonlat(lon, lat, t):

topo = earth.topos(lat, lon)

alt, az, dist = topo.at(trise).observe(sun).apparent().altaz() ## apparent() args for atmospheric refraction

return alt.degrees

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as spo

# your example:  '2014-08-27', sunset = '10:00 PM', sunrise = '7:30 AM'

trise  = ts.utc(2014, 8, 27,  7, 30, 0)
tset   = ts.utc(2014, 8, 27, 22,  0, 0)

earth = data['earth']
sun   = data['sun']

zerozero = earth.topos(0.0, 0.0)   # gotta start looking somewhere!

alt, az, dist = zerozero.at(trise).observe(sun).apparent().altaz() ## apparent() args for atmospheric refraction

print "at trise, JD = ", trise.tt
print "at (0N, 0E) Sun's altitude: ", alt.degrees, "azimuth: ", az.degrees
print "at (0N, 0E) Sun's distance (km): ", dist.km

# Find points on equator where sun is on horizon (rise or set) at t=trise

limits   = ((0, 180.), (180, 360.))
lonzeros = []

for a, b in limits:

answer, info = spo.brentq(alt_lonlat, a, b,
args=(0.0, trise),
full_output = True )

if info.converged:
print "limits ", a, b, " converged! Found longitude (deg): ", answer
else:
print "limits ", a, b, "whaaaa?"
lonzeros.append(None)

# make some curves

lats = np.linspace(-60, 60, 13)

longis = []
for lon0 in lonzeros:
lons = []
for lat in lats:

answer, info = spo.brentq(alt_lonlat, lon0-90, lon0+90,
args=(lat, trise),
full_output = True )
if info.converged:
else:
lons.append(None)

lons = [(lon+180)%360.-180 for lon in lons]  # wraparound at +/- 180

longis.append(lons)

plt.figure()

for lons in longis:
plt.plot(lons, lats)

for lons in longis:
plt.plot(lons, lats, 'ok')

plt.xlim(-180, 180)
plt.ylim(-90, 90)

plt.title("at trise, JD = " + str(trise.tt))

plt.show()


As pointed out by @barrycarter in this comment below this answer, the JPL SPICE Kernels are available. I'm not familliar with them but it's what NASA uses so it must be pretty good :)

Appendix:

Here are some screenshots for London, UK (Gaisma from here):

• @barrycarter the question comes from a programmer - if you can improve the wording of the Answer 3: SPICE Kernels section please do! – uhoh Aug 24 '16 at 4:23
• This answer is great! I had not considered the intricacies of time zones, so UTC is the way to go. Answer 2 looks very promising and I just so happen to know Python pretty well so I will try it as soon as I can – xanxerus Aug 24 '16 at 14:56
• OK that's great! If you have python questions you can ask them on Stackoverflow - don't forget to add the skyfield tag. – uhoh Aug 24 '16 at 14:59
• But you might consider up voting and/or accepting this answer if you found it useful, then you can ask a new question if you'd like to go deeper. – uhoh Aug 24 '16 at 15:10

Due to the complexity of the problem (see my comment above) and the fact that you only want a rough answer why not try calculating the times of Sun rise / set at say one degree intervals for the whole Earth and then compare these with the times you have. If one degree is not near enough you could use a smaller interval and also use a range in the comparison (say +/- 5 mins).

Exactly how fine your location spacing and time range is will depend on how accurate you want your answers and both could easily be tweaked until you are satisfied with the results.

This approach may not be the most efficient, but with the speed of modern computer shouldn't take very long.

• This kind of brute force approach is what I was hoping to avoid. I'm more interested in finding a mathematical solution or an efficient approach than I am with what the locations actually are. – xanxerus Jul 8 '16 at 7:44
• I wouldn't really consider this a brute force approach. In the strictest sense yes, but using a lookup table for such a tiny amount of values (even every tenth of a degree only gives you 3600 values) is well within the capabilities of modern computers and used quite often for many tasks (see flipping bit algorithms for example). – zephyr Aug 23 '16 at 12:46
• @zephyr A lookup table would be a great idea for finding a lot of values on the same date, but calculating one value for each date where the data changes from day to day does not seem like something that uses lookup tables efficiently. Am I misunderstanding something? – xanxerus Aug 25 '16 at 18:29

As @Andy says, at any time there is a great circle around the Earth where the Sun is currently at the horizon, either rising or setting (or both at two points in the (ant)arctics). This great circle can be described by the coordinates of the point where the Sun is directly overhead (which is the vector pointing from the Earth's centre to that location on the surface, or to the Sun).

Such a great circle exists for your moment for sunset, and another one exists for your time for sunrise. You are interested in one of the two intersection points of these two great circles.

The two great circles lie in two planes through the Earth's centre, described by two of the vectors mentioned before. The intersection of the two planes is a line through the centre, whose direction is another vector (the normalised cross product of the two normal vectors), which can be converted back to two antipode positions on the Earth's surface. This solution would be rather exact if the Earth were a sphere.