Úloha 5

Kdyby se Země najednou zastavila na své dráze, za jakou dobu $τ$ by dopadla do Slunce?

import numpy as np
import matplotlib.pyplot as plt

konstanty: hmotnost měříme v jednotkách hmotnosti Země, čas měříme v letech

Ms=3.33166e5 #hmotnost Slunce v jednotkach hmotnosti Zeme
G=1.18e-4 #gravitacni konstanta AU^3/rok Mz
Rs=0.00465 #poloměr Slunce v AU

počáteční podmínky

dt=0.001 #časový krok 0.01 roku
t=[] #pole časů
x=[] #pole souřadnic (je to jednorozměrný problém)
v=[] #pole rychlostí
t=np.append(t,0) #počáteční čas
x=np.append(x,1) #počáteční poloha Země v AU
v=np.append(v,-G*Ms/x[0]**2*dt/2) #rychlost Země čase dt/2

iterativní řešení pohybové rovnice
$\ddot{x}=-G\frac{M_s}{x^2}$

i=0
while True:
    t=np.append(t,t[i]+dt)
    x=np.append(x,x[i]+v[i]*dt)
    if np.min(x)<=Rs:
        break
    v=np.append(v,v[i]-G*Ms/x[i+1]**2*dt)
    i+=1
print(f"Země, která se zastaví, narazí do Slunce za {365*t[i]:6.3f} dní")    

Země, která se zastaví, narazí do Slunce za 64.605 dní

plt.plot(365*t,x) #čas vyneseme ve dnech
plt.xlim(0,)
plt.ylim(0,)
plt.xlabel("čas (dny)")
plt.ylabel("vzdálenost od Slunce (AU)")
plt.show()

plt.plot(365*t[0:np.size(t)-1],abs(v))
plt.xlim(0,)
plt.ylim(0,)
plt.xlabel("čas (dny)")
plt.ylabel("rychlost Země (AU/rok)")
plt.show()