Вычисление координат планет

[Sparrow]

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

Я решил показать рассчет координат планет в некий момент времени. Сперва думал использовать ресурсы:
http://www.gregmiller.net/astro/PlanetCalc.txt
http://home.att.net/~srschmitt/zenosamples/zs_planetorbits.html
но там довольно серьезные вычисления. Мне же для демонстрации не нужно учитывать изменение кеплероских элементов через многие года. В идеале мне неоходимы несколько формул, которые по значению параметров орбиты вычисляют положение объекта в координатах (x, y, z) относительно Солнца на плюс / минус несколько лет от сегодняшнего дня. Погрешности не столь страшны, главное что бы вместе система планет выглядела более - менее реалистично.

Возможно кто-либо имеет простые алгоритмы, или может помочь мне советом / ресурсом.

Заранее благодарен.

[Tau]

Что вас так напугало? Эти вычисления занимают 20 строчек кода, если не учитывать возмущения и прецессию. 

[Sparrow]

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

[taurus]

Ничего не поделаешь: уровнение Кеплера все равно решать придется, придется перед этим время считать в юлианских столетиях и затем придется переводить координаты в XYZ. Вот, собственно, и все....

[Tau]

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

[Sparrow]

Порылся тут: http://www.gregmiller.net/astro/PlanetCalc.txt

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

     xh = r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i) )
     yh = r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i) )
     zh = r * ( sin(v+w) * sin(i) )

     r = sqrt( xv*xv + yv*yv )
     v = atan2( yv, xv )

     xv = a * ( cos(E) - e )
     yv = a * ( sqrt(1.0 - e*e) * sin(E) )

     E = M + e * sin(M) * ( 1.0 + e * cos(M) )

Значения e, M, i, N, w, a известны для планет.
Я нигде не ошибся ?
Судя по моим предположениям это получатся координаты с Солнцем в центре.

Меня больше смущает то, что данные даны ввиде:
N =  74.0005 + 1.3978E-5 * d
i = 0.7733 + 1.9E-8 * d
w =  96.6612 + 3.0565E-5 * d
a = 19.18171 - 1.55E-8 * d
e = 0.047318 + 7.45E-9 * d
M = 142.5905 + 0.011725806 * d

где
d = 367*y - 7 * ( y + (m+9)/12 ) / 4 + 275*m/9 + D - 730530

Сойдет ли для моих целей не учитывать слагаемое содержащее d, везде, кроме значения M. насколько сильные будут погрешности ?

[Tau]

Использвать первое приближение в уравнении Кеплера?

E = M + e * sin(M) * ( 1.0 + e * cos(M) )

Сойдет . Ошибки будут на уровне сотых долей а.е., то есть в несколько миллионов километров.

Только заметьте:

M = 142.5905 + 0.011725806 * d

здесь M в градусах
а здесь

E = M + e * sin(M) * ( 1.0 + e * cos(M) )

в радианах.

v, w, N, i тоже даются в градусах, а в тригонометрических функциях используются в радианах

[Sparrow]

Извините, но я не знаю что такое первое приближение

Насчет градусов и радиан я знаю. Просто пример данных был скопирован.

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

add. Значения по ссылке для планет указаны в градусах и а.е. как я понял ?)

[Tau]

Первое приближение в данном случае - это подстановка, с которой начинается решение трасцендентного уравнения Кеплера.

В каком временном промежутке вычисления будут давать адекватные значения ?

Плюс-минус 100 лет , вы же отбросили вековые члены везде, кроме M.

ЗЫ. Не перепутайте пременные e и E.

add. Значения по ссылке для планет указаны в градусах и а.е. как я понял ?)

Большая полуось A в а.е, все остальное, кроме e, в градусах. e размерности не имеет  .
XYZ будут в а.е.

[Sparrow]

Добрый день еще раз =)
А время в формуле d считается по какому календарю ?

[Tau]

По текущему, григорианскому. y- год, m - месяц, D - день.

[Tau]

Поправка: у Сатурна, Урана, Нептуна на интервале +-100 лет ошибки будут уже на уровне десятых долей а.е., так что интервал +-10 лет будет более честным.

[TMonax]

Если тема актуальна, то вот мой вариант:
класс планеты:
TPlanet = class {планеты}
  public
    C_RADIUS,
    C_MASS,
    C_SID,   // сидерический
    C_SIN,  // синодический
    C_AFEL,  // афелий
    C_PERIG,  // перигелий
    C_EX        : extended;  // эксцентриситетик)
    NamePlanet  : string;
    MaxX, MaxY,  // максимальная и минимальная полуоси эллипсоидной орбиты
    TurnTime, OrbiTime  : single; // периоды оборота вокрут своей оси и Солнца
    PrevPos,
    NowPos,
    FuturePos   : TAffineVector; // позиции относительно цены деления времени(скорости) предыдущая, текущая и буд.
  end;

все параметры планет задавал как константы и использовал при создании классов планет относительно имени указанной планеты, например:
  MARS_RADIUS = 3376.2e3; {метр}
  MARS_MASS = 6.4185e23; {килограмм}
  MARS_SID = 1.88 * EARTH_SID;
  MARS_SIN = 779.94 * 24 * 60;   {минут} {синодический - вокруг своей оси}
  MARS_AFEL = 249.23e9; {метр}
  MARS_PERIG = 206.62e9; {метр}
  MARS_EX = 0.0935;
или наприер Церера:
  CERES_RADIUS = 940000; {примерно} {метр}
  CERES_MASS = 9.5e20; {килограмм}
  CERES_SID = 4.6 * EARTH_SID; {так же относится ко всему поясу астероидов}
  CERES_SIN = 0;              {минут} {синодический - вокруг своей оси}
  CERES_AFEL = 2.85 * 149597870610; {метр}
  CERES_PERIG = 2.54 * 149597870610;
  CERES_EX = 0.08;
а так же такие параметры для координат в программе как масштаб:
  MASHTAB = 1e6; // физический масштаб

затем создавал все планеты в списке и проматывал их каждый шаг времени, при этом как планета проходила полный период вокруг себя или солнца, то я сбрасывал параметры её собственного времени(относительно него я отталкиваюсь в следующем вычислении и т.д.)
кусок кода расчёта координат:
  // движение небесных тел
  Sun.TurnAngle := MyTime.hour * 360/84960; // Солнце вокруг своей оси  
  for i:=1 to PlanetList.Count-1 do
  begin
    // сначало высчитываю личное новое время оборота вокруг себя и на орбите
    TPlanet(PlanetList.Items).TurnTime := TPlanet(PlanetList.Items).TurnTime + IncTime  -
      trunc( TPlanet(PlanetList.Items).TurnTime/TPlanet(PlanetList.Items).C_SIN ) *
      TPlanet(PlanetList.Items).C_SIN;
    TPlanet(PlanetList.Items).OrbiTime := TPlanet(PlanetList.Items).OrbiTime + IncTime -
      trunc( TPlanet(PlanetList.Items).OrbiTime/TPlanet(PlanetList.Items).C_SID ) *
      TPlanet(PlanetList.Items).C_SID;
    // а затем реально поворачиваю планету вокруг себя и ставлю в новую координату
    TPlanet(PlanetList.Items).Planet.TurnAngle :=  
      TPlanet(PlanetList.Items).TurnTime * 360/TPlanet(PlanetList.Items).C_SIN;
    // сохраняю текущее положение как прошлое
    TPlanet(PlanetList.Items).PrevPos[0] := TPlanet(PlanetList.Items).Planet.AbsolutePosition[0];
    TPlanet(PlanetList.Items).PrevPos[1] := TPlanet(PlanetList.Items).Planet.AbsolutePosition[1];
    TPlanet(PlanetList.Items).PrevPos[2] := TPlanet(PlanetList.Items).Planet.AbsolutePosition[2];
    // х
    TPlanet(PlanetList.Items).PlanetOwner.Position.X := TPlanet(PlanetList.Items).MaxX *
      Cos(GToRad( TPlanet(PlanetList.Items).OrbiTime * 360/TPlanet(PlanetList.Items).C_SID ));
    // y
    TPlanet(PlanetList.Items).PlanetOwner.Position.Z := -TPlanet(PlanetList.Items).MaxY *
      Sin(GToRad( TPlanet(PlanetList.Items).OrbiTime * 360/TPlanet(PlanetList.Items).C_SID ));
    // а z = 0 потому что я относительную систему координат планеты уже наклонил на угол её эклиптики относительно главной, т.е. она вращается по 2D эллипссу в наклоненой плоскости)
    TPlanet(PlanetList.Items).PlanetOwner.Position.Y := 0;

    TPlanet(PlanetList.Items).NowPos[0] := TPlanet(PlanetList.Items).Planet.AbsolutePosition[0];
    TPlanet(PlanetList.Items).NowPos[1] := TPlanet(PlanetList.Items).Planet.AbsolutePosition[1];
    TPlanet(PlanetList.Items).NowPos[2] := TPlanet(PlanetList.Items).Planet.AbsolutePosition[2];
    // высчитываю будущие координаты, с помощью разницы расчётов текущих и будущих, можно лучше заметить смысл формул) ибо я их сам выводил.
    TPlanet(PlanetList.Items).FuturePos[0] :=
      TPlanet(PlanetList.Items).MaxX *
      Cos(GToRad( ( TPlanet(PlanetList.Items).OrbiTime + IncTime -
      trunc( TPlanet(PlanetList.Items).OrbiTime/TPlanet(PlanetList.Items).C_SID ) *
      TPlanet(PlanetList.Items).C_SID ) * 360/TPlanet(PlanetList.Items).C_SID ));
    TPlanet(PlanetList.Items).FuturePos[1] := 0;
    TPlanet(PlanetList.Items).FuturePos[2] :=
      -TPlanet(PlanetList.Items).MaxY *
      Sin(GToRad( (TPlanet(PlanetList.Items).OrbiTime + IncTime -
      trunc( TPlanet(PlanetList.Items).OrbiTime/TPlanet(PlanetList.Items).C_SID ) *
      TPlanet(PlanetList.Items).C_SID) * 360/TPlanet(PlanetList.Items).C_SID ));
  end;

изначальное положение я находил так:
узнал что за главную космическую ось была выбрана ось пересекающая Землю в осеннее равноденствие в 2000 году, я решил расставить все планеты в эту дату в программе как в начальное положение.
на сайте наса я ввёл дату в их калькулятор и потом транспортиром вымерял углы к оси других планет, вот так и парился) а потом относительно тех положений высчитываю положения в любую дату по вышеприведённым формулам.
погрешности наверное на пару десятков тысяч километров, но не страшно для примерного отображения)
потом сверял взаимное расположение с тем же калькулятором наса - всё ровно даже в 2800 году.

вот наглядный пример работы:
http://95.129.163.190:23433/sys/sys.rar (доступно для скачивания в основном по вечерам)

PS:
т.к. сейчас делаю совершенно улучшенный вариант системы, то приведённые здесь ссылки расчётов для меня просто находка, которую я искал 3 дня - парился, поэтоу спасибо.)