UProLa

Неокріпші думки

Мандруючи Всесвітом. Розрахунок орбіти планети

with 5 comments

Пролог.

Коперник вважав, що орбіта планети – ідельне коло. Якби так і було насправді, то і проблем з позиціюванням планети в просторі навколо Сонця не було б(шкільна тригонометрія). Але з”явився Кеплер… Він сформулював три закони планетарних орбіт, ввів параметри орбіти і він просто крутий чел. Зараз спробую популярно пояснити розрахунок орбіти по Кеплеру і Ньютону.

Система координат

  1. Будь-яка планета рухається по орбіті навколо Зірки (центру мас кількох Зірок). Орбіта – еліпс, в одному з його фокусів знаходиться Зірка. В зв”язку з тим, що Зірка нерухома (в нашому випадку), можна вважати її центр (центр мас) як початок усіх координат. Уявіть собі точку в чорному просторі – це буде початок координат.
  2. Від початку координат відходять три ортогональні промені – вектори.
    Вектор, спрямований вверх – Y. Вправо – Х. До нас – Z.
    Екліптика (в нашому випадку) – площина XZ.
    Початковий вектор – вектор X. Уявіть собі ці вектори і площини, оскільки від них буде все відштовхуватись.
  3. Орієнтація орбіти в просторі задається трьома кутамиi, \Omega, \omega.
    i – нахил орбіти до площини екліптики. Напрямок – проти годинникової стрілки в нашій системі(ми дивимось на XY, вектор Z прямує до нас).
    \Omega – орієнтація орбіти відносно вектора X. Точніше, кут повороту лінії перетину площини орбіти та екліптики відносно вектора X в площині XZ, проти годинниковою стрілкою (ми дивимось на XZ, вектор Y прямує до нас). Цими кутами задається площина орбіти.
    \omega – кут повороту орбіти в площині орбіти. Відлік йде від лінії перетину екліптики та площини орбіти до великої півосі еліпса орбіти, яка проходить через центр мас (Зірку). Відлік – проти годинникової стрілки (якщо лінія перетину спрямована вліво від нас, а велика піввісь – вверх).
    Ці три параметри точно задають положення орбіти в просторі відносно центру мас (Зірки).

Слова тяжко дають зрозуміти картину, тому покажу картинку з вікіпедії.
Положення орбіти

Параметри орбіти.

Оскільки положення уже визначене, приступимо до конкретно розрахунку орбіти. Наше завдання – знайти функції x(t), y(t), де t – час, який пройшов від початку руху планети з найближчої до центру мас точки до теперішнього її положення. Маючи координати, задані параметрично, ми зможемо знайти всі потрібні нам величини.

  1. Параметри орбіти:
    a – довжина великої півосі еліпса
    \varepsilon – ексцентриситет еліпса
    Ms – маса Зірки (Зірок навколо центру мас)
    m – маса планети
  2. Константи
    G – гравітаційна константа
    \pi – відношення довжини кола до його діаметру

    Ці чотири змінні повністю детермінують орбіту тому є основними в розрахунках. Всі інші параметри та змінні будуть допоміжними.

  3. Знайдемо період обертання планети на орбіти з модифікованого Ньютоном Третього закону Кеплера:
    P=\sqrt\frac{4\pi^2a^3}{G(Ms+m)}
  4. Знайдемо середню аномалію (характеризує положення планети на орбіті):
    M=\frac{2\pi}{P}t
  5. Розв”яжемо трансцендентне рівняння відносно E:
    M=E-\varepsilon\sin E
  6. Обчислимо істинну аномалію \theta (фактично є кутом між великою піввіссю еліпса і прямою, яка з”єднує центр мас і планету, True Anomaly на малюнку):
    \theta=2arctg(\sqrt\frac{1+\varepsilon}{1-\varepsilon} tg\frac{E}{2})
  7. Обчислимо відстань з центру мас до планети \rho:
    \rho=\frac{a(1-\varepsilon^2)}{1+\varepsilon\cos\theta}
    Фактично, \theta і \rho – полярні координати планети відносно Зірки (центру мас)
  8. Обчислимо декартові координати:
    x = \rho\cos\theta
    y = \rho\sin\theta

Ось і все. Не такий вже й складний алгоритм.

Програма на Python

import math

#import math required
#MeanAnomaly = E - Ecc*sin(E)
def Solve_E(MeanAnomaly, Ecc):
    if Ecc == 0:
        return MeanAnomaly
    E = MeanAnomaly
    for i in range(0,10):
        E = E - (E - MeanAnomaly - Ecc*math.sin(E))/(1 - Ecc*math.cos(E))
    return E

#import math reqired
#MeanAnomality = 2*pi/P
def MeanAnomality(P, t):
    return 2*math.pi/P*t

#import math reqired
#Period of an object mass m travelling around object mass Mass
def Period(a, Mass, m):
    return 2*math.pi*a*math.sqrt(a/(6.67428*10**(-11)*(Mass+m)))

#import math reqired
#True anomality
def Theta(Ecc, E):
    return 2*math.atan(math.sqrt((1+Ecc)/(1-Ecc))*math.tan(E/2))

#import math reqired
#Distance to planet from Sun
def Ro(a, Ecc, theta):
    return a*(1-Ecc**2)/(1+Ecc*math.cos(theta))

def X(a, Ecc, Mass, m, t):
    theta = Theta(Ecc,Solve_E(MeanAnomality(Period(a,Mass,m),t),Ecc))
    return Ro(a,Ecc,theta)*math.cos(theta)

def Y(a, Ecc, Mass, m, t):
    theta = Theta(Ecc,Solve_E(MeanAnomality(Period(a,Mass,m),t),Ecc))
    return Ro(a,Ecc,theta)*math.sin(theta)

Додатково

Очевидно, що алгоритм дещо неоптимізований і для орбіт, наближених до кола (\varepsilon\approx0), нераціональний. Але оптимізація лежить на плечах програміста, а не на моїх =)
Ще один нюанс – період орбіти вираховується динамічно з маси Зірки, планети і довжини великої півосі. В деяких випадках справжній період може відрізнятись від даної динамічної величини, саме тому в астрономії використовується середня аномалія M, як один з основних параметрів орбіти, задаючий положення планети на орбіті. В такому випадку алгоритм треба почати з пункту 4.

About these ads

Written by danbst

Жовтень 17, 2009 at 21:23

5 Responses

Subscribe to comments with RSS.

  1. 1. Ти будеш перекладати функції на C, чи це лежить на плечах програміста? :)

    2. Що таке площина екліптики я вияснив сам. Хоча читачам зручно було б додати посилання на вікіпедію.

    3. Чому тільки x(t) і y(t)? Чи це тільки в площині орбіти. Певне треба дописати для 3d.

    4. І що таке істинна аномалія? Те що ти вже кілька тижднів вивчаєш небесну механіку це ясно, але інші люди взагалі не в курсі діла.

    bunyk

    Жовтень 17, 2009 at 22:54

  2. 1. Так, це зробить програміст =) Пітон – інтуітивно зрозуміла мова програмування, тому перекласти проблем не буде.

    2. Площина екліптики – це площина орбіти Землі. Очевидно, що в інших сонячних системах Землі немає, тому я ввів нове поняття. У мене – це площина XZ.

    3. Так, це мабуть треба дописати – в площині орбіти. Я просто вважаю, що при програмуванні будуть використовуватись повороти матриці, тому і викладаю в такій послідовності – поворот матриці – розрахунок орбіти. В OpenGL це буде просто зробити.

    4. Істинна аномалія – це кут True Anomaly на графіку і обраховується по формулі з пункту 5.

    danbst

    Жовтень 17, 2009 at 23:17

  3. До речі. Колись Ейнтштейн намагвся визначити орбіту Меркурія, вона виявилась дещо дивною, оскільки кут між екліптикою та еліптичною площиною Меркурія був…ем…(див. вікіпедію), і його еліптична орбіта була сильно витягнута, якщо пам’ять мені не зраджує, то тоді була висунута теорія, що на всі небесні тіла діє певний тиск який спричиняється всесвітом, тобто “пустота” має вагу…з ціє теорії випливає, що всесвіт майже обмежений, тобто як циліндр з відкритими кінцями, а якщо вважати це за вірним (те що я написав про всесвіт), то справджується ось цей чудовий парадокс http://ru.wikipedia.org/wiki/Фотометрический_парадокс

    Андрій

    Січень 5, 2010 at 02:46

    • пустота має вагу – це сучасне вчення про Темну Енергію. Щоправда, з цієї теорії не випливає обмеженість Всесвіту (точніше, я не знаю як вивести це твердження).

      danbst

      Січень 8, 2010 at 23:14


Напишіть відгук

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Змінити )

Twitter picture

You are commenting using your Twitter account. Log Out / Змінити )

Facebook photo

You are commenting using your Facebook account. Log Out / Змінити )

Google+ photo

You are commenting using your Google+ account. Log Out / Змінити )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

%d блогерам подобається це: