Keplerova úloha

výpočet trajektorie Země obíhající okolo Slunce

import numpy as np
import matplotlib.pyplot as plt

Počátek soustavy souřadnic umístíme do Slunce a zvolíme vhodné jednotky:

• hmotnost budeme měřit v jednotkách hmotnosti Země, $M_Z = 1$
• čas budeme měřit v rocích, tj. perioda oběhu Země okolo slunce, $T_Z = 1$
• vzdálenost budeme měřit v Astronimuckých jednotkách (AU), tj. za jednotkovou vzdálenost bereme střední vzdálenost Země od Slunce

konstanty (ve výše uvedených jednotkách)

Mz = 1.0 #hmotnost Země
Ms = 333166.0 #hmotnost Slunce (v jednotkach Mz)
G = 1.18e-4 #gravitační konstanta (Mz^-1 AU^3 rok^-2)
T = 1 #doba obehu (1 rok)
v0 = 6.166 #pocatecni rychlost Země v aféliu  (AU/rok)

časový krok pro numerické řešení pohybové rovnice

dt = T/10000

Vytvoříme pole pro časy, souřadnice, vzdálenost Země od Slunce, složky rychlosti, velikost rychlosti a plošnou rychlost

t=np.arange(0,T+dt,dt) #pole časů
x=np.empty(np.size(t)) #x-ové souřadnice
y=np.empty(np.size(t)) #y-ové souřadnice
r=np.empty(np.size(t)) #vzdálenost Země od Slunce
vx=np.empty(np.size(t)) #x-ová složka rychlosti
vy=np.empty(np.size(t)) #y-ová složka rychlosti
v=np.empty(np.size(t)) #velikost rychlosti
vp=np.empty(np.size(t)) #plošná rychlost

Řešíme pohybovou rovnici danou Newtonovým gravitačním zákonem $$M_Z \vec{a} = -G\frac{M_S M_Z}{r^2}\frac{\vec{r}}{r}$$ Rozepsáno do složek je to následující soustava diferenciálních rovnic: \begin{align} \ddot{x} &= -G\frac{M_S}{r^3}x,\\ \ddot{y} &= -G\frac{M_S}{r^3}y,\\ \end{align} kde $$r = \sqrt{x^2 + y^2}$$

počáteční podmínky

x[0] = 1.0167 #zaciname v afeliu
y[0] = 0.0
r[0] = np.sqrt(x[0]**2+y[0]**2) #vzdálenost Země od Slunce
vx[0] = -G*Ms/r[0]**3*x[0]*dt/2
vy[0] = v0-G*Ms/r[0]**3*y[0]*dt/2
v[0] = np.sqrt(vx[0]**2+vy[0]**2)
vp[0] = 0.5*r[0]*v[0]

Numerické řešení pohybové rovnice

rmin = x[0] #do této proměnné uložíme nejmenší vzdálenost Země od Slunce
rmax = 0    #do této proměnné uložíme největší vzdálenost Země od Slunce
for i in range(np.size(t)-1):
    x[i+1] = x[i]+vx[i]*dt #updatování x-ové souřadnice t -> t+dt 
    y[i+1] = y[i]+vy[i]*dt #updatování y-ové souřadnice t -> t+dt 
    r[i+1] = np.sqrt(x[i+1]**2+y[i+1]**2) #aktuální vzdálenost Země od Slunce
    vx[i+1] = vx[i]-G*Ms/r[i+1]**3*x[i+1]*dt #updatování x-ové složky rychlosti t -> t+dt 
    vy[i+1] = vy[i]-G*Ms/r[i+1]**3*y[i+1]*dt #updatování y-ové složky rychlosti t -> t+dt 
    v[i+1] =np.sqrt(vx[i+1]**2+vy[i+1]**2) #aktuální velikost rychlosti Země
    vp[i+1] = 0.5*(v[i]+v[i+1])/2*r[i+1] #aktuální velikost plošné rychlosti Země
    if r[i+1]<rmin: #najdeme nejmenší vzdálenost Země od Slunce a uložíme ji do proměnné rmin    
        rmin = r[i+1]
    if r[i+1]>rmax: #najdeme největší vzdálenost Země od Slunce a uložíme ji do proměnné rmax 
        rmax = r[i+1]        

Nakreslíme graf trajektorie

#trajektorie
fig,ax=plt.subplots(figsize=(5,5))
plt.title('trajektorie',fontsize=14)
ax.plot(x,y)
ax.scatter(0,0,s=50,marker='+',c='red')
plt.xlim(-1.1,1.1)
plt.ylim(-1.1,1.1)
ax.set_xlabel('x (AU)',fontsize=14)
ax.set_ylabel('y (AU)',fontsize=14)
plt.show()

Vypočítáme hlavní a vedlejší poloosu a excentricitu oběžné dráhy Země okolo Slunce

e = (rmax-rmin)/2 #excentricita
a = (rmax+rmin)/2 #hlavní poloosa
b = np.sqrt(a**2-e**2) #vedlejší poloosa
epsilon = e/a #numerická excentricita
print (f"hlavní poloosa: {a:8.6f} AU,    vedlejší poloosa: {b:8.6f} AU") 
print (f"excentricita:   {e:5.3f} AU, numerická excentricita: {epsilon:5.4f}") 

hlavní poloosa: 0.999936 AU,    vedlejší poloosa: 0.999795 AU
excentricita:   0.017 AU, numerická excentricita: 0.0168

časový průběh vzdálenosti Země od Slunce

#vzdalenost od slunce
fig,ax=plt.subplots(figsize=(5,5))
plt.title('vzdálenost od Slunce',fontsize=14)
plt.xlim(0,1)
ax.plot(t,r,c='green')
ax.set_xlabel('t (rok)',fontsize=14)
ax.set_ylabel('r (AU)',fontsize=14)
plt.show()

časový průběh velikosti rychlosti

#velikost rychlosti
fig,ax=plt.subplots(figsize=(5,5))
plt.title('velikost rychlosti',fontsize=14)
ax.plot(t,v,c='red')
plt.xlim(0,1)
ax.set_xlabel('t (rok)',fontsize=14)
ax.set_ylabel('v (AU/rok)',fontsize=14)
plt.show()

časový průběh plošné rychlosti

#plosna rychlost
fig,ax=plt.subplots(figsize=(5,5))
plt.title('plošná rychlost',fontsize=14)
plt.xlim(0,1)
plt.ylim(3,3.3)
ax.plot(t,vp,c='violet')
ax.set_xlabel('t (rok)',fontsize=14)
ax.set_ylabel('vp (AU${}^2$/rok)',fontsize=14)
plt.show()