The tenth homework of computational physics about $T^{2}/a^{3}$ of noncircular orbits of 4.8

defeng kong 2014301020081

Abstract

Background

4.8

Verify Kepler's third law for elliptical orbits.Run the planetary motion program with initial conditions chosen to give orbits that are noncircular.Calculate and compare with the values $T^{2}/a^{3}$ given in Table 4.2.

The main body and conclusion

the code is:

import matplotlib.pyplot as plt 
import numpy as np 
class noncircular(): 
    def __init__(self,e,a,y0,vx0,M_P,M_S,dt): 
        self.e=e
        self.a=a
        self.x = e*a+a 
        self.y = y0 
        self.vx =vx0
        self.M_P = M_P
        self.M_S = M_S
        self.vy = np.sqrt(4*np.pi**2*(1-self.e)*(1+self.M_P/self.M_S)/(self.a*(1+self.e)))
        self.t=0 
        self.dt=dt
        self.tra_x = [a] 
        self.tra_y = [y0]
        self.W=0
        self.b=1
    def calculate(self):
       for i in  range(30000):
           self.r=np.sqrt(self.x**2+self.y**2)
           self.vx=self.vx-4*np.pi**2*self.x*self.dt/self.r**3
           self.vy=self.vy-4*np.pi**2*self.y*self.dt/self.r**3
           self.x=self.x+self.vx*self.dt
           self.y=self.y+self.vy*self.dt
           self.tra_x.append(self.x) 
           self.tra_y.append(self.y)
           self.t=self.t+self.dt    
    def run_trajectory(self):
        self.W = (0.001*2*np.pi*max(self.tra_x)*max(self.tra_y)\
        /((self.e*self.a+self.a)*self.tra_y[1]))**2/self.a**3
        print self.W
        plt.axis('equal')
        plt.ylim(-max(self.tra_y),max(self.tra_y))
        plt.xlabel('x(AU)') 
        plt.ylabel('y(AU)') 
        #plt.title('Simulation of circular orbit ')
        plt.title('Simulation of elliptical orbit ')
        plt.text(0,0.5,'T^2/a^3=1.11520692458')
        plt.plot(self.tra_x,self.tra_y,'.')
        plt.show()
A=noncircular(0.056,9.54,0,0,5.7*10**26,2.0*10**30,0.001)
A.calculate()
A.run_trajectory()

For circular orbit:
for Venus:

for Earth:

for Mars(Venus in plot is wrong):

for Jupiter:

for Saturn:

For elliptical orbit
for Venus:

for Earth(Venus in plot is wrong):

for Mars(Venus in plot is wrong):

for Jupiter:

for Saturn:

so for circular orbit(the value of T^2/a^3):
Venus:0.9898
Earth:0.9938
Mars:0.9967
Jupiter:1.0004
Saturn:1.00007
For elliptical orbit:(the value of T^2/a^3)
Venus:1.0039
Earth:1.0280
Mars:1.1913
Jupiter:1.0987
Saturn:1.1152
We can see the value is close to 1 for circular orbit and elliptical orbit.

Thanks