Местоположение спутника

[DmYakov]

Благодарю Вас за оперативный ответ.

Как вы думаете, с чем может быть связано это расхождение в значениях?

[Tau]

Расхождение в сантиметрах, возможно точность где-то теряется.

[DmYakov]

Tau, здравствуйте.
Такой вопрос: в алгоритме пересчета эфемерид НКА на текущий момент времени, приведенном в ИКД ГЛОНАСС редакция 5.1, есть такое выражение:

Ek=qk+ek*sin(Ek) (стр. 59 после системы (3)).

Подскажите, верное ли это выражение? Или допущена ошибка?

[Tau]

Tau, здравствуйте.
Такой вопрос: в алгоритме пересчета эфемерид НКА на текущий момент времени, приведенном в ИКД ГЛОНАСС редакция 5.1, есть такое выражение:

Ek=qk+ek*sin(Ek) (стр. 59 после системы (3)).

Подскажите, верное ли это выражение? Или допущена ошибка?

Это так называемое уравнение Кеплера. Очень симпатичное уравнение, связывающее среднюю аномалию с эксцентрической аномалией.
E - эксцентрическая аномалия, e - эксцентриситет, q - средняя аномалия, она там ниже на странице вычисляется: q =  q0 + q1*T.
Известно q и e, нужно найти E из этого уравнения. Уравнение трансцендетное, решается итерационно. В  алгоритме без учета влияния Луны и Солнца это уравнение тоже фигурировало (стр.71).

Или вас смутило не то, что уравнение в лоб нерешаемо, а что-то другое?

[DmYakov]

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

Исходя из аналогичного примера (стр. 71) как я понял, в качестве такового (E0) нужно взять qk. Верно?

[Tau]

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

Исходя из аналогичного примера (стр. 71) как я понял, в качестве такового (E0) нужно взять qk. Верно?

Да, можно взять среднюю аномалию. Для случая малых эксцентриситетов пойдет. Только не забудьте в радианы перевести.

[DmYakov]

Tau, подскажите.
Вы реализовывали алгоритм пересчета эфемерид НКА на текущий момент времени, приведенный в ИКД ГЛОНАСС редакция 5.1?
Я бы хотел с кем-нибудь свериться, так как мои результаты расчетов не сходятся с приведенными в примере к этому алгоритму.

[Tau]

Tau, подскажите.
Вы реализовывали алгоритм пересчета эфемерид НКА на текущий момент времени, приведенный в ИКД ГЛОНАСС редакция 5.1?
Я бы хотел с кем-нибудь свериться, так как мои результаты расчетов не сходятся с приведенными в примере к этому алгоритму.

DmYakov
По этому примеру сразу вылезает ошибка: если перевести начальный вектор состояния из гринвической системы координат в ECI и разложить на орбитальные элементы, то орбита на глонассовскую не похожа совершенно. Причина в неправильных формулах перевода вектора скорости из гринвича в eci (стр 60).
X0=x*cos(S)-y*sin(S)
Y0=x*sin(S)+y*cos(S)
Это правильно. А теперь дифференцируем по времени:
Vx0=x'*cos(S)-x*sin(S)*S'-y'*sin(S)-y*cos(S)*S'=Vx*cos(S)-Vy*sin(S)-w*(x*sin(S)+y*cos(S))
Vy0=x'*sin(S)+x*cos(S)*S'+y'*cos(S)-y*sin(S)*S'=Vx*sin(S)+Vy*cos(S)+w*(x*cos(S)-y*sin(S))

[DmYakov]

Tau

Здравствуйте.

Система (1):
X0=x*cos(S)-y*sin(S), (1.1)
Y0=x*sin(S)+y*cos(S); (1.2)

Система (2):
Vx0=x'*cos(S)-x*sin(S)*S'-y'*sin(S)-y*cos(S)*S'=Vx*cos(S)-Vy*sin(S)-w*(x*sin(S)+y*cos(S))=|с учетом (1.2)|=Vx*cos(S)-Vy*sin(S)-w*Y0;
Vy0=x'*sin(S)+x*cos(S)*S'+y'*cos(S)-y*sin(S)*S'=Vx*sin(S)+Vy*cos(S)+w*(x*cos(S)-y*sin(S))=|с учетом (1.1)|=Vx*sin(S)+Vy*cos(S)+w*X0;

Вроде бы соответствует формулам в ИКД

[Tau]

Действительно, это я лоханулась, не заметила, что там X0 и Y0. Обычно эти формулы полностью через компоненты исходного вектора пишутся.
Сразу говорю, что этот геометрический ужас я не программировала, у меня тупо возмущения от Луны и Солнца вычисляются через их вектора, и расхождение с примером в ИКД в долях метров (что вообщем-то для 15-минутного интервала интегрирования не очень хорошо). Поменять возмущающую функцию в интеграторе - дело одного часа, но если у вас опять расхождения в сантиметрах и вас это не удовлетворяет, то я умываю весла .
Собственно, вопрос: какова степень несовпадения вашего результата с ИКД-шным?

[DmYakov]

 

Разность моих результатов и приведенных в примере (в метрах и метрах в секунду соответственно):

dx0(ti)=4036.635928;
dy0(ti)=14672.211381;
dz0(ti)=30659.088585;
dvx0(ti)=4.576060;
dvy0(ti)=1.458397;
dvz0(ti)=2.384721

[Tau]

Так, это действительно нужно лечить .
Сейчас пойду рыться, но предварительный вопрос такой - интегратор у вас хорошо работает? Там рекомендуется метод Рунге-Кутта 4-го порядка. Вы его запрограммировали и проверили например на случае невозмущенного движения или на любом другом дифференциальном уравнении?

[DmYakov]

Это сложный вопрос

Я, наверное, ошибся, но систему (1) я свел к следующим уравнениям:

x(k+1) = x(k) + h * vx(k),
y(k+1) = y(k) + h * vy(k),
z(k+1) = z(k) + h * vz(k),
vx(k+1) = vx(k) + h * (...),
vy(k+1) = vy(k) + h * (...),
vz(k+1) = vz(k) + h * (...),

где h - шаг, сек;
(...) - выражение, стоящее в правой части соответствующих диф. уравнений скоростей из системы (1) для момента k;

[Tau]

Я, наверное, ошибся, но систему (1) я свел к следующим уравнениям:

Могу поспорить, что этот вариант даст заметные ошибки . Но каждый должен постичь своё ДАО сам . Скачайте в сети дивную книжку "Справочное руководство по небесной механике и астродинамике" (под редакцией Г.Н. Дубошина) и начиная со страницы 477 ... 
Возвращаясь к ИКД:
Посмотрите здесь http://www.sit-com.ru/www/content/view/60/10/
Добрые люди уже выловили опечатки, их там три штуки:
1) начальная долгота перигея лунной орбиты Г0' без минуса (334 19 46.4)
2-3) при вычислении направляющих косинусов Луны  в eta11 и zeta12.
После исправления этих прелестей направляющие косинусы Луны становятся похожи на правду (с некоторой точностью). Хотя на 15-минутном интервале влияние лунно-солнечных возмущений настолько копеечное, что никак не могло привести к ошибке в 30 км.

В архиве три файла:
j2sunmoon - с учетом второй зональной гармоники (С20), Луны, Солнца.
j2 - с учетом только второй зональной гармоники.
not_pert - интегрирование уравнений невозмущенного движения. Советую начать плясать отсюда для проверки вашего интегратора.
Примечание: шаг интегрирования 1 сек, допуск решения 0.00000001 км. Метод Рунге-Кутта с переходом на метод Адамса.

Сведете ошибку к метрам - считайте, что достигли финиша. Полное совпадение с примером из ИКД все равно не получиться, потому что параметры интегрирования авторы не указали. Типа засунули в черный ящик такие числа, получили на выходе этакие.  А, может, и получите совпадение. Новичкам везет .

[DmYakov]

Tau, здравствуйте.

Я добрался-таки до эфемерид. Огромная благодарность Вам за выложенные файлы.

Реализация алгоритма у меня получилась. Я сверял результаты с Вашими и с данными из файла эфемерид ГЛОНАСС на 30-минутном отрезке.
Что заметил интересного: при использовании для расчета ускорений от лун. и сол. возм. формул из ИКД наибольшая ошибка составляла на 30-мин. интервале около 6 метров. шаг - 1 мин. (у Вас, если не ошибаюсь, в представленных файлах аналогичная ситуация - с учетом возмущений ошибка больше).

А при использовании ускорений из файла эфемерид - около 2.3 м на 30-мин. интервале. шаг - 1 мин.

Я думаю, что ошибка возникает из-за не совсем точного метода вычисления звездного времени, т. е. с погрешностью 10^-6.

В "РД 50-25645.325-89 Методические указания. Спутники Земли искусственные. Основные системы координат" приводится более точный метод с использованием данных из Астрономического ежегодника.

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

[Tau]

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

http://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html?__nnn=true
В Bulletin A даются поправки UT1-UTC.
Хотя у меня создалось впечатление, что в ИКД звездное время без этой поправки вычисляется.

[DmYakov]

Здравствуйте.
Подскажите, где использовать поправки UT1-UTC, о которых Вы говорите.

А я имел ввиду следующее: В РД 50-25645.325-89 в формулах 23 и 24 вычисляются значения нутации в долготе и нутации в наклоне с точностью 10^-6, а до этого данные нутации определяются через разложения, в которые входят различные коэффициенты, значения которых (их там, если я правильно понял, 106 штук на каждый коэффициент) приводятся, по словам из этого дока, в Астрономическом ежегоднике.

[Tau]

Подскажите, где использовать поправки UT1-UTC, о которых Вы говорите.

В РД формулы 6 и 5. Обратите внимание на переменную М формуле (5).
Можете еще почитать про алгоритмы: http://vadimchazov.narod.ru/text_pdf/comalg.pdf

А я имел ввиду следующее: В РД 50-25645.325-89 в формулах 23 и 24 вычисляются значения нутации в долготе и нутации в наклоне с точностью 10^-6, а до этого данные нутации определяются через разложения, в которые входят различные коэффициенты, значения которых (их там, если я правильно понял, 106 штук на каждый коэффициент) приводятся, по словам из этого дока, в Астрономическом ежегоднике.

В РД есть код программы на фортране, в который эти коэффициенты вбиты.
Можете их посмотреть тут (модель нутации на 1980 г):
http://earth-info.nga.mil/GandG/publications/tr8350.2/tr8350.2-a/Appendix.pdf
или тут (модель нутации на 2000 г):
http://maia.usno.navy.mil/conv2010/convupdt/chapter5/tab5.3a.txt
Но короткого разложения обычно хватает, а для него  (формулы 23 и 24) таблицы коэффициентов не нужны, самые значимые коэффициенты в формулах в явном виде присутствуют.

[ont_07]

В РД 50-25645.325-89 программа приведена с ошибками (имена переменных, формат записи чисел), впрочем прошу простить, если это только в моем варианте...
Вопрос кому удалось все правильно собрать и посчитать с точностью 6-го знака по приведенным тестам: Каким образом проверить работу модуля NUTAC? И почему после преобразования координат с использованием матрицы Rpol высота орбиты изменяется в несколько раз (300 - 2300 км, реально должно быть 800 км +/- 10 км)?? Вероятно смещения полюса не должны быть столь большими.

[Tau]

И почему после преобразования координат с использованием матрицы Rpol высота орбиты изменяется в несколько раз (300 - 2300 км, реально должно быть 800 км +/- 10 км)?? Вероятно смещения полюса не должны быть столь большими.

Матрица Rpol почти единичная (если xp и yp выражены правильно - в радианах), а почти единичная матрица вносит малые изменения. Примерчик приведите, интересно стало.