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ě, MZ=1
- čas budeme měřit v rocích, tj. perioda oběhu Země okolo slunce, TZ=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 MZ→a=−GMSMZr2→rr Rozepsáno do složek je to následující soustava diferenciálních rovnic: ¨x=−GMSr3x,¨y=−GMSr3y, kde r=√x2+y2
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()