Побудова орбіт небесних тіл засобами Python

Системи відліку для визначення орбіти

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

Використовуючи це припущення, розглянемо рух однієї і тієї ж точки в двох різних системах відліку К і К’, з яких друга рухається відносно першої з довільною швидкістю — радіус-вектор, що описує положення точки початку системи координат До’ відносно системи відліку До).

Будемо описувати рух точки в системі К’ радіус-вектором , спрямованим з початку координат До системи’ в поточне положення точки. Тоді рух розглянутої точки відносно системи відліку До описується радіус-вектором :

(1)

а відносна швидкість

(2)

де — швидкість точки відносно системи К’; -швидкість руху системи відліку К’ відносно системи відліку К.

Таким чином, для знаходження закону руху точки в довільній системі відліку К необхідно:

1) визначити закон руху точки відносно системи відліку К’ (функцію ;
2) визначити закон руху системи відліку К’ відносно системи відліку До (функцію;
3) визначити закон руху точки відносно системи відліку До, відповідно до (1).

Побудова орбіти Місяця в геліоцентричній системі відліку

У геліоцентричної системі відліку (система) Земля рухається по колу радіуса
R1 = 1.496*10**8 км (період обертання Т1= 3.156*10**7 с.). Місяць, у свою чергу, рухається навколо Землі (система К’) по колу радіуса R2 = 3.844*10**5 км. (період обертання Т2= 2.36*10**6 с. Як відомо [1,2], при русі матеріальної точки по колу радіуса R з постійною кутовою швидкістю координати радіус-вектора, проведеного з початку координат до поточного положення точки, що змінюються за законом:

(3)

де — початкова фаза, що характеризує положення частинки в момент часу t= 0, яку в подальшому ми будемо вважати рівною нулю. Замінюючи в (3) R R1 і R2 і підставляючи в (1), отримуємо залежність радіус-вектора Місяця в геліоцентричної системі координат від часу:

(4)

Вираз (4) задає орбіту Місяця в параметричній формі, де параметром є час. Для побудови шуканої орбіти засобами Python, задамо радіуси орбіт і періоди обертання Землі і Місяця:

З урахуванням (4) визначимо функції залежності координат від часу:

(5)

Використовуючи (5), отримаємо пару координат для орбіти Місяця:

(6)

Задамо кількість точок, в яких обчислюються координати N=1000 і дискретне час на інтервалі періоду обертання Землі dt=T1/N. Напишемо програму і побудуємо графік для позитивної області зміни координат:

Визначення орбіт Землі і Місяця

from numpy import*
from matplotlib.pyplot import*
R1=1.496*10**8
T1=3.156*10**7
R2=3.844*10**5
T2=2.36*10**6
N=1000.0
def X(t):
 return R1*cos(2*pi*t/T1)
def Y(t):
 return R1*sin(2*pi*t/T1)
def x(t):
 return R2*cos(2*pi*t/T2)
def y(t):
 return R2*sin(2*pi*t/T2)
k=100
t=[T1*i/N for i in arange(0,k,1)]
X=array([X(w) for w in t])
Y=array([Y(w) for w in t])
x=array([x(w) for w in t])
y=array([y(w) for w in t])
XG=X+x
YG=Y+y
figure()
title("Траєкторія орбіт Землі і Місяця.n Для позитивних значень координат")
xlabel('X(t),XG(t)')
ylabel('Y(t),YG(t)')
axis([1.2*10**8,1.5*10**8,0,1*10**8])
plot(X,Y,label='Орбіта Землі')
plot(XG,YG,label='Орбіта Місяця')
legend(loc='best')
grid(True)
show()

Отримаємо:

Рис.1

Створений графік дозволяє розширити завдання і подивитися який буде орбіта місяця, якщо радіус орбіти Місяця буде дорівнює R2=3.844*10e7. Внесемо відповідні зміни в програму:

Визначення орбіт Землі і Місяця

 from numpy import*
from matplotlib.pyplot import*
R1=1.496*10**8
T1=3.156*10**7
R2=3.844*10**7
T2=2.36*10**6
N=1000.0
def X(t):
 return R1*cos(2*pi*t/T1)
def Y(t):
 return R1*sin(2*pi*t/T1)
def x(t):
 return R2*cos(2*pi*t/T2)
def y(t):
 return R2*sin(2*pi*t/T2)
t=[T1*i/N for i in arange(0,N,1)]
X=array([X(w) for w in t])
Y=array([Y(w) for w in t])
x=array([x(w) for w in t])
y=array([y(w) for w in t])
XG=X+x
YG=Y+y
figure()
title("Геліоцентрична орбіта Землі і Місяця")
xlabel('X(t),XG(t)')
ylabel('Y(t),YG(t)')
axis([-2.0*10**8,2.0*10**8,-2.0*10**8,2.0*10**8])
plot(X,Y,label='Орбіта Землі')
plot(XG,YG,label='Орбіта Місяця')
legend(loc='best')
grid(True)
show()

Отримаємо:

Рис.2

Порівнюючи орбіти Місяця, представлені на рис. 1 і 2, виявляємо їх суттєві відмінності. Для пояснення причини цих відмінностей необхідно порівняти лінійні швидкості руху Місяця в першому і в другому випадку і лінійну швидкість руху Землі.

Так як напрямок лінійної швидкості руху Землі відносно Сонця, як і напрямок лінійної швидкості руху Місяця відносно Землі, змінюється у часі, а швидкість залишається постійною за величиною.

В якості кількісної характеристики співвідношення лінійних швидкостей руху Місяця і Землі в геліоцентричної системі координат слід вибрати різниця між модулем лінійної швидкості руху Землі і проекцією лінійної швидкості Місяця на напрям вектора лінійної швидкості Землі:

(7)

Визначимо функції, що описують закони зміни складових швидкості Землі і Місяця:

(8)

Щоб визначити результуючу швидкість з урахуванням проекції, скористаємося співвідношенням:

(9)

Напишемо програму з урахуванням(5), (8), (9) і радіуса орбіти Місяця R2=3.844*10**5 км.:

Місяць і Земля рухаються в одному напрямку

from numpy import*
from matplotlib.pyplot import*
R1=1.496*10**8
T1=3.156*10**7
R2=3.844*10**5
T2=2.36*10**6
N=1000.0
k1=2*pi/T1
k2=2*pi/T2
def Vx(t):
 return -k1*R1*sin(k1*t)
def Vy(t):
 return k1*R1*cos(k1*t)
def vx(t):
 return -k2*R2*sin(k2*t)
def vy(t):
 return k2*R2*cos(k2*t)
def D(t):
 return sqrt(Vx(t)**2+Vy(t)**2)-sqrt(vx(t)**2+vy(t)**2)*(Vx(t)*vx(t)+Vy(t)*vy(t))/((sqrt(Vx(t)**2+Vy(t)**2))*(sqrt(vx(t)**2+vy(t)**2)))
x=[T1*i/N for i in arange(0,N,1)]
y=[D(t) для t in x]
title("Місяць рухається в одному напрямку з Землею n Радіус орбіти Місяця R2=3.844*10**5 км.")
xlabel('t')
ylabel('D(t)')
plot(x,y)
show()

Отримаємо:

Рис.3.

Напишемо програму з урахуванням (5), (8), (9) і радіуса орбіти Місяця R2=3.844*10**7 км:

Місяць періодично рухається в протилежному напрямку до Землі

from matplotlib.pyplot import*
R1=1.496*10**8
T1=3.156*10**7
R2=3.844*10**7
T2=2.36*10**6
N=1000.0
k1=2*pi/T1
k2=2*pi/T2
def Vx(t):
 return -k1*R1*sin(k1*t)
def Vy(t):
 return k1*R1*cos(k1*t)
def vx(t):
 return -k2*R2*sin(k2*t)
def vy(t):
 return k2*R2*cos(k2*t)
def D(t):
 return sqrt(Vx(t)**2+Vy(t)**2)-sqrt(vx(t)**2+vy(t)**2)*(Vx(t)*vx(t)+Vy(t)*vy(t))/((sqrt(Vx(t)**2+Vy(t)**2))*(sqrt(vx(t)**2+vy(t)**2)))
x=[T1*i/N for i in arange(0,N,1)]
y=[D(t) для t in x]
title(" Періодично Місяць рухається в протилежному до Землі n напрямку. Радіус орбіти Місяця R2=3.844*10**7 км.")
xlabel('t')
ylabel('D(t)')
plot(x,y)
show()

Отримаємо:

Рис.4.

Аналіз залежностей дозволяє пояснити причину відмінностей орбіт. Функція D(t) при R2 =3,844*10**5 км завжди позитивна, тобто Місяць завжди рухається в напрямку руху Землі і петлі не утворюються. При R2 = = 3,844*10**7 км величина D(t) приймає негативні значення, тобто існують моменти часу, в які Місяць рухається в напрямку, протилежному напрямку руху Землі, а тому орбіта має петлі.

Побудова орбіти Марса в системі відліку, пов’язаної з Землею

.

У геліоцентричної системі відліку (система) Земля рухається по колу радіуса R1= 1.496*10**8 км, період обертання Т1= 365.24 доби, Марс рухається по еліпсу, велика піввісь якого ам = 2.28*10**8 км, період обертання Марса Тм = 689.98 сут., ексцентриситет орбіти е = 0.093 [3]. Рух Землі описується радіус-вектором R(t), що задається виразом (3). У зв’язку з тим, що орбіта Марса є еліпсом, залежності x=x(t), y=y(t) від часу задаються параметрично [4]:

(10)

(11)

(12)

Повного обороту по еліпсу відповідає зміна параметра від 0 до. Для побудови орбіти Марса необхідно обчислити в одні і ті ж моменти часу, координати радіус-векторів, що описують положення Землі та Марса в геліоцентричної системі відліку, потім відповідно до співвідношенням обчислити координати Марса в системі відліку, пов’язаної з Землею.

Для побудови орбіти Марса в системі відліку, пов’язаної з Землею скористаємося раніше наведеними параметрами орбіти Землі і Марса, співвідношеннями (10)-(12), а також співвідношеннями для координат Землі:

(13)

(14)

Слід врахувати, що число періодів обертання Марса навколо Сонця дорівнює K=9, тоді кількість точок, в яких слід зробити розрахунок і відстань між ними, будуть визначатися з співвідношень:

(15)

Орбіта Марса в системі відліку Землі

from numpy import*
from matplotlib.pyplot import*
R1=1.496*10e8
T1=365.24
am=2.28*10e8
Tm=689.98
ee=0.093
N=36000
def x(g):
 return am*(cos(g)-ee)
def y(g):
 return am*sqrt(1-ee**2)*sin(g)
def t(g):
 return Tm*(g-ee*sin(g))/2*pi
def X(g):
 return R1*cos(2*pi*t(g)/T1)
def Y(g):
 return R1*sin(2*pi*t(g)/T1)
y=array([y(2*pi*i/N) for i in arange(0,N,1)])
x=array([x(2*pi*i/N) for i in arange(0,N,1)])
X=array([X(2*pi*i/N) for i in arange(0,N,1)])
Y=array([Y(2*pi*i/N) for i in arange(0,N,1)])
t=array([t(2*pi*i/N) for i in arange(0,N,1)])
figure()
title("Гелиоцентрические орбіти Землі і Марса")
xlabel('x(g),X(g)')
ylabel('y(g),Y(g)')
plot(x,y,label='ОрбитаМарса')
plot(X,Y,label='ОрбитаЗемли')
legend(loc='best')
figure()
title("Положення Марса в системі відліку, пов'язаної з Землею")
xlabel('x1/10e8')
ylabel('y1(g/10e8')
x1=(x-X)
y1=(y-Y)
plot(x1/10e8,y1/10e8)
figure()
title("Залежність відстані між Землею і Марсом n від часу у роках)
xlabel('t/365.24')
ylabel('sqrt(x1**2+y1**2)/10e8')
y2=sqrt(x1**2+y1**2)/10e8
x2=t/365.24
plot(x2,y2)
show()

Отримаємо:

Рис.5

Обчислимо координати радіус-вектора, що описує положення Марса в системі відліку, пов’язаної з Землею, і побудуємо орбіти (Мал.6), використовуючи співвідношення:

(16)

Рис.6

Ще однією важливою характеристикою руху Марса (в першу чергу для міжпланетних космічних польотів) є відстань між Землею і Марсом s(t), яка визначається модулем радіус-вектора, що описує положення Марса в системі відліку, пов’язаної з Землею. Залежність відстані між Землею і Марсом від часу, вимірюваного в земних роках, представлена на рис.7.

Рис.7

Аналіз залежності, представленої на рис.7, показує, що відстань між Землею і Марсом є складною періодичною функцією часу. Якщо скористатися термінологією теорії сигналів [5], то про залежності s(t) можна сказати, що вона являє собою амплітудно-модульований сигнал, який прийнято представляти у вигляді добутку двох функцій високочастотної (несучої) і низькочастотної функції, що задає амплітудну модуляцію (огинаючої):

(17)

де — постійна складова функції u(t); а — амплітуда сигналу; — частота несучої; — амплітуда функції, що задає глибину амплітудної модуляції; — частота модулюючим функції.

З рис.7 видно, що період несучої складає приблизно 2 роки, період модулюючим функції приблизно 17 років]6].

Побудова геліоцентричної орбіти комети Галлея

В останній раз комета Галлея проходила через свій перигелій (найближча до Сонця точка орбіти) 9 лютого 1986 року. (Саме Сонце вважається розташованим в початку координат.)

Координати і компоненти швидкості комети Галлея в той момент були рівні P0 = (0.325514, 0.459460, 0.166229) і v0 = (-9.096111, -6.916686, -1.305721) відповідно, причому відстань тут виражено в астрономічних одиницях довжини – а.е.д., або просто а.е. (астрономічна одиниця, тобто довжина великої головною півосі земної орбіти), а час – у роках. В цих одиницях виміру тривимірні рівняння руху комети мають вигляд:

(18)

де:,

Побудова геліоцентричної орбіти комети Галлея

from numpy import*
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.patches import Circle 
def f(y, t):
 y1, y2, y3, y4,y5,y6 = y
 return [y2, -(4*pi*pi*y1)/(y1**2+y3**2 +y5**2)**(3/2),y4,-(4*pi*pi*y3)/(y1**2+y3**2 +y5**2)**(3/2),y6,-(4*pi*pi*y5)/(y1**2+y3**2 +y5**2)**(3/2)] 
t = linspace(0,300,10001)
y0 = [0.325514,-9.096111, -0.459460,-6.916686,0.166229,-1.305721]
[y1,y2, y3, y4,y5,y6]=odeint(f, y0, t, full_output=False).T
fig, ax = plt.subplots()
plt.title("Орбіта комети Галлея(відстань а.е., час в роках) n Сонце в центрі координат")
plt.xlabel('x(t)')
plt.ylabel('y(t)')
fig.set_facecolor('white')
ax.plot(y1,y3,linewidth=1)
circle = Circle((0, 0), 0.2, facecolor='orange') 
ax.add_patch(circle)
plt.axis([1,-21,-1,29])
plt.grid(True)
fig, ax = plt.subplots()
plt.title("Орбіта комети Галлея n Сонце в центрі координат")
plt.xlabel('x(t)')
plt.ylabel('z(t)')
fig.set_facecolor('white')
ax.plot(y1,y5,linewidth=1)
circle = Circle((0, 0), 0.1, facecolor='orange') 
ax.add_patch(circle)
plt.axis([1,-21,1,-11])
plt.grid(True)
fig, ax = plt.subplots()
plt.title("Орбіта комети Галлея n Сонце в центрі координат")
plt.xlabel('y(t)')
plt.ylabel('z(t)')
fig.set_facecolor('white')
ax.plot(y3,y5,linewidth=1)
circle = Circle((0, 0), 0.2, facecolor='orange') 
ax.add_patch(circle)
plt.axis([-1,29,1,-11])
plt.grid(True)
fig, ax = plt.subplots()
plt.title("Проекція швидкості руху комети Галлея n на площині ZOX і ZOY ")
ax.plot(t,y1,linewidth=1)
ax.plot(t,y3,linewidth=1)
plt.show()

Отримаємо:

Ваша власна комета

Спробуйте провести експеримент. У ніч ви встановите свій телескоп на вершині, розташованої недалеко від вашого будинку височини. Ніч повинна бути ясною, безхмарним, зоряної і, якщо вам посміхнулася фортуна: 0 годин 30 хвилин ночі ви помітите нову комету.

Після повторних спостережень в наступні ночі вам вдасться вирахувати її координати в ту першу ніч. Координати геліоцентричної системі координат: P0= (x0, y0, z0) і вектор швидкості v0 = (vx0, vy0, vz0).

Використовуючи ці дані, визначте:

• відстань комети від Сонця в перигелії (найближча до Сонця точка орбіти) і в афелії (найдальша від Сонця точка орбіти);
• швидкості комети під час проходження через перигелій і через афелій;
• період обертання комети навколо Сонця;
• наступні дві дати проходження комети через перигелій.

Якщо вимірювати відстань у астрономічних одиницях, а час — у роках, то рівняння руху комети приймуть вигляд (18). Для вашої власної комети виберіть довільні початкові координати і швидкості того ж порядку, що і у комети Галлея.

В разі потреби, повторно здійснюйте довільний вибір початкового положення вектора швидкості до тих пір, поки не отримаєте правдоподібну ексцентричну орбіту, що виходить за межі орбіти Землі (як у більшості цих комет).

Посилання:

 

1. Фейнман Р., Лейтон Р., Сендс М. Фейнмановские лекції з фізики. 3-е изд. Т. 1.-2. М.: Світ, 1977.
2. Матвєєв А. Н. Механіки і теорії відносності. М.: Высш. шк., 1986.
3. Фізична енциклопедія. Т. 3. М.: Велика російська енциклопедія, 1992.
4. Ландау Л. Д., Ліфшиц Е. М. Курс теоретичної фізики. Механіка. М.: Фю-матгиз, 1958.
5. Баскаков С. В. Радіотехнічні ланцюги і сигнали. М.: Высш. шк., 1988.
6. Поршнєв C. Ст. Комп’ютерне моделювання фізичних процесів з використанням пакету mathcad.