Using pyephem to calculate when a satellite crosses a Longitude

Question

I am having a hard time figuring out how to calculate when a satellite crosses a specific Longitude. It would be nice to able to provide a time period and a TLE and be able to return all the times at which the satellite crosses a given longitude during the specified time period. Does pyephem support something like this?

Answer 1

There are so many possible circumstances that users might ask about — when a satellite crosses a specific longitude; when it reaches a specific latitude; when it reaches a certain height or descends to its lowest altitude; when its velocity is greatest or least — that PyEphem does not try to provide built-in functions for all of them. Instead, it provides a newton() function that lets you find the zero-crossing of whatever comparison you want to make between a satellite attribute and a pre-determined value of that attribute that you want to search for.

Note that the SciPy Python library contains several very careful search functions that are much more sophisticated than PyEphem's newton() function, in case you are dealing with a particularly poorly-behaved function:

http://docs.scipy.org/doc/scipy/reference/optimize.html

Here is how you might search for when a satellite — in this example, the ISS — passes a particular longitude, to show the general technique. This is not the fastest possible approach — the minute-by-minute search, in particular, could be sped up if we were very careful — but it is written to be very general and very safe, in case there are other values besides longitude that you also want to search for. I have tried to add documentation and comments to explain what is going on, and why I use znorm instead of returning the simple difference. Let me know if this script works for you, and explains its approach clearly enough!

import ephem

line0 = 'ISS (ZARYA)             '
line1 = '1 25544U 98067A   13110.27262069  .00008419  00000-0  14271-3 0  6447'
line2 = '2 25544  51.6474  35.7007 0010356 160.4171 304.1803 15.52381363825715'

sat = ephem.readtle(line0, line1, line2)
target_long = ephem.degrees('-83.8889')

def longitude_difference(t):
    '''Return how far the satellite is from the target longitude.

    Note carefully that this function does not simply return the
    difference of the two longitudes, since that would produce a
    terrible jagged discontinuity from 2pi to 0 when the satellite
    crosses from -180 to 180 degrees longitude, which could happen to be
    a point close to the target longitude.  So after computing the
    difference in the two angles we run degrees.znorm on it, so that the
    result is smooth around the point of zero difference, and the
    discontinuity sits as far away from the target position as possible.

    '''
    sat.compute(t)
    return ephem.degrees(sat.sublong - target_long).znorm

t = ephem.date('2013/4/20')

# How did I know to make jumps by minute here?  I experimented: a
# `print` statement in the loop showing the difference showed huge jumps
# when looping by a day or hour at a time, but minute-by-minute results
# were small enough steps to bring the satellite gradually closer to the
# target longitude at a rate slow enough that we could stop near it.
#
# The direction that the ISS travels makes the longitude difference
# increase with time; `print` statements at one-minute increments show a
# series like this:
#
# -25:16:40.9
# -19:47:17.3
# -14:03:34.0
# -8:09:21.0
# -2:09:27.0
# 3:50:44.9
# 9:45:50.0
# 15:30:54.7
#
# So the first `while` loop detects if we are in the rising, positive
# region of this negative-positive pattern and skips the positive
# region, since if the difference is positive then the ISS has already
# passed the target longitude and is on its way around the rest of
# the planet.

d = longitude_difference(t)

while d > 0:
    t += ephem.minute
    sat.compute(t)
    d = longitude_difference(t)

# We now know that we are on the negative-valued portion of the cycle,
# and that the ISS is closing in on our longitude.  So we keep going
# only as long as the difference is negative, since once it jumps to
# positive the ISS has passed the target longitude, as in the sample
# data series above when the difference goes from -2:09:27.0 to
# 3:50:44.9.

while d < 0:
    t += ephem.minute
    sat.compute(t)
    d = longitude_difference(t)

# We are now sitting at a point in time when the ISS has just passed the
# target longitude.  The znorm of the longitude difference ought to be a
# gently sloping zero-crossing curve in this region, so it should be
# safe to set Newton's method to work on it!

tn = ephem.newton(longitude_difference, t - ephem.minute, t)

# This should be the answer!  So we print it, and also double-check
# ourselves by printing the longitude to see how closely it matches.

print 'When did ISS cross this longitude?', target_long
print 'At this specific date and time:', ephem.date(tn)

sat.compute(tn)

print 'To double-check, at that time, sublong =', sat.sublong

The output that I get when running this script suggests that it has indeed found the moment (within reasonable tolerance) when the ISS reaches the target longitude:

When did ISS cross this longitude? -83:53:20.0
At this specific date and time: 2013/4/20 00:18:21
To double-check, at that time, sublong = -83:53:20.1