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

Question

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.

Answer 1

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.

Answer 1: Astronomical Algorithms

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


from skyfield.api import load
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as spo

data  = load('de421.bsp')
ts    = load.timescale()

# 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:
        lonzeros.append(answer)
        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:
            lons.append(answer)
        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()

Answer 3: SPICE Kernels

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

Answer 2

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.

Answer 3

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.