separation() returns huge angle

Question

I am writing a program which goes through FITS files with photometry and looks for stars given in a .dat file. One of the steps is computing distances between two given stars using ephem.separation()

It works well. However, from time to time separation returns angles like 1389660529:33:00.8

import ephem
import math

star = ['21:45:15.00', '65:49:24.0']
first_coo = ['21:45:15.00', '65:49:24.0']

check = ephem.FixedBody()
check._ra = ephem.hours(star[0])
check._dec = ephem.degrees(star[1])
check.compute()

# star is a list with coordinates, strings in form %s:%s:%s

first = ephem.FixedBody()
first._ra = ephem.hours(first_coo[0])
first._dec = ephem.degrees(first_coo[1])
first.compute()

sep = math.degrees(float(ephem.separation(check,first)))
print sep

It occurs randomly. Have anybody encountered such behaviour?

I search for 18 stars in 212 files, which makes 3816 cycles. Might have something to do with it?

Answer 1

UPDATE: I have released a new PyEphem 3.7.5.2 that fixes this special case of comparing an angle to itself.

Your code sample has two interesting features: first, it contains a slight bug that I thought at first might be behind the problem; and, second, I was wrong that your code was the problem because your code does indeed expose a flaw in the separation() function when it is asked how far a position is from itself!

The bug in your own code is that calling compute() and asking about .ra and .dec returns those coordinates in the coordinate system of the very moment that you are calling compute() for — so your two compute() calls are returning coordinates in two different coordinate systems that are very slightly different, and so the resulting positions cannot be meaningfully compared with separation() because separation() requires two coordinates that are in the same coordinate system.

To fix this problem, chose a single now moment to use as your equinox epoch:

now = ephem.now()
...
check.compute(epoch=now)
...
first.compute(epoch=now)

That will give you coordinates that can be meaningfully compared.

Now, on to the problem in PyEphem itself!

When presented with two copies of the same position are provided to separation() it goes ahead and tries to find a distance between them anyway, and winds up doing a calculation that amounts to:

acos(sin(angle) ** 2.0 + cos(angle) ** 2.0)

which should remind us of the standard Euclidean distance formula but with an acos() around it instead of a sqrt(). The problem is that while in theory the value inside of the parens should always be 1.0, the rounding inside of IEEE floating point math does not always produce squares that sum to exactly 1.0. Instead, it sometimes produces a value that is a bit high or a bit low.

When the result is a bit below 1.0, separation() will return a very slight separation for the coordinates, even though they are actually “the same coordinate.”

When the result exceeds 1.0 by a little bit, separation() will return not-a-number (nan) on my system — and, I will bet, returns that huge return value that you are seeing printed out on yours — because cos() cannot, by definition, return a number greater than 1.0, so there is no answer acos() can return when asked “what angle returns this value that is greater than one?”

I have created a bug report in GitHub and will fix this for the next version of PyEphem:

https://github.com/brandon-rhodes/pyephem/issues/31

Meanwhile, your code will have to avoid calling separation() for two positions that are actually the same position — could you use an if statement with two == comparisons between ra and dec to detect such cases and just use the value 0.0 for the separation instead?