Missing rise and set events in rising_and_settings?

[lucaberti]

Hello, I post here what I posted also on StackExcange, since this could be me doing something wrong, but it could be a bug as well.

To compute the Andromeda Galaxy Risings and Settings I essentially use this code (in a more colplex code, but this serves the purpose):

from skyfield import api, almanac
from skyfield.api import N,S,E,W, wgs84, load, Star

ts = api.load.timescale()
eph = api.load('de421.bsp')
t0 = ts.utc(2021, 11, 1)
t1 = ts.utc(2021, 11, 13)

#The object
andromeda = Star(ra_hours=(00, 42, 44.330), dec_degrees=(41, 16, 7.50), parallax_mas=6.0, radial_km_per_s=-300, names=("Andromeda"))

#The Location
location = api.wgs84.latlon(46.1927800 * N, 9.0170300 * E, elevation_m=200 )

rise_set = almanac.risings_and_settings(eph, andromeda, location)
riseset_time, isrise = almanac.find_discrete(t0, t1, rise_set)
for riseset_time_u, isrise_u in zip(riseset_time, isrise):
    print(riseset_time_u.utc_iso(), 'Rise' if isrise_u else 'Set')

I do optain this results

2021-11-04T07:51:31Z Set
2021-11-04T10:33:33Z Rise
2021-11-05T07:47:35Z Set
2021-11-05T10:29:37Z Rise
2021-11-06T07:43:39Z Set
2021-11-06T10:25:41Z Rise
2021-11-07T07:39:43Z Set
2021-11-07T10:21:45Z Rise
2021-11-08T07:35:47Z Set
2021-11-08T10:17:49Z Rise
2021-11-09T07:31:51Z Set
2021-11-09T10:13:53Z Rise

Now the problem is that t0 is set on November the 1st, but the rises_and_settings function does not return any result before Novembr 4th.

P.S Transits seem to be correct, though

[brandon-rhodes]

Thanks for opening the issue, as this is something we'll want to improve! In the meantime, to get your application working, here’s what I answered over on the Stack Exchange website:

--

The rising-setting search routine is currently optimized for objects near the celestial equator, like the Sun, Moon, and planets, that spent about half the day in the sky. This assumption lets it search in ¼-day increments, making it more efficient for such objects.

But the Andromeda galaxy and your observer position are both far to the north of the ecliptic, such that the galaxy spends only (if I'm doing the math in my head correctly, given the output you posted) about 4 hours below the horizon each day. It's easy for Skyfield to miss that 4 hour window when it’s stepping forward in ¼-day = 6 hour steps looking for the galaxy to be below the horizon.

Maybe Skyfield should someday have more intelligence around risings and settings! For now, you can tell it to search a tighter collection of points:

rise_set.step_days = 0.1

If you put that before your call to almanac.find_discrete(), I believe you will find that it computes all of the galaxy’s risings and settings!

[lucaberti]

Thank you for the answer. This actually helped. Sorry for the queston: I'm still trying to grab how Skyfield works in its core.

[davidmikolas]

Just some thoughts, without appreciation of other problems it may cause:

If step_days were an explicit argument of risings_and_settings() like horizon_degrees and radius_degree then readers of the documentation (or those using help()) would know this was there and the default was 0.25. Right now it's a hidden assumption.

The curse of writing such an amazing package is that we assume everything is "right" and "perfect" if not explicitly warned otherwise :-)

[bzhaas]

FYI, I encountered this issue when trying compute lunar nodes. The code below dropped all nodes between 2022-05-29 and 2022-09-28:

eph = load('de440s.bsp') 0 = ts.utc(2022, 1, 1) t1 = ts.utc(2023, 1, 1) t, n = almanac.find_discrete(t0, t1, almanac.moon_nodes(eph))

2022-01-13 04:19 UTC ascending 60.6 2022-01-27 06:15 UTC descending 239.2 2022-02-09 06:12 UTC ascending 57.9 2022-02-23 06:54 UTC descending 236.2 2022-03-08 08:21 UTC ascending 54.8 2022-03-22 08:11 UTC descending 233.7 2022-04-04 13:05 UTC ascending 52.9 2022-04-18 14:01 UTC descending 232.6 2022-05-01 19:53 UTC ascending 52.5 2022-05-15 23:44 UTC descending 232.5 2022-05-29 02:34 UTC ascending 52.6 2022-09-28 23:43 UTC descending 223.8 2022-10-11 21:49 UTC ascending 43.5 2022-10-26 06:31 UTC descending 223.4 2022-11-08 06:08 UTC ascending 43.4 2022-11-22 16:23 UTC descending 223.4 2022-12-05 12:38 UTC ascending 43.2 2022-12-20 01:35 UTC descending 222.5

Making the change to step_days described above fixed this problem as well.

[brandon-rhodes]

Here are updates on these two issues:

• @bzhaas — Oh, no! That's a very serious problem, and not really related to the sunrise/sunset problem. Moon nodes don't involve edge cases like a town right on the edge of the Arctic, or a body that just barely crests the horizon for a few minutes, both of which are very hard on rising/setting routines. Instead, Moon nodes are big plain obvious great-circle events where a body riding one circle intersects another. The problem here is that I simply misunderstood how close together two nodes can fall — given that the lunar month is 29½ days, I thought I could just divide that by two, round down, and be safe. But it turns out that lunar nodes can be much closer than 14 days! As is the case with these dates that were skipped in 2022, when some nodes came only about 12.38 days apart. So I'll adjust the constant step size down to 12.0 to fix; watch for a commit mentioning this issue.
• Sunrises and sunsets are different entirely, because of the edge cases I mentioned in the previous item. A step size simply is a bad idea for that kind of search because, what if the Sun rises at your Arctic location for only half a second on a given day? So late last year I came up with a better approach and tested it out on my old PyEphem library. As soon as I get time I'm planning on writing it up for Skyfield, at which point I'll be able to close this issue. Stay tuned!

[edose]

I believe this can be solved by: (1) finding the times of altitude extrema, (2) returning None or empty list if all altitudes are above or all are below the user's horizon altitude, else (3, normal case) inserting those extrema times to the list of window edges then letting find_discrete() continue as normally.

If one relies on smaller step size alone, one will always miss some extrema. The way not to miss extrema is first to find extrema, for example, by finding changes in (sign of d(altitude)/d(time)), which is quite a smooth function, even for solar system bodies. Altitude extrema found, one cannot then miss risings and settings, however close together they are, or however polar the latitude.

[brandon-rhodes]

I've finally had time to sit down and finish the experiments I did in late 2021 that were inspired by this issue, and so the next version of Skyfield should include a rising/setting finder that is much faster, and that never misses a rising or setting.

Currently I've only documented it for the Sun:

https://github.com/skyfielders/python-skyfield/blob/1ac697dca28c0feb23b9947360d827d9e191caeb/skyfield/documentation/almanac.rst#sunrise-and-sunset

But within a few days I'll also rewrite the rising-and-settings docs, a couple of sections down, to also use this new routine. If you want to go ahead and try it out, you'll need to install the current development version of Skyfield:

pip install -U https://github.com/skyfielders/python-skyfield/archive/master.zip

If you pass something besides the Sun, try providing a 5th argument with the value -0.5667, which I'll explain when I write the docs!

[brandon-rhodes]

After much fussing with them when I've had time over the past few months, the new faster and more accurate rising and setting routines have now been released!

https://rhodesmill.org/skyfield/almanac.html#risings-and-settings

Thanks for suggesting we need routines that never miss risings or settings.