Телескопы покупают здесь


A A A A Автор Тема: Определение элементов орбиты  (Прочитано 1153 раз)

0 Пользователей и 1 Гость просматривают эту тему.

Оффлайн LilimАвтор темы

  • Новичок
  • *
  • Сообщений: 9
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Lilim
Определение элементов орбиты
« : 24 Сен 2018 [13:04:00] »
Доброго дня!

Дано: r(x,y,z) и v(x,y,z) - радиус вектор и скорость искусственного спутника в текущий момент времени. Также известны все параметры центрального тела.
Нужно: "прикинуть" орбиту ИС, т.е. найти её элементы. Для упрощения - расчет будет производится только в случае эллиптической орбиты.

Меня хватило только на определение нормали плоскости орбиты через скалярное произведение нормированных r и v.

Вы мне не поможите с расчетами? Хотелось бы получить пошаговую инструкцию с формулами для определения эксцентриситета, величин большой и малой осей эллипса ну и всего остального что можно найти.

Оффлайн xd

  • *****
  • Сообщений: 17 977
  • Благодарностей: 378
    • Skype - deimos.belastro.net
  • Награды Открытие комет, астероидов, сверхновых звезд, научно значимые исследования.
    • Сообщения от xd
    • Белорусская любительская астрономическая сеть
Re: Определение элементов орбиты
« Ответ #1 : 24 Сен 2018 [13:15:50] »
Астрономический календарь, постоянная часть. Глава I, параграф 20 "Определение орбит".
У природы нет плохой погоды, у неё просто на нас аллергия.

Учение без размышления бесполезно, но и размышление без учения опасно /Конфуций/
Слово есть поступок. /Л. Толстой/

Оффлайн Astro_Coder

  • Новичок
  • *
  • Сообщений: 31
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Astro_Coder
Re: Определение элементов орбиты
« Ответ #2 : 24 Сен 2018 [14:30:22] »
Еще в этой есть - Монтенбрук, Пфлегер "Астрономия на персональном компьютере".

Оффлайн xd

  • *****
  • Сообщений: 17 977
  • Благодарностей: 378
    • Skype - deimos.belastro.net
  • Награды Открытие комет, астероидов, сверхновых звезд, научно значимые исследования.
    • Сообщения от xd
    • Белорусская любительская астрономическая сеть
Re: Определение элементов орбиты
« Ответ #3 : 24 Сен 2018 [15:22:33] »
Кстати, если для расчёта использовать некоторый произвольный радиус-вектор как единичный, а первую космическую скорость на этом единичном расстоянии как единичную скорость, и пронормировать радиус и скорость по ним, то уравнения примут очень простой вид.

Исходные данные: вектор положения \(\vec{r}\), вектор скорости \( \vec{v} \), момент времени \(t\).

В этих же нормированных величинах большая полуось: \(\displaystyle a = \frac{1}{2r-v^2} \)
Период обращения спутника \(\displaystyle T = 2\pi a^{3/2} \)
Вводим вспомогательный параметр \(\displaystyle d = \vec{r}\cdot \vec{v} \) (скалярное произведение)
Эксцентриситет и эксцентрическая аномалия считаются из пары уравнений:

\(\displaystyle  e \sin E = \frac{d}{\sqrt{a}} \)
\(\displaystyle  e cos E = 1 - \frac{r}{a} \)

Средняя аномалия из уравнения Кеплера: \(\displaystyle M = E - e \sin E \)

Дальше с помощью дополнительных вспомогательных параметров \(\displaystyle  q = r - \frac{d^2}{2} \) и \(\displaystyle  \theta = \frac{d}{\sqrt{2q}} \) находим момент прохождения перицентра \(\tau\):

\(\displaystyle  \frac{t - \tau}{T} = \sqrt{\frac{q^3}{2\pi^2}}\left( \theta + \frac{\theta^3}{3} \right) \)

Долгота восходящего узла, наклонение орбиты и перицентрическое расстояние выражается матричным уравнением:

\( p \left( \sin i \sin \Omega ; \sin i \cos \Omega ; \cos i \right)^T = \vec{r}\times\vec{v} \) (справа векторное произведение)

Истинная аномалия считается так:

\( \kappa \cos \nu = p - r \)
\( \kappa \sin \nu = d \sqrt{r} \)

(величина \(\kappa\) не нужна )

Вспомогательная величина \(\upsilon \)

\( \rho \cos \upsilon = (x \cos \Omega + y \sin \Omega) \sin i \)
\( \rho \sin \upsilon = z \)

Здесь x,y,z - компоненты вектора \(\vec{r}\)

Долгота восходящего узла:

\(\omega = \upsilon - \nu \)
У природы нет плохой погоды, у неё просто на нас аллергия.

Учение без размышления бесполезно, но и размышление без учения опасно /Конфуций/
Слово есть поступок. /Л. Толстой/

Онлайн Toth

  • *****
  • Сообщений: 2 581
  • Благодарностей: 174
    • Сообщения от Toth
Re: Определение элементов орбиты
« Ответ #4 : 24 Сен 2018 [21:41:48] »
Вот, если надо , код на Паскаль. Проверял, работает при e<1. Правда, если орбита вдруг окажется гипер/параболической, будет переполнение.
unit modDetOr;

interface
uses Math;
// в систеие СИ, углы - радианы
const
 GM_EA=3.9860043289694E+14; 

type tVector=record
  X,Y,Z:Double;
end;
procedure deterOrbit(vP,vV:tVector;var a,e,i,node,w,M:Double);  // e<1
implementation
//========
function modVector(v:tVector):Double;  // = |v|
begin
 Result:=Sqrt(v.X*v.X+v.Y*v.Y+v.Z*v.Z);
end;
function proScal(v1,v2:tVector):Double; // =v1*v2  скалярное
begin
 Result:=v1.x*v2.x+v1.y*v2.y+v1.z*v2.z;
end;
function to2PI(x:Double):Double;  // угол x -> 0..2пи
begin
 if x<0
  then Result:=x+2*PI
  else if x>2*PI
    then Result:=x-2*PI
    else Result:=x;
end;
procedure deterOrbit(vP,vV:tVector;var a,e,i,node,w,M:Double);
var v0,r0,v2c,v0v2c,cosFi,e2,v_p,p,nu,EE:Double;
begin
  v0:=modVector(vV);
  r0:=modVector(vP);
  v2c:=GM_EA/r0;
  v0v2c:=v0*v0/v2c;
  a:=r0/(2-v0v2c);
  v_p:=proScal(vP,vV);
  cosFi:=v_p/(r0*v0);
  e2:=(v0v2c-1)*(v0v2c-1)+(r0/a)*v0v2c*cosFi*cosFi;
  e:=Sqrt(e2);
  p:=a*(1-e2);
  nu:=ArcTan2(v_p*Sqrt(p/GM_EA),p-r0);
  i:=ArcCos((vP.X*vV.Y-vP.Y*vV.X)/Sqrt(p*GM_EA));
  node:=to2PI(ArcTan2(vP.Y*vV.Z-vP.Z*vV.Y,vP.X*vV.Z-vP.Z*vV.X));
  w:=to2PI(ArcTan2(vP.Z*Cosecant(i),vP.X*Cos(node)+vP.Y*Sin(node))-nu);
  EE:=2*ArcTan(Tan(nu/2)*Sqrt((1-e)/(1+e)));
  M:=to2PI(EE-e*Sin(EE));
end;
end.

Если надо, могу и для e>1 сделать.

Онлайн Toth

  • *****
  • Сообщений: 2 581
  • Благодарностей: 174
    • Сообщения от Toth
Re: Определение элементов орбиты
« Ответ #5 : 24 Сен 2018 [22:05:28] »
Не стал исправлять прошлый, чтобы не запутались.
Сделал для элл. и гипер. орбит. Вот сама процедура . Если e>1 , то a получится < 0.
procedure deterOrbitAll(vK,vV:tVector;var a,e,i,node,w,M:Double);
var r0,V0,r0V0,t,E0,pp,eSinE:Double;
begin
 r0:=ModVector(vK);
 V0:=ModVector(vV);
 t:=2/r0-V0*V0/GM_EA;
 a:=1/t;
 r0V0:=proScal(vK,vV);
 eSinE:=r0V0/(Sqrt(GM_EA*Abs(a)));
 t:=1-r0/a;
 if a>0 then
  begin    // эллипс
   E0:=ArcTan2(eSinE,t);
   e:=t/Cos(E0);
   M:=E0-e*Sin(E0);
   if M<0 then M:=M+2*PI;
  end   else
  begin    // гипербола
   E0:=ArcTanH(eSinE/t);
   e:=t/CosH(E0);
   M:=-E0+e*SinH(E0);
  end;
 pp:=a*(1-e*e);
 t:=(vK.X*vV.Y-vK.Y*vV.X)/Sqrt(pp*GM_EA);
 i:=ArcCos(t);
 Node:=ArcTan2(vK.Y*vV.Z-vK.Z*vV.Y,vK.X*vV.Z-vK.Z*vV.X);
 if Node<0 then Node:=Node+2*PI;
 w:=ArcTan2(vK.Z/Sin(i),vK.X*Cos(Node)+vK.Y*Sin(Node))-ArcTan2(Sqrt(pp)*r0V0,Sqrt(GM_EA)*(pp-r0));
 if w<0 then w:=w+2*PI;
 if w>(2*PI) then w:=w-2*PI;
end;

Оффлайн xd

  • *****
  • Сообщений: 17 977
  • Благодарностей: 378
    • Skype - deimos.belastro.net
  • Награды Открытие комет, астероидов, сверхновых звезд, научно значимые исследования.
    • Сообщения от xd
    • Белорусская любительская астрономическая сеть
Re: Определение элементов орбиты
« Ответ #6 : 24 Сен 2018 [22:18:14] »
Паскаль... Кровь из глаз :(
У природы нет плохой погоды, у неё просто на нас аллергия.

Учение без размышления бесполезно, но и размышление без учения опасно /Конфуций/
Слово есть поступок. /Л. Толстой/

Оффлайн LilimАвтор темы

  • Новичок
  • *
  • Сообщений: 9
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Lilim
Re: Определение элементов орбиты
« Ответ #7 : 25 Сен 2018 [12:11:15] »
Всем спасибо! Особо благодарю товарища Toth за практически готовый код. Я лишь перевел его на С++.
// a - большая полуось
// e - эксцентриситет
// i - наклонение орбиты
// node - долгота восходящего узла
// w - аргумент перицентра
// M - средняя аномалия
void allObdit(const vec3d &P, const vec3d &V, double GM, double &a, double &e, double &i, double &node, double &w, double &M)
{
static const double _2PI = 6.283185307179586476925286766559;

double v0 = V.length();
double r0 = P.length();
double t = 2 / r0 - v0*v0 / GM;
a = 1 / t;

double r0v0 = P * V;
double eSinE = r0v0 / sqrt(GM*abs(a));
t = 1 - r0 / a;

double E0;
if (a > 0) // эллипс
{
E0 = atan2(eSinE, t);
e = t / cos(E0);
M = E0 - e*sin(E0);
if (M < 0) M = M + _2PI;
}
else      // гипербола
{
E0 = atanh(eSinE / t);
e = t / cosh(E0);
M = -E0 + e*sinh(E0);
}

double pp = a*(1 - e*e);
t = (P.x()*V.y() - P.y()*V.x()) / sqrt(pp*GM);
i = acos(t);

node = atan2(P.y()*V.z() - P.z()*V.y(), P.x()*V.z() - P.z()*V.x());
if (node<0) node = node + _2PI;

w = atan2(P.z() / sin(i), P.x()*cos(node) + P.y()*sin(node)) - atan2(sqrt(pp)*r0v0, sqrt(GM)*(pp - r0));
if (w<0) w = w + _2PI;
if (w>_2PI) w = w - _2PI;
}


Онлайн Ulmo

  • *****
  • Сообщений: 1 887
  • Благодарностей: 67
    • Сообщения от Ulmo
Re: Определение элементов орбиты
« Ответ #8 : 25 Сен 2018 [15:21:35] »
Паскаль... Кровь из глаз :(
Паскаль отличный язык. На нем же код читается без ломания мозга. Да и написано на нем много чего полезного.
Я до сих пор с содроганием впоминаю конструкции С++ вида а+++b :)
Вот, если надо , код
Спасибо!

Оффлайн xd

  • *****
  • Сообщений: 17 977
  • Благодарностей: 378
    • Skype - deimos.belastro.net
  • Награды Открытие комет, астероидов, сверхновых звезд, научно значимые исследования.
    • Сообщения от xd
    • Белорусская любительская астрономическая сеть
Re: Определение элементов орбиты
« Ответ #9 : 25 Сен 2018 [15:29:19] »
Я до сих пор с содроганием впоминаю конструкции С++ вида а+++b :)
За такие конструкции в приличном обществе руки обрывают по самые шорты ;D
У природы нет плохой погоды, у неё просто на нас аллергия.

Учение без размышления бесполезно, но и размышление без учения опасно /Конфуций/
Слово есть поступок. /Л. Толстой/

Оффлайн XiToi

  • Новичок
  • *
  • Сообщений: 18
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от XiToi
Re: Определение элементов орбиты
« Ответ #10 : 14 Окт 2018 [00:19:45] »
Это будет работать если, объект будет перемещаться по орбите и вектор скорости так же будет меняться?

Онлайн Toth

  • *****
  • Сообщений: 2 581
  • Благодарностей: 174
    • Сообщения от Toth
Re: Определение элементов орбиты
« Ответ #11 : 14 Окт 2018 [00:27:00] »
Это будет работать если, объект будет перемещаться по орбите и вектор скорости так же будет меняться?
Конечно. А как иначе? В задаче 2-х тел не бывает решений  V=const .

Оффлайн LilimАвтор темы

  • Новичок
  • *
  • Сообщений: 9
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Lilim
Re: Определение элементов орбиты
« Ответ #12 : 30 Окт 2018 [16:23:02] »
В опытах наблюдается несоответствие траекторий.


procedure deterOrbitAll(vK,vV:tVector;var a,e,i,node,w,M:Double);

Я так понемаю vK это вектор начального положения. Это вектор положения относительно центра притяжения или центра эллипса? Просто у меня получается, по всей видимости, что относительно центра эллипса.

Онлайн Toth

  • *****
  • Сообщений: 2 581
  • Благодарностей: 174
    • Сообщения от Toth
Re: Определение элементов орбиты
« Ответ #13 : 30 Окт 2018 [16:46:16] »
Это вектор положения относительно центра притяжения или центра эллипса?
Относительно центра Земли (или иного центрального тела ), то есть фокуса эллипса ( гиперболы ).

Оффлайн LilimАвтор темы

  • Новичок
  • *
  • Сообщений: 9
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Lilim
Re: Определение элементов орбиты
« Ответ #14 : 30 Окт 2018 [17:26:35] »
Давайте разбираться.
Вот исходные данные задачи -
    Декартовые координаты (у - вверх, x - вправо), в центре тело тяготения массой 5e10;
    GM = 6.672e-11 * 5e10 = 3,336;
    Pos(0,-10,0);
    Vel(0.7 , 0, 0);

В симуляторе, каждую секунду я высчитываю вектор гравитации, изменяю скорость, добовляя к ней вектор гравитации и полученный вектор перемещения добавляю к текущиму положению тела. Вот траектория полученная таким путем:


https://drive.google.com/open?id=1P9xzI-qx2Pf_rFrdmJKgJo2XM4GWGaL9

Оффлайн LilimАвтор темы

  • Новичок
  • *
  • Сообщений: 9
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Lilim
Re: Определение элементов орбиты
« Ответ #15 : 30 Окт 2018 [17:33:25] »
А здесь прогнозируемая deterOrbitAll орбита в техже условиях:

https://drive.google.com/open?id=1SuIC40Lzb04VghHThnXFPLVvuZMOdPVG

a = 18.826185101580123
e = 0.46882494004796127
i = 0
node = -0
w = NaN
M = 0

Оффлайн LilimАвтор темы

  • Новичок
  • *
  • Сообщений: 9
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Lilim
Re: Определение элементов орбиты
« Ответ #16 : 30 Окт 2018 [17:55:20] »
Toth не мог бы ты прокомментировать все действия в функции deterOrbitAll

Онлайн Toth

  • *****
  • Сообщений: 2 581
  • Благодарностей: 174
    • Сообщения от Toth
Re: Определение элементов орбиты
« Ответ #17 : 30 Окт 2018 [17:55:53] »
У меня небольшой глюк при нулевых наклонениях, поэтому я взял скорость (0.7;0;0.0000001)
У меня получилось при замене GM_EA на 3.336:
a := 18.8261851015812
e := 0.468824940047991
i := 8E-6
Node := 270
w := 0
M := 0  то есть начинается с перицентра

w+node = долгота перицентра = 270 -то есть перицентр внизу.
Вроде все правильно.

PS Надо что-то исправить , чтоб при нулевых наклонениях считало.

Оффлайн LilimАвтор темы

  • Новичок
  • *
  • Сообщений: 9
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Lilim
Re: Определение элементов орбиты
« Ответ #18 : 30 Окт 2018 [18:06:02] »
Наверное я неправильно нахожу малую полуось.

b = 1 - a * e

Оффлайн LilimАвтор темы

  • Новичок
  • *
  • Сообщений: 9
  • Благодарностей: 0
  • Мне нравится этот форум!
    • Сообщения от Lilim
Re: Определение элементов орбиты
« Ответ #19 : 30 Окт 2018 [18:17:46] »
Так и есть, нужно:

b = Math.sqrt(1 - orbit.e * orbit.e) * a;