Using a list of refraction values for calculating sunrises and sunsets

[roe-dl]

It is a speed issue. I want to calculate sunrise and sunset for a whole year, using measurements of temperature and air pressure to adjust refraction.

I did that with PyEphem with a for loop over the days of the year. First I calculated the rising (setting) without adjustment. Then I looked up temperature and pressure in a database. Then I calculated the rising (setting) again, this time using temperature and pressure. Unfortunately, this approach does not work with Skyfield for speed reasons. It is about 125 times slower than with PyEphem. A small test script (see below for the script) illustrates that:

Skyfield range(365)      : 4.647 s, CPU 4.632 s
Skyfield single statement: 0.034 s, CPU 0.034 s
PyEphem  range(365)      : 0.072 s, CPU 0.072 s

How to solve that issue?

The first step of my calculation I can put into one single Skyfield statement. Then I can look up the database for all the timestamps calculated. That all works in reasonable time. I am sure I can calculate a list of horizons out of the temperatures and pressures using earthlib.refraction within one single statement.

Then I have:

• a Time holding a list of rising (setting) times
• an array of refraction values
• and of course the array of y values

Code snippet in principle:

...
t, y = almanac.find_risings(obs, sun, t0, t1)
temp = []
pressure = []
for ti in t.ut1:
    reply = databaselookup(ti)
    temp.append(reply.temperature)
    pressure.append(reply.pressure)
temp = np.array(temp)
pressure = np.array(pressure)
refr = earthlib.refraction(-body_radius_degrees, temp, pressure)
...

But then?

---

This is the script showing the time consumption:

#!/usr/bin/python3

LAT = 50.0
LON = 15.0

import time
import ephem
from skyfield import almanac
from skyfield.api import load, wgs84

ts = load.timescale()
eph = load('de440.bsp')

def observer():
    """ location of the observer """
    here = wgs84.latlon(LAT,LON)
    return eph['Earth'] + here

def loop1():
    time_t0 = time.time()
    time_t2 = time.thread_time_ns()*0.000000001
    for i in range(365):
        t0 = ts.utc(2025,3,21)
        t1 = ts.utc(2025,3,22)
        obs = observer()
        sun = eph['Sun']
        tri, yr = almanac.find_risings(obs, sun, t0, t1)
    time_t3 = time.thread_time_ns()*0.000000001
    time_t1 = time.time()
    return time_t1-time_t0,time_t3-time_t2

def loop2():
    time_t0 = time.time()
    time_t2 = time.thread_time_ns()*0.000000001
    t0 = ts.utc(2025,1,1)
    t1 = ts.utc(2026,1,1)
    obs = observer()
    sun = eph['Sun']
    tri, yr = almanac.find_risings(obs, sun, t0, t1)
    time_t3 = time.thread_time_ns()*0.000000001
    time_t1 = time.time()
    return time_t1-time_t0,time_t3-time_t2

def loop3():
    time_t0 = time.time()
    time_t2 = time.thread_time_ns()*0.000000001
    for i in range(365):
        obs = ephem.Observer()
        obs.date = '2025/3/21'
        obs.lat = LAT
        obs.lon = LON
        sun = ephem.Sun()
        tr = obs.next_rising(sun)
    time_t3 = time.thread_time_ns()*0.000000001
    time_t1 = time.time()
    return time_t1-time_t0,time_t3-time_t2

print('Skyfield range(365)      : %.3f s, CPU %.3f s' % loop1())
print('Skyfield single statement: %.3f s, CPU %.3f s' % loop2())
print('PyEphem  range(365)      : %.3f s, CPU %.3f s' % loop3())

[brandon-rhodes]

That's an interesting question! Which means, alas, that I don't have a quick answer; but I figured I should go ahead and at least add some kind of comment here, so you at least know I did see the question.

My guess is that you should compute the angle of the horizon over all of the mornings of the year; then maybe create a function that uses NumPy.interp() to deliver those horizon values very quickly when given an array of input times; then write your own version of almanac.risings_and_settings() that uses that variable horizon value for its h.

Or maybe you could do something else: write your own version of risings_and_settings(), but one that passes the temperature and pressure to altaz(), so that you don't even need to adjust the horizon differently for every day.

Or you could ask for all sunrises and sunsets using the fast Skyfield routine, then figure out how far above / below the horizon each one is when you apply temperature and pressure, then adjust the time of sunrise by computing the rate at which the Sun is moving vertically.

Have you measured how much difference the temperature and pressure make for a typical sunrise you are looking at — is the difference a few dozen seconds, or a few minutes, typically?