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

[Valentin]

Всем доброго времени суток!
Очень нужна помощь. Не могли бы Вы привести формулы для расчета положения нка в инерциальной системе координат по данным tle файла на заданный момент времени заданной календарной даты при невозмущенном движении. Спасибо.

[Tau]

Файл в формате TLE прочитали?
В невозмущенном движении изменяется только один параметр (из шести) - средняя аномалия (или эксцентрическая аномалия, или истинная, они взаимовыражаемы).
1. Вычисляете среднюю аномалию M на момент t:
M=M0 + 2*PI*n*(t-t0)
t-t0 - промежуток времени между нужным моментом и моментов времени из TLE (в сутках).
M0 - средняя аномалия из TLE
n - среднее движение из TLE
PI - собственно "пи"
2. Вычисляете эксцентрическую аномалию E на момент t из уравнения:
E - e*sin E = M
e - эксцентриситет из TLE
3. Вычисляете истинную аномалию v на момент t:
tg(v/2)=((1+e)/(1-e))^1/2*tg(E/2)
4. Вычисляете аргумент широты u на момент t:
u = v + w
w - аргумент перигея из TLE
5. Вычисляете геоцентрическое расстояние r до ИСЗ на момент t:
r = a*(1-e)/(1+e*cos(v))
e - эксцентриситет
a - большая полуось, вычисляетеся из среднего движения n
a=(fm/nn²)^1/3
fm = 398600.4418 гравитационная постоянная Земли
nn = n*2*PI/86400  среднее движение в радианах/cек
6. Вычисляете пространственные координаты:
x = r*(cos(u) * cos(Om) - sin(u) * sin(Om) * cos(i))
y = r*(cos(u) * sin(Om) + sin(u) * cos(Om) * cos(i))
z = r * sin(u)*sin(i)
Om - долгота восходящего узла из TLE
i - наклон орбиты из TLE

Это был пересказ азбуки.

[Valentin]

Да, файл прочитал.
Большое спасибо за подробный алгоритм.
У меня еще пару вопросов:

- можно ли для нахождения эксцентрической аномалии воспользоваться рядом
 E=M+e*sin(M)+((e^2)/2)*sin(2*M)+((e^3)/8)*(3*sin(3*M)-sin(M)   ?

-хотел бы уточнить по поводу промежутка времени t-t0. Из времени эпохи мне брать целую часть(номер дня в году) ?

[Tau]

- можно ли для нахождения эксцентрической аномалии воспользоваться рядом
 E=M+e*sin(M)+((e^2)/2)*sin(2*M)+((e^3)/8)*(3*sin(3*M)-sin(M)   ?

При больших эксцентриситетах ошибка будет заметной, если пользоваться этим разложением. Возьмите в качестве первого приближения:
E1 = M + e * Sin(M) * (1 + e * Cos(M))
и итерационно
Ek+1 = Ek - (Ek - e * Sin(Ek) - M) / (1 - e * Cos(Ek))
пока Ek+1 не перестанет изменяться.

-хотел бы уточнить по поводу промежутка времени t-t0. Из времени эпохи мне брать целую часть(номер дня в году) ?

Если у вас эпоха TLE и момент расчета  в одном году, то просто дни отнимать. И не целую часть, а вместе с долями суток. А если в разных годах, то... Лучше сделать универсально - вычислить юлианскую дату обоих моментов времени и взять их разность.

[Valentin]

Спасибо!
Вы мне очень помогли 

[Веселый гном]

Итак после небольшого перерыва я возвращаюсь к этой проблеме.
В моем примере я использую исходные данные из интерфейсного документа. Я не могу понять на каком этапе у меня ошибка. Может кто проверить?

1. Методом последовательного приближения находится большая полуось орбиты
a = 25509,8

2. Рассчитываются  момент  прохождения  восходящего  узла  орбиты  t_lamda_k  на витке, к которому принадлежит момент ti  , и долгота восходящего узла на этом витке lamda_k
t_lamda_k  = 27122,1
Omega = 6,61738

3. Вычисляются константы интегрирования на момент t_lamda_k
d_a1 = 2,13375
d_h1 = 5,89748E-8
d_l1 = 4,58295E-5
d_Omega1 = 2,07052E-7
d_i1 = 1,93652E-5
d_lamda_1 = 4,49324E-8

Тоже самое при _lamda=M+omega+n*tau;
-0,730522
-2,39878E-5
-4,88949E-5
-2,06789E-5
-6,60955E-6
2,21208E-5

4. Вычисляются поправки к элементам орбиты НКА на момент времени ti за счет влияния второй зональной гармоники С20

[Tau]

Веселый гном, всё верно.

[Веселый гном]

4 Вычисляются поправки к элементам орбиты НКА на момент времени ti за счет влияния второй зональной гармоники С20:
d_a = -2,86427
d_h = -2,40468E-5
d_l = -9,47244E-5
d_Omega = -2,0886E-5
d_i = -2,59747E-5
d_lamda2 = 2,20759E-5

5 Вычисляются возмущенные элементы орбиты НКА на момент времени t_i
Eps_i = 0,00144344
a_i = 25506,9
i_i = 1,13701
Omega_i = Omega_i
lamda_00 = 1,02173
M_i = -0,422841

6 Вычисление координат и сосстовляющих вектора скорости НКА в системе координат OXoYoZo на момент времени t_i
E_i = -0,423434
V_i = -0,424027
u_i = 1,02055
r_i = 25473,3
Vr_i = -0,00234768
Vu_i = 3,95833

Я подозреваю ошибку в E_i
E_i2 = M_i;
error2=0;
do{
    error2++;
    E_i = E_i2;
    E_i2 = M_i + Eps_i*sin(E_i);
} while ((abs(E_i - E_i2) > 1E-8)&&(error2 < 10));

[Tau]

5 Вычисляются возмущенные элементы орбиты НКА на момент времени t_i
Eps_i = 0,00144344
a_i = 25506,9
i_i = 1,13701
Omega_i = Omega_i        ?
lamda_00 = 1,02173       ?
M_i = -0,422841              ?

Правильные значения:
w_i = 1.4445749
Omega_i = 6.61736696
L* = 0.9603304
M_i  = -0.484244547

[Веселый гном]

Спасибо! Ошибку нашел. Координаты совпали. Едем далее.
Что за цифры в ИКД: 1572, 8287, 3362?

Для перевода координат в гринвичевскую систему делаю так
X = Xoi * cos(S) + Yoi * sin(S)
Y = -Xoi * sin(S) + Yoi * cos(S)
Z = Zi
где S - из пункта 2 ИКД
Получаю X, Y, Z:
17027,9
-974,274
18922,1

Далее перевожу в широту, долготу, высоту
48,0174
-3,27469
1,9108E+7

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

Может что-то с переводом времени?
10.01.2012,   05:03:30
JD = 2455936,71076 - юлианская дата
S0 = 1,90483
N0 = 10
t_i = 18210

[Tau]

Спасибо! Ошибку нашел. Координаты совпали. Едем далее.
Что за цифры в ИКД: 1572, 8287, 3362?

Это такой перенос строки. Следует читать: 10947.021572 и т.п.

Для перевода координат в гринвичевскую систему делаю так
X = Xoi * cos(S) + Yoi * sin(S)
Y = -Xoi * sin(S) + Yoi * cos(S)
Z = Zi
где S - из пункта 2 ИКД

S из п.2 - это звездное время на момент прохождения восходящего узла. А нужно звездное время на момент вычисления координат Ti.
То есть S = S0 + w*(Ti - 10800)
НО эта формула будет работать только внутри суток, в ином случае нужно к Ti прибавлять/отнимать 86400*количество суток. Это делает алгоритм заморочным. Мой вам совет: уходите от подобного вычисления звездного времени. Звездное время правильнее вычислять независимо на любой момент времени любого дня.
На момент вычисления S = 1.381556 радиан.
X= 14904.809   
Y= -8291.176   
Z= 18922.063 

Может что-то с переводом времени?
10.01.2012,   05:03:30
JD = 2455936,71076 - юлианская дата
S0 = 1,90483
N0 = 10
t_i = 18210

Не поняла, причем тут 10 января 2012? Но юлианская дата вычислена верно в случае, если 05:03:30 - всемирное время, и неверно в случае, если 05:03:30 - декретное московское время, которым оперируют в ИКД. В этом случае  JD= 2455936.58576. Всемирное время меньше декретного на 3 часа.

[m-d-k]

Здравствуйте!
Для прогноза невозмущенного движения КА даны:
дата и время эфемерид;
большая полуось;
эксцентриситет;
долгота восходящего узла;
наклонение;
аргумент перигея;
аргумент широты КА .

В литературе вместо аргумента широты КА используется время прохождения перигея.

Подскажите, пожалуйста, как расчитать время прохождения перигея из аргумента широты КА.

[Tau]

Здравствуйте!
Для прогноза невозмущенного движения КА даны:
дата и время эфемерид;
большая полуось;
эксцентриситет;
долгота восходящего узла;
наклонение;
аргумент перигея;
аргумент широты КА .

В литературе вместо аргумента широты КА используется время прохождения перигея.

Подскажите, пожалуйста, как расчитать время прохождения перигея из аргумента широты КА.

Берете любую книжку про невозмущенное движение. Вам нужны четыре формулы: зависимость средней аномалии M от времени, связь между средней аномалией М и эксцентрической аномалией  E, связь между эксцентрической аномалией E и истинной аномалией v, связь между истинной аномалией v и аргументом широты u. Если не найдете книжку - в этой теме в ответе 141 смотрите эти формулы.
Из аргумента широты, данного на некий момент времени, вычисляете истинную аномалию, эксцентрическую аномалию, среднюю аномалию на этот момент времени.
В перигее средняя аномалия равна 0 (как и все остальные).
Из формулы зависимости средней аномалии от времени вычисляете промежуток времени, который прошел от момента, на который дан аргумент широты и на который вы вычислили среднюю аномалию, до момента прохождения перигея.
Аккуратнее с размерностями, особенно у среднего движения.

[ont_07]

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

[m-d-k]

Tau:
Спасибо за подсказку!
У меня еще вопрос - для эллиптической невозмущенной орбиты эксцентрическая аномалия меняется пропорционально времени?

ont_07:
Исходные коды программ (много ссылок):
http://www.kiam1.rssi.ru/~den/links.html

[xd]

Пропорционально времени меняется средняя аномалия, которая через уравнение Кеплера приводится к эксцентрической аномалии, а затем и к истинной.
В википедии по уравнению Кеплера кстати очень неплохая статья.

[m-d-k]

Спасибо, там я нашел анимацию, которая все наглядно показала!
Как я понял, между тремя аномалиями существует чисто геометрическая зависимость.
Подскажите, пожалуйста, как аналитически выразить среднюю аномалию через истинную или эксцентрическую, зная большую полуось и эксцентриситет орбиты?

[m-d-k]

Так: E - e*sin(E)=M ?

[xd]

Астрономический календарь, постоянная часть, страница 61, формулы (1.82), (1.83).
Начинайте уже читать мануалы наконец!

[m-d-k]

Спасибо за ссылку на первоисточник - там такая же формула!
Получилась такая цепочка:
1. Считаю разность юлианских дат эфемерид и прогноза.
2. Определяю время прогноза t в звездных секундах - делю разность на (86400*0,997269566329).
3. По эфемеридам и временному интервалу t определяю положение спутника в плоскости орбиты.
4. Уточняю положение восходящего узла и аргумента перигея с учетом вековых прецессий.
5. Определяю широту и долготу спутника относительно линии весны.
6. Определяю звездное время на дату прогноза (угол поворота Земли) и географическую долготу спутника.
7. Преобразую геоцентрическую широту в геодезическую (для отображения на карте).
Подскажите, пожалуйста, такая последовательность правильная?