search

brandon-rhodes / pyephem

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

Comet with eccentricity close to 1.0 causing a warning #239

Closed Bernmeister closed 2 years ago

Bernmeister commented 2 years ago

When running this test code:

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

import datetime, ephem

city = ephem.city( "London" )
city.date = datetime.datetime.utcnow()
line = "C/1980 Y1 (Bradfield),e,138.5850,115.3515,358.2941,945.0557,0.0000339,0.999725,359.9999,12/27.0/1980,2000,g  9.0,4.0"
body = ephem.readdb( line )
body.compute( city )

# Prints nothing, but get the warning:
#    Near-parabolic orbit: inaccurate result.
#    e = 0.999725, lambda = 0.000138, w = 143.033560
print( "All", body.name, str( body.az ), str( body.alt ), str( body.ra ), str( body.dec ), str( body.mag ) )

# These all print.
# Same warning as above emitted by all the lines below EXCEPT for body.name.
#    Near-parabolic orbit: inaccurate result.
#    e = 0.999725, lambda = 0.000138, w = 143.033560
print( "Name", body.name )
print( "Az", str( body.az ) )
print( "Alt", str( body.alt ), str( body.ra ), str( body.dec ), str( body.mag ) )
print( "RA", str( body.ra ) )
print( "Dec", str( body.dec ) )
print( "App Mag", str( body.mag ) )

I get the following warning:

Near-parabolic orbit: inaccurate result.
e = 0.999725, lambda = 0.000138, w = 143.033560

I searched for the warning text on both the PyEphem and XEphem sites but no such luck (I also added searching .h files to no avail.

In the data line, e = 0.999725, the eccentricity, is very close to 1.0. From the conversion tool, line 83, the eccentricity should be p for parabolic, not e for elliptical, but that's a data issue, which is ultimately and, correctly I presume, tripping the warning.

Now that I know the warning exists, I can preemptively do something like:

lineSplit = line.split( ',' )
if lineSplit[ 1 ] == 'e' and float( lineSplit[ 7 ] ) < 0.99:
    ...safe to compute...

Is it possible to catch this warning?

brandon-rhodes commented 2 years ago

Searching for the words "Near-parabolic" in this repository will get you close—to within a page or two of the error message. But to search for the error message itself, ask GitHub for the phrase nNear-parabolic and it will jump right to it:

https://github.com/brandon-rhodes/pyephem/blob/7cd1754e30e4da0f2604d6f6cd05d807459690da/libastro/twobody.c#L180

Alas, I guess that GitHub search doesn't understand that \n isn't part of the words that follow!

It is indeed a bit awkward that the underlying C code is talking directly to the user rather than raising the issue where Python can see it. Let's trace: after printing the message, vrc() returns -1, which is detected in the caller obj_elliptical() and triggers a single line of code:

https://github.com/brandon-rhodes/pyephem/blob/2b695debd8fac56531750157be76c9dde0460c59/libastro/circum.c#L351

It looks like it's setting a flag I've never heard of. Let's look for its definition. Here it is:

/* mark an object as circum calculation failed */
#define NOCIRCUM    FUSER7

Well, then. I wonder if PyEphem should be checking that flag every time it calls obj_cir()? I guess I'll go add a check for it, and see if the tests still pass!

brandon-rhodes commented 2 years ago

@Bernmeister — As you'll see in the commit shown just above, I have removed the old printed warning from the underlying C library, and have instead added detection of the error flag and raised an exception for this circumstance. If you get the chance to try it out with a local install of the development version, I'll be interested to hear whether the change fixes things on your end too!

To install the development version of this project, run:

pip install -U https://github.com/brandon-rhodes/pyephem/archive/master.zip
Bernmeister commented 2 years ago

Ran the test script and got:

Traceback (most recent call last):
  File "test.py", line 18, in <module>
    print( "All", body.name, str( body.az ), str( body.alt ), str( body.ra ), str( body.dec ), str( body.mag ) )
RuntimeError: cannot compute the body's position at 2022/7/16 03:36:45 with any accuracy because its orbit is nearly parabolic and it is very far from the Sun

Winner winner, chicken dinner!

Thank you @brandon-rhodes

One point: Running sudo pip3 install -U https://github.com/brandon-rhodes/pyephem/archive/master.zip and then running the test script and gave the old/original error. It seems (and I'm guessing here) the ephem version I had installed of 4.1.3 and the version I just tested against are the same and so pip didn't actually overwrite despite downloading. So I did sudo pip3 uninstall ephem and then did the install and bingo.