Помогите найти ошибку в расчетах ПСТ по TLE

[ersar]

Здравствуйте, друзья. Обильно посыпая голову пеплом, смиренно прошу Вашей помощи.
Задумал я сделать свой проект и повторить МКС-треккер наподобие этих:
http://www.isstracker.com/
http://www.lizard-tail.com/isana/tracking/
Из документации ко второму узнал, что автор использовал модель sgp4 для расчета положения МКС. К счастью для меня, какой-то добрый человек (здоровья и счастья тебе, Brandon Rhodes) реализовал sgp4 на родном для меня Python (https://pypi.python.org/pypi/sgp4).
Достаточно быстро найдя двустрочник для МКС, я принялся за эксперименты:

Код: [Выделить]

line1 = '1 25544U 98067A   13091.91587330  .00007763  00000-0  13253-3 0  5255'
line2 = '2 25544  51.6457 126.9400 0012456  95.7479 324.5807 15.52293063822862'

year = 2013
month = 4
day = 2
hour = 11
minute = 23
second = 00

satellite = twoline2rv(line1, line2, wgs72)
position, velocity = satellite.propagate(year, month, day, hour, minute, second)
x,y,z = position

print position

Что дает мне вектор:
[854.47201141952939, 5049.4295755923677, -4473.6843863335607]

Из статьи Келсо http://space.kursknet.ru/ts_kelso/russian/v02n01/v02n01.sht я понял, что это координаты в инерциальной геоцентрической системе координат, которую можно перевести в связанную с вращением Земли  через поворот на угол звёздного времени.

Формулу расчёта звёздного времени нашёл тут: http://www2.arnes.si/~gljsentvid10/sidereal.htm

Добавлено позже: в расчёте звёздного времени, приведённом ниже, найдена ошибка (спасибо, Deimos). В коде следует исправить следующую строчку:
GMST = 281.46061837+360.98564736629*d на GMST = 280.46061837+360.98564736629*d

Код: [Выделить]

def days(year,month,day,hour,minute,second):
    year = float(year)
    month = float(month)
    day = float(day)
    hour = float(hour)
    minute = float(minute)
    second = float(second)
    dwhole = 367*year-int(7*(year+int((month+9)/12))/4)+int(275*month/9)+day-730531.5
    dfrac = (hour+minute/60+second/3600)/24
    d = dwhole+dfrac
    print d
    return d

def GMST(year,month,day,hour,minute,second):
    d = days(year,month,day,hour,minute,second)
    GMST = 281.46061837+360.98564736629*d
    GMST = GMST-360*int(GMST/360)
    if GMST<0:
        GMST = 360.0+GMST
    return GMST

print GMST(year,month,day,hours,minutes,seconds)
print GMST(year,month,day,hour,minute,second)*pi/180.0

Получаю гринвичское звёздное время в градусах: 2.71854555188
в радианах: 0.0474475707457

Поворот вокруг оси Z:

Код: [Выделить]


cz = cos(alfa) #alfa - гринвичское звёздное время в радианах
sz = sin(alfa)

Rz = matrix([
[cz,sz,0],
[-sz,cz,0],
[0,0,1]])

p =  dot(position,Rz)

После поворота вектор положения станции содержит следующие значения:
[[  614.01708475  5084.27423304 -4473.68438633]]

Для проверки результата воспользовался первым приближением и пересчитал координаты из декартовой в сферическую систему координат:

Код: [Выделить]

r = sqrt(x*x+y*y+z*z)
theta = acos(z/sqrt(x*x+y*y+z*z))
phi = atan(y/x)

print theta*180/pi,phi*180/pi

Что дало мне:
тета (широта): 131.1391
фи (долгота): 83.1138

Сверяюсь с трекерами (переведя часы на 2 апреля 2013 года 11:23:00):
Широта: 51
Долгота: 84

Вот такие пироги - помогите разобраться - что я делаю не так. Буду рад любым комментариям. Заранее спасибо!

[1212Lupus]

Что дало мне:
тета (широта): 131.1391
фи (долгота): 83.1138

Сверяюсь с трекерами (переведя часы на 2 апреля 2013 года 11:23:00):
Широта: 51
Долгота: 84

Широты отличаются почти ровно на 180^о. 

[ersar]

Спасибо за поддержку, но боюсь, что мне просто повезло с цифрой в первом разряде.

Добавлено позже: это сообщение целиком состоит из моих ошибок

Треккер на сегодняшнее утро:
5 апреля 2013 года, 7:12:00 GMT показывает:
широта: 38.608
долгота: -91.724

Моё чудо-юдо:
Координаты в инерциальной СК:
[-37270.258940999469, 19818.899315317911, 3327.9207917840313]
Звёздное время:
в углах: 193.457995957
в радианах: 3.37647899376
Координаты в вращающейся СК:
[ 40859.3484028  -10600.68958597   3327.92079178]

широта: 85.49
долгота: -14.54

[xd]

Для вектора
[[  614.01708475  5084.27423304 -4473.68438633]]
Z<0, поэтому широта меньше 0.

[-37270.258940999469, 19818.899315317911, 3327.9207917840313]
Какие единицы измерения? Вектор явно неверный, слишком далеко улетела станция.

[ersar]

Прошу прощения - случайно использовал другой TLE.

Прошу прощения повторно - сегодня 4 апреля, а не 5. Исправил это сообщение в 11:22. Чувствую, что сегодня расчётами точно заниматься не стоит

Пересчитал на то же время для правильного TLE:

Треккер на сегодняшнее утро:
4 апреля 2013 года, 7:12:00 GMT показывает:
широта: -30.5
долгота: 82.8

Моё чудо-юдо:
Координаты в инерциальной СК:
[-4501.079   164.54   5060.012]
Звёздное время:
в углах: 192.47
в радианах: 3.35
Координаты в вращающейся СК:
[ -2161.20  3951.70   5060.012] исправлено в 14:20

широта: 41.67
долгота: -61.32 исправлено в 14:20

Согласно документации к sgp4 https://pypi.python.org/pypi/sgp4 координаты в километрах от центра Земли

[xd]

Перевод выглядит правдоподобным.

Можете свериться, к примеру, с сервисом http://heavens-above.com/

[ersar]

Спасибо за ссылку. Координаты МКС по ней получаются такими же, как и на других треккерах. Информации о положении в декартовой геоцентрической системе координат я там не нашёл.

Зато благодаря часам:
http://www.heavens-above.com/whattime.aspx?lat=0&lng=0&loc=Unspecified&alt=0&tz=UCT

Смог проверить алгоритм расчёта звёздного времени и даже найти ошибку в 4 минуты.

GMST = 280.46061837+360.98564736629*d

 К сожалению, принципиально расчёт все-ещё даёт неверные результаты 

[xd]

Не допускаете ошибку при пересчёте в 1 сутки? Как раз 4 минуты набегут.

[ersar]

Ошибку в 4 минуты я локализовал и исправил - спасибо за советы. К сожалению, на общий результат это сказалось слабо. Всё еще "бью" даже не в те полушария 

Никто из форумчан не использует TLE для расчётов? Мне бы проверить - на выходе sgp4 у меня правильные данные получаются или нет.

[1212Lupus]

Лена Tau, к сожалению, разрегистрировалась -- она самый знающий в этих вопросах человек на форуме была. Ещё один спец -- Александр SleepWalker, автор программы "Heavensat" по расчёта параметров пролёта ИСЗ, может он подскажет.

[ersar]

Для интересующихся темой нашел еще одно решение проблемы: пакет pyEphem http://rhodesmill.org/pyephem/index.html

Помимо прочего с его помощью можно организовать треккинг за станцией:

Код: [Выделить]

iss = ephem.readtle('ISS',
                    '1 25544U 98067A   13095.51060533  .00016717  00000-0  10270-3 0  9008',
                    '2 25544  51.6448 109.0796 0010440 108.5171 251.7116 15.51989879 23420')
while True:
    date = datetime.utcnow()
    iss.compute(date)
    print date,iss.sublat,iss.sublong
    time.sleep(1.0)

P.S. первоначальный вопрос остаётся. Должен же мой треккер чем-то отличаться от аналогов. Будет выводить координаты в декартовой СК  Да и дело чести как бы.

[Sleepwalker]

line1 = '1 25544U 98067A   13091.91587330  .00007763  00000-0  13253-3 0  5255'
line2 = '2 25544  51.6457 126.9400 0012456  95.7479 324.5807 15.52293063822862'

Для интересующихся темой нашел еще одно решение проблемы: пакет pyEphem http://rhodesmill.org/pyephem/index.html

Помимо прочего с его помощью можно организовать треккинг за станцией:

Код: [Выделить]

iss = ephem.readtle('ISS',
                    '1 25544U 98067A   13095.51060533  .00016717  00000-0  10270-3 0  9008',
                    '2 25544  51.6448 109.0796 0010440 108.5171 251.7116 15.51989879 23420')
while True:
    date = datetime.utcnow()
    iss.compute(date)
    print date,iss.sublat,iss.sublong
    time.sleep(1.0)

P.S. первоначальный вопрос остаётся. Должен же мой треккер чем-то отличаться от аналогов. Будет выводить координаты в декартовой СК  Да и дело чести как бы.

Ошибка где-то при переводе в Гринвич, т.к. координаты В METE верные.

И еще вызывают сомнения вот эти данные

Код: [Выделить]

Треккер на сегодняшнее утро:
4 апреля 2013 года, 7:12:00 GMT показывает:
широта: -30.5
долгота: 82.8

04/04/2013 07:12:00.000 UTC

METE:
-4501.0799
164.5407
5060.0128

Гринвич:
-2443.8099
-3783.4744
5060.0038

широта:48.506
долгота:-122.859

[Scorpion_ak107]

Здравствуйте.
У меня похожая задача(только на С), тоже разбираюсь, сложности с размерностями.

Давайте разбираться вместе ?

Вот кстати вопрос - если эти координаты в км, от центра Земли, то получается станция Мир находится где-то в глубине Земли  :

author=ersar link=topic=106104.msg2368624#msg2368624 date=1365004069]Что дает мне вектор:
[854.47201141952939, 5049.4295755923677, -4473.6843863335607]

.

[xd]

Радиус-вектор 6800,05 км, что соответствует высоте 422 км, что даже высоковато будет, и уж точно не под поверхностью мирового океана.

[ersar]

04/04/2013 07:12:00.000 UTC

METE:
-4501.0799
164.5407
5060.0128

Гринвич:
-2443.8099
-3783.4744
5060.0038

Понятия не имею, как это произошло, но у меня сошлось всё вплоть до Гринвича. Над широтой и долготой еще надо поработать, но то, что уже есть воодушевляет - большое спасибо!