brandon-rhodes / pyephem

Scientific-grade astronomy routines for Python
Other
783 stars 121 forks source link

Endless loop in computing moon rise/set #232

Open j0nes2k opened 2 years ago

j0nes2k commented 2 years ago

I am using ephem==4.1.3 together with pylunar==0.6.0. When trying to compute moon rise and set times in Norway, ephem hangs in an endless loop.

Stacktrace from a test run:

self = <ephem.Observer date='2022/4/25 07:57:26' epoch='2000/1/1 12:00:00' lon='15:28:12.0' lat='78:15:00.0' elevation=0.0m horizon=-0:34:00.0 temp=15.0C pressure=0.0mBar>
body = <ephem.Moon "Moon" at 0x7fbd719d0e40>, start = None, use_center = False

    @describe_riset_search
    def next_rising(self, body, start=None, use_center=False):
        """Search for the given body's next rising"""
>       return self._find_rise_or_set(body, start, use_center, +1, True)

/workspace/virtualenvs/api_app/lib/python3.6/site-packages/ephem/__init__.py:439: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = <ephem.Observer date='2022/4/25 07:57:26' epoch='2000/1/1 12:00:00' lon='15:28:12.0' lat='78:15:00.0' elevation=0.0m horizon=-0:34:00.0 temp=15.0C pressure=0.0mBar>
body = <ephem.Moon "Moon" at 0x7fbd719d0e40>, start = None, use_center = False, direction = 1, do_rising = True

    def _find_rise_or_set(self, body, start, use_center, direction, do_rising):
        if isinstance(body, EarthSatellite):
            raise TypeError(
                'the rising and settings methods do not'
                ' support earth satellites because of their speed;'
                ' please use the higher-resolution next_pass() method'
                )

        original_pressure = self.pressure
        original_date = self.date
        try:
            self.pressure = 0.0  # otherwise geometry doesn't work
            if start is not None:
                self.date = start
            prev_ha = None
            while True:
                body.compute(self)
                horizon = self.horizon
                if not use_center:
>                   horizon -= body.radius
E                   KeyboardInterrupt
Bernmeister commented 2 years ago

If you could provide a complete test script which reproduces the issue, that would help to start diagnosing...

j0nes2k commented 2 years ago

Here is a sample script. Please note that I am using pylunar as an entry point to ephem, but the error is occuring in ephem:

from datetime import datetime
from typing import Tuple

from pylunar import MoonInfo

def dd2dms(dd: float) -> Tuple:
    """Decimal degrees to degrees minutes seconds converter"""

    is_positive = dd >= 0
    dd = abs(dd)
    minutes, seconds = divmod(dd * 3600, 60)
    degrees, minutes = divmod(minutes, 60)
    degrees = degrees if is_positive else -degrees

    return degrees, minutes, seconds

latitude = 78.2500
longitude = 15.4700
olson_id = "Arctic/Longyearbyen"
datetime_start = datetime.strptime("2022-04-26", "%Y-%m-%d")

latitude_dms = dd2dms(latitude)
longitude_dms = dd2dms(longitude)

moon_info = MoonInfo(latitude_dms, longitude_dms)
moon_info.MAIN_PHASE_CUTOFF = 12
moon_info.update(datetime_start)
moon_rise_set = moon_info.rise_set_times(olson_id)
Bernmeister commented 2 years ago

I commented out the last section of your script which referred to moon_info as I don't have pylunar installed and then ran your script without issue. I then added a couple of lines to print the latitude/longitude and got:

(78.0, 15.0, 0.0)
(15.0, 28.0, 12.0)

Your script contains no pyephem calls...

j0nes2k commented 2 years ago

Yes, using an example including pylunar makes it much easier for me to generate a reproducible test case. MoonInfo is a part of the pylunar package, pylunar calls pyephem in this method:

https://github.com/mareuter/pylunar/blob/master/pylunar/moon_info.py#L397

Are you able to install pylunar and see if the error occurs on your machine, too?

Bernmeister commented 2 years ago

Are you able to install pylunar and see if the error occurs on your machine, too?

Done using pip3. Took your script, added #!/usr/bin/env python3 to the top, ran it and stuck in a loop. Hit CTRL + C and got:

  File "test.py", line 33, in <module>
    moon_rise_set = moon_info.rise_set_times(olson_id)
  File "/usr/local/lib/python3.8/dist-packages/pylunar/moon_info.py", line 397, in rise_set_times
    mjd_time = getattr(self.observer,
  File "/usr/local/lib/python3.8/dist-packages/ephem/__init__.py", line 439, in next_rising
    return self._find_rise_or_set(body, start, use_center, +1, True)
  File "/usr/local/lib/python3.8/dist-packages/ephem/__init__.py", line 465, in _find_rise_or_set
    horizon -= body.radius

Next step I'd suggest is to isolate: is the issue in pylunar or pyephem or your code...try to call pyephem directly, bypassing pylunar.

j0nes2k commented 2 years ago

Thanks - as you can see, when you hit CTRL + C, the line where the endless loop occurs is in /usr/local/lib/python3.8/dist-packages/ephem/__init__.py - therefore this issue is related to `pyephem?

Bernmeister commented 2 years ago

Not sure if your development environment has a debugger, but that could be the next step; maybe see what happens inside of _find_rise_or_set...hopefully too @brandon-rhodes might have some more insight.

aendie commented 2 years ago

Problem verified...

A user has raised an issue regarding an endless loop using Skyalmanac, which uses Ephem for moonrise/moonset times. The same happens in Pyalmanac (which only uses Ephem). For example, this error appears in the Nautical Almanac when calculating between 23.06.2020 and 25.06.2020 inclusive.

I do not have this error with Ephem 4.1 (installed on my PC). It hangs in versions 4.1.1, 4.1.2 and 4.1.3 (installed on my laptop).

I'll try to give more specific details (with further testing).

aendie commented 2 years ago

Here is a minimal test case that loops forever when using Ephem >= 4.1.1 ....

#!/usr/bin/env python3
# -*- coding: UTF-8 -*-

import datetime
import ephem

def testcase():
    first_day = datetime.date(2020, 6, 25)  # for testing a specific date
    date = ephem.Date(first_day)      # convert date to float
    lat = 72
    print("moonrise_set {} {}".format(date,lat))

    obs = ephem.Observer()
    latitude = ephem.degrees('{}:00:00.0'.format(lat))
    obs.lat = latitude
    obs.pressure = 0
    obs.horizon = '-0:34'       # 34' (atmospheric refraction)
    m = ephem.Moon(obs)

    out  = ['--:--','--:--','--:--','--:--','--:--','--:--']    # first event

#-----------------------------------------------------------
    # Moonrise/Moonset on 1st. day ...
    d = ephem.date(date) - 30 * ephem.second
    obs.date = d
    m.compute(d)

    obs.date = d
    try:
        firstsetting = obs.next_setting(m)
        if firstsetting-obs.date >= 1:
            raise ValueError('event next day')
        out[5] = date2time(firstsetting)  # note: overflow to 00:00 next day is correct here
    except Exception:              # includes NeverUpError and AlwaysUpError
        out[5] = '--:--'

## MAIN ##
testcase()

I hope this is useful to you.

aendie commented 2 years ago

And here as a reference are correct moonrise/moonset times (hh:mm) from Ephem 4.1 for June 23 (Tue) to June 25 (Thu) 2020 from Pyalmanac for a few sample latitudes:

Moondata 2020-06-25

Change the date or latitude in my testcase and it will not loop indefinitely.

5nomads commented 2 years ago

I also ran into this issue using version 4.1.3 at high latitudes. I discovered this only happens when setting observer latitude as a string, it does not happen when I use floats or ints. It only happens with some dates, not all:

The following hangs:

def moon_rise(obs_lat, obs_datetime, ):
    body = ephem.Moon()
    observer = ephem.Observer()

    observer.date = ephem.Date(obs_datetime.strftime("%Y/%m/%d 00:00:00"))
    observer.lat = '88.5'
    observer.pressure = 0
    observer.horizon = '-0:34'
    body.compute(observer)
    try:
        return observer.next_rising(body)
    except Exception as error:
        return error

The following does not get stuck in a loop, but does NOT return the right answer:

def moon_rise(obs_lat, obs_datetime, ):
    body = ephem.Moon()
    observer = ephem.Observer()

    observer.date = ephem.Date(obs_datetime.strftime("%Y/%m/%d 00:00:00"))
    observer.lat = 88.5
    observer.pressure = 0
    observer.horizon = '-0:34'
    body.compute(observer)
    try:
        return observer.next_rising(body)
    except Exception as error:
        return error