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?

Kiushu

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

Embank

Sorry, if I understand, in the time "2013/04/20 00:18:21" the satellite Altitude from the longitude -83.8889 is 90º?, I mean, is "exactly" above of the longitude -83.8889? – Kikelia

The satellite is somewhere above that longitude. Of course, that could be anywhere from the north pole to the south, and only for one point along that line of longitude will the satellite be straight up. For the rest of the sites along that line of longitude the satellite will be lower than 90° or even below the horizon. – Embank

Answer 2

There is a difference of time between the time the program calculates the passes over the longitude and the real time. I've checked it with the LIS system ( that it's inside the ISS ) of the Nasa to find lightnings. And I have discovered that in Europe in some orbits the time that the program calculates the pass, it's 30 seconds in advanced than the real time. And in Colombia in some orbits the time in advanced is about 3 minutes ( perhaps because 1 degree of longitud in Colombia is bigger in amount of Km than 1 degree of longitude in Europe ). But this problem only happens in 2 particular orbits ! The one that pass over France and goes down in Sicilia. And the one that pass over USA, and goes down in Cuba. Why could this be possible ? In my pinion I think maybe there's some mistake in the ephem.newton algorithm or maybe with the TLE, that normally it reads the one created at 00:00:00 at night when it changes of day (and not the actual, because the ISS creates 3-4 TLE per day ) or maybe with the sat.sublong function that calculates a wrong Nadir of the satellite. Does anyone has an idea or an explanation for this problem ? Why it happens ?. PS: I need to checked it for sure because I need to know when the ISS crosses an area ( for detecting the lightnings inside the area ). And if the time that the program calculates in some orbits it's in advanced than the real time, then the sat.sublong function calculates it's outside the area ( it calculates it hasn't arrived to the area yet ) but the program shows it's inside the area. So the real time doesn't match with the one that the program calculates, in some occasions. Thanks a lot for your time !

Aftermost