infinite loop when calculating earth satellite with circular orbit

[cogenda]

I'm trying to define a satellite with user supplied perigee and apogee. I found that when the orbit is completely circular, ( _e = 0 ), pyephem is stuck in an infinite loop. On the other hand, even with a very small eccentricity (1e-15) pyephem runs fine.

The backtrace is

#0  0x00002b8604402618 in refract (pr=0, tr=15, ta=-nan(0x8000000000000), aa=0x7fffdf10d710) at libastro-3.7.5/refract.c:79
#1  0x00002b8604419ac9 in obj_earthsat (np=0x2a13680, op=0x2a136d0) at libastro-3.7.5/earthsat.c:188
#2  0x00002b8604421545 in obj_cir (np=0x2a13680, op=0xfff8000000000000) at libastro-3.7.5/circum.c:50

A test script is shown below.

import ephem

def hH_To_ne(h, H):
    '''
    @param h: Perigee (km), the minimum altitude above the Earth's surface attained
                during the spacecraft's orbit.
    @param H: Apogee (km), the maximum altitude above the Earth's surface attained
                during the spacecraft's orbit.
    '''
    # earthAvgRadius (km): The average radius of the earth.
    earthAvgRadius = 6371.2

    # a (km): Semi-major axis
    # b (km): Semi-minor axis
    a = H + earthAvgRadius
    b = h + earthAvgRadius

    # ************************
    # l : Semi-latus rectum
    # e : Eccentricity
    # al=b*b && l=a(1-e*e) => e*e = 1-b*b/a*a
    # ************************
    e = (1-(b**2)/(a**2))**(1/2.0)

    # ************************
    # n**2 = (GM/a**3)
    #   G (Nm**2kg**(-2)): The gravitational constant
    #   M (kg): The mass of central body
    #   a (m)
    #   n (radian/s)
    #
    # ==>
    #
    # _n =(approximately equal to) 8.6808 * 10**6 * a**(-3/2)
    #   a  (km)
    #   _n (turns/day)
    # ************************
    _n = 8.6808 * 10**6 * a**(-3.0/2.0)

    return _n, e

#_n, _e = hH_To_ne(499., 500.)
_n, _e = hH_To_ne(500., 500.)

sat = ephem.EarthSatellite()
sat._epoch = ephem.date("2012/9/9 06:07:00")
sat._n = _n
sat._e = _e
sat._inc = 30.
sat._raan = 0.
sat._ap = 0.
sat._M = 0.

sat.compute(ephem.date("2012/9/3 06:07:00"))

print "sat.sublat:", sat.sublat

[brandon-rhodes]

Oh, how interesting — I figured your example would be trying to compute the perigee or apogee, and complaining that the search never completes. But you are just trying to do a simple compute()? Well, drat, I had better fire up the C compiler and have a look!

[brandon-rhodes]

The error seems to occur deep in the sgp4() function in sgp4.c which is a rather old implementation of the algorithm. As you can see by reading the following report, much work has been done on tidying up the loose ends of this numerical recipe over the past few decades:

http://www.celestrak.com/publications/AIAA/2006-6753/

And I note that on Page 14 of the most recent revision of that report, it notes that one of the improvements they have made involves a test case described as, “Near Earth normal drag case but with low eccentricity (0.000 088 4) so certain drag terms are set to zero to avoid math errors / loss of precision.” This makes it sound like low eccentricity was not handled well by old versions of the algorithm, but they have fixed it recently.

It just so happens that I converted the new SGP4 algorithm to Python just in the past few weeks:

http://pypi.python.org/pypi/sgp4/

And I have near-term plans to replace the old C-language SGP4 implementation in PyEphem with this new Python version that I have released, but I am several weeks away from having this work complete. Feel free to play with the SGP4 module itself to see whether your satellite kills the algorithm like this or not. Meanwhile, I will try to continue to make progress on integrating my new SGP4 implementation into PyEphem itself, and I'll circle back around on this issue once that is working!

[brandon-rhodes]

Thanks to a recent commit, the infinite loop no longer occurs. The SGP4 error is now correctly detected by PyEphem and an exception is raised:

sat.sublat:
Traceback (most recent call last):
  File "tmp4.py", line 56, in <module>
    print "sat.sublat:", sat.sublat
RuntimeError: cannot compute the body's position at 2012/9/3 06:07:00

This improvement will be in the new version of PyEphem that I will be releasing next week.

If you want to try a more modern version of the SGP4 algorithm than the one bundled with PyEphem, I invite you to try the same calculation with Skyfield, a more recent library I have been working on — here is how its Earth satellite logic can be invoked:

http://rhodesmill.org/skyfield/earth-satellites.html