Побудова орбіт небесних тіл засобами 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). Для вашої власної комети виберіть довільні початкові координати і швидкості того ж порядку, що і у комети Галлея.
В разі потреби, повторно здійснюйте довільний вибір початкового положення вектора швидкості до тих пір, поки не отримаєте правдоподібну ексцентричну орбіту, що виходить за межі орбіти Землі (як у більшості цих комет).
Посилання:
- Фейнман Р., Лейтон Р., Сендс М. Фейнмановские лекції з фізики. 3-е изд. Т. 1.-2. М.: Світ, 1977.
- Матвєєв А. Н. Механіки і теорії відносності. М.: Высш. шк., 1986.
- Фізична енциклопедія. Т. 3. М.: Велика російська енциклопедія, 1992.
- Ландау Л. Д., Ліфшиц Е. М. Курс теоретичної фізики. Механіка. М.: Фю-матгиз, 1958.
- Баскаков С. В. Радіотехнічні ланцюги і сигнали. М.: Высш. шк., 1988.
- Поршнєв C. Ст. Комп’ютерне моделювання фізичних процесів з використанням пакету mathcad.
