Nature:Physics

克普勒(Kepler)行星運動定律

太陽系的行星數據

行星繞日的軌道是橢圓形

克卜勒的行星運動第一定律表示：任何行星繞日運動時的軌道是橢圓軌道，而太陽就在橢圓的兩個焦點其中之一，因此行星繞日運動時有所謂的遠日點和近日點，橢圓的離心率就在標示出這個軌道和圓的軌道有多大的差異。大部分的行星繞日的軌道都非常接近圓，也就是離心率接近於0，對於彗星的軌道而言離心率就非常的大，我們所熟知的哈雷彗星其離心率約為0.967，週期約為75.3年，其軌道離開太陽的距離(遠日點=35.1 AU, AU=astronomical unit=天文單位，地球至太陽的距離)可以達到海王星(30 AU)與冥王星(39 AU)之間。我們可以利用程式模擬的方式來顯現出行星繞日運動的橢圓軌道，行星運動的動力就是牛頓的重力定律，而運動的法則就是牛頓第二運動定律：\[\vec{F}=m\vec{a}, \frac{d\vec{v}}{dt}=\frac{\vec{F}}{m}=-\frac{GM_s}{r^3}\vec{r}, \frac{d\vec{r}}{dt}=\vec{v} \]
運動的差分近似 \[ \vec{v}(t+dt) \simeq \vec{v}(t)+\vec{a}(t)dt \] \[ \vec{r}(t+dt) \simeq \vec{r}(t)+\vec{v}(t)dt \]

程式樣本1, 以Vpython模擬行星運動

GlowScript 2.9 VPython
orb=['mercury','venus','earth','mars','jupiter','saturn','uranus','neptune','pluto','halley']
ra=[0.467,0.728,1.017,1.666,5.459,10.124,20.078,30.386,49.304,35.082]
rp=[0.307,0.719,0.983,1.381,4.95,9.041,18.324,29.709,29.658,0.586]
a=[0.387,0.724,1,1.523,5.204,9.582,19.201,30.047,39.481,17.834]
va=[1.303,1.168,0.983,0.738,0.417,0.305,0.218,0.18,0.123,0.031]
e=[0.207,0.006,0.017,0.094,0.049,0.057,0.046,0.011,0.249,0.967]
T=[0.241,0.615,1,1.88,11.86,29.46,84.01,164.8,247.7,75.32]

NP=len(rp)-2; NP=10
RE=1.5e11; TE=365*86400; VE=2.*pi*RE/TE
MS=1.989e30; ME=5.972e24; G=6.67e-11; V=[]
print(' i      planets      ra         rp          a         va          e          T')
for i in range(NP):
    ai=(ra[i]+rp[i])/2.
    rat=rp[i]/ra[i]/((rp[i]+ra[i])/2.)
    vir=sqrt(rat)
    ei=abs(1-((vir*2*pi)**2*(ra[i])/4/pi**2))
    PL='{:4d}'.format(i)+'{:10s}'.format(orb[i])+'{:9.3f}'.format(ra[i]) \
+'{:9.3f}'.format(rp[i])+'{:9.3f}'.format(ai)+'{:9.3f}'.format(vir) \
+'{:9.3f}'.format(ei)+'{:9.3f}'.format(T[i])
    print(PL)
    vi=VE*vir
    V.append(vi)
ran=random

scene = canvas(width=800, height=700, center=vector(0,0,0),
                background=vector(0,0,0))
            
SUN = sphere(radius=0.1*RE,pos=vector(0,0,0),color=color.yellow)
planets=[sphere(radius=0.2*RE, pos=vector(ra[i]*RE,0,0), make_trail=True) for i in range(NP)]
planets[0].radius=0.02*RE
planets[1].radius=0.05*RE
planets[2].radius=0.1*RE
planets[3].radius=0.05*RE
planets[4].radius=0.3*RE

planets[0].color=color.cyan
planets[1].color=color.orange
planets[2].color=color.blue
planets[3].color=color.red
planets[4].color=color.orange
planets[5].color=color.green
planets[6].color=color.white
planets[7].color=color.purple
planets[8].color=color.white
planets[9].color=color.yellow
for i in range(NP):
    planets[i].trail_color=planets[i].color
    planets[i].v=vector(0,V[i],0)
t=0.
time='{:4.0f}'.format(t/TE)
T=label(text=time,pos=vector(40*RE,10*RE,0))
Ns=10000
dt = TE/Ns; NR=Ns*100
scene.waitfor('click')
while True:
    rate(NR)
    time='{:4.0f}'.format(t/TE)
    T.text=time
    for i in range(2,NP):
        planets[i].v += -dt*G*MS/mag(planets[i].pos)**3 * planets[i].pos
        planets[i].pos += planets[i].v * dt    
    t=t+dt
    if(t/TE > 80): break

等面積定律:

行星至太陽間的連心線在相同的時間內掃過相同的面積

程式樣本2, 哈雷彗星與等面積定律

GlowScript 2.9 VPython
orb=['earth','halley']
ra=[1.017,35.082]
rp=[0.983,0.586]
a=[1,17.834]
va=[0.983,0.031]
e=[0.017,0.967]
T=[1,75.32]
NP=2; RE=1.5e11; TE=365*86400; VE=2.*pi*RE/TE
MS=1.989e30; ME=5.972e24; G=6.67e-11; V=[]
print(' i      planets      ra         rp          a         va          e          T')
for i in range(NP):
    ai=(ra[i]+rp[i])/2.
    rat=rp[i]/ra[i]/((rp[i]+ra[i])/2.)
    vir=sqrt(rat)
    ei=abs(1-((vir*2*pi)**2*(ra[i])/4/pi**2))
    PL='{:4d}'.format(i)+'{:10s}'.format(orb[i])+'{:9.3f}'.format(ra[i])+'{:9.3f}'.format(rp[i])+'{:9.3f}'.format(ai) \
    +'{:9.3f}'.format(vir)+'{:9.3f}'.format(ei)+'{:9.3f}'.format(T[i])
    print(PL)
    vi=VE*vir
    V.append(vi)
scene = canvas(width=1200, height=600, center=vector(ra[1]*RE/2.,0,0),
                background=vector(0,0,0))    
SUN = sphere(radius=0.3*RE,pos=vector(0,0,0),color=color.yellow)
planets=[sphere(radius=0.3*RE, pos=vector(ra[i]*RE,0,0), make_trail=True) for i in range(NP)]
planets[0].radius=0.1*RE
planets[1].radius=0.1*RE
planets[0].color=vec(0.5,0.5,1)
planets[1].color=color.red
for i in range(NP):
    planets[i].trail_color=planets[i].color
    planets[i].v=vector(0,V[i],0)

t=0.
time='{:02.0f}'.format(t/TE)
T=label(text=time,pos=vector(20*RE,7*RE,0))
arrow(pos=vec(0,0,0),axis=planets[1].pos,shaftwidth=0.02*RE,headwidth=0.05*RE,color=color.yellow)
Ns=20000
dt = TE/Ns; NR=Ns*1000; Y1=0
scene.waitfor('click')
while True:
    rate(NR/4)
    time='{:02.0f}'.format(t/TE)
    Y=int(t/(4*TE))
    if(Y > Y1):
        Y1=Y
        arrow(pos=vec(0,0,0),axis=planets[1].pos,shaftwidth=0.02*RE,headwidth=0.05*RE,color=color.yellow)
    T.text=time
    for i in range(NP):
        planets[i].v += -dt*G*MS/mag(planets[i].pos)**3 * planets[i].pos
        planets[i].pos += planets[i].v * dt    
    t=t+dt
    if(t/TE > 80): break

克卜勒的行星運動第三定律：週期的平方正比於距離的三次方(自檔案讀入數據範例)

在下面這個例子當中，我們要從外部的檔案名稱為(solar.txt)的檔案讀進我們太陽系8個行星的數據，透過這些數據我們可以檢驗克卜勒的行星運動第三定律：週期的平方正比於距離的三次方。

程式樣本3, 克卜勒的行星運動第三定律

from visual import *
#f=open('solar.txt','r')
#a=f.read()
a='\
Mercury	0.39	0.24 \
Venus	0.72	0.62 \
Earth	1	1 \
Mars	1.52	1.88 \
Jupiter	5.2	11.86 \
Saturn	9.54	29.46 \
Uranus	19.22	84.01 \
Neptune	30.06	164.8 \
'
print a
b=a.split()
lenb=len(b)
np=8
print b,lenb
dat1=[]
for i in range(len(b)):
    if(i%3 ==0):  dat1.append(b[i])
    else:   dat1.append(float(b[i]))
planet=[dat1[i*3] for i in range(np)]
R=[dat1[i*3+1] for i in range(np)]
T=[dat1[i*3+2] for i in range(np)]
s3="%10.3f"*3
sqx=[]
for i in range(8):
    sqx.append(T[i]**2/R[i]**3)
    print (s3 %(T[i],R[i],sqx[i]))

太陽系的行星運動

克普勒(Kepler)行星運動定律

行星繞日的軌道是橢圓形

等面積定律:

克卜勒的行星運動第三定律：週期的平方正比於距離的三次方(自檔案讀入數據範例)