Использование pyephem для расчета, когда спутник пересекает долготу

Есть так много возможных обстоятельств, о которых могут спросить пользователи — когда спутник пересекает определенную долготу; когда он достигает определенной широты; когда он достигает определенной высоты или опускается до самой низкой высоты; когда его скорость наибольшая или наименьшая — PyEphem не пытается предоставить встроенные функции для всех из них. Вместо этого он предоставляет функцию newton(), которая позволяет вам найти переход через нуль любого сравнения, которое вы хотите провести между атрибутом спутника и заранее определенным значением этого атрибута, который вы хотите найти.

Обратите внимание, что библиотека SciPy Python содержит несколько очень тщательных функций поиска, которые намного сложнее, чем функция newton() PyEphem, на случай, если вы имеете дело с особенно плохо работающей функцией:

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

Вот как вы можете искать, когда спутник — в этом примере, МКС — проходит определенную долготу, чтобы показать общий метод. Это не самый быстрый возможный подход — в частности, поминутный поиск можно было бы ускорить, если бы мы были очень осторожны — но он написан как очень общий и очень безопасный, на случай, если есть другие значения, кроме долготы, которые вы также хотите искать. Я попытался добавить документацию и комментарии, чтобы объяснить, что происходит, и почему я использую znorm вместо того, чтобы возвращать простую разницу. Дайте мне знать, если этот сценарий работает для вас, и достаточно ясно объясняет его подход!

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

Вывод, который я получаю при запуске этого скрипта, предполагает, что он действительно нашел момент (в пределах разумного допуска), когда МКС достигает целевой долготы:

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

Существует разница во времени между временем, когда программа рассчитывает проходы по долготе, и реальным временем. Я проверил это с помощью системы LIS (она находится внутри МКС) НАСА, чтобы найти молнии. И я обнаружил, что в Европе на некоторых орбитах время, когда программа рассчитывает проход, опережает реальное время на 30 секунд. А в Колумбии на некоторых орбитах время опережения составляет около 3 минут (возможно, потому, что 1 градус долготы в Колумбии по количеству км больше, чем 1 градус долготы в Европе). Но эта проблема возникает только на двух конкретных орбитах! Тот, что проходит над Францией и уходит в Сицилию. И тот, что пролетает над США и падает на Кубу. Почему это могло быть возможно? По моему мнению, может быть, есть какая-то ошибка в алгоритме ephem.newton или, может быть, с TLE, который обычно читает созданный в 00:00:00 ночью, когда меняется день (а не фактический, потому что МКС создает 3-4 TLE в день) или, возможно, с функцией sat.sublong, которая вычисляет неправильный надир спутника. У кого-нибудь есть идея или объяснение этой проблемы? Почему так происходит ?. PS: Мне нужно проверить это наверняка, потому что мне нужно знать, когда МКС пересекает область (для обнаружения молний внутри области). И если время, которое программа вычисляет на некоторых орбитах, опережает реальное время, то функция sat.sublong вычисляет, что оно вне области (она вычисляет, что еще не прибыло в область), но программа показывает, что оно внутри область. Таким образом, в некоторых случаях реальное время не совпадает с тем, которое вычисляет программа. Большое спасибо за ваше время!