Created on Sat Nov 26 23:04:26 2016
@author: xiongyiheng
"""
import math
import matplotlib.pyplot as plt
GM=4*(math.pi**2)
perihelion=0.39*(1-0.206)
class Precession:
def __init__(self,e=0.206,time=1,dt=0.0001,vcoefficient=1.1,beta=2.1,alpha=0.01):
self.alpha=alpha
self.beta=beta
self.vcoefficient=vcoefficient
self.e=e
self.a=perihelion/(1-e)
self.x0=self.a*(1+e)
#self.x0=1
self.y0=0
self.vx0=0
self.vy0=self.vcoefficient*math.sqrt((GM*(1-self.e))/(self.a*(1+self.e)))
#self.vy0=2*math.pi
self.X=[self.x0]
self.Y=[self.y0]
self.Vx=[self.vx0]
self.Vy=[self.vy0]
self.T=[0]
self.dt=dt
self.time=time
self.ThetaPrecession=[]
self.TimePrecession=[]
def run(self):
while self.T[-1]<self.time:
r=math.sqrt(self.X[-1]**2+self.Y[-1]**2)
newVx=self.Vx[-1]-(GM*(1+self.alpha/r**2)*self.X[-1]/r**(1+self.beta))*self.dt
newX=self.X[-1]+newVx*self.dt
newVy=self.Vy[-1]-(GM*(1+self.alpha/r**2)*self.Y[-1]/r**(1+self.beta))*self.dt
newY=self.Y[-1]+newVy*self.dt
if abs(newX*newVx+newY*newVy)<0.0014 and r<self.a:
theta=math.acos(self.X[-1]/r)*180/math.pi
if (self.Y[-1]/r)<0:
theta=360-theta
theta=abs(theta-180)
self.ThetaPrecession.append(theta)
self.TimePrecession.append(self.T[-1])
self.Vx.append(newVx)
self.Vy.append(newVy)
self.X.append(newX)
self.Y.append(newY)
self.T.append(self.T[-1]+self.dt)
def show_results(self):
plt.plot(self.X,self.Y,'g')
plt.legend(loc='upper right',frameon=False,fontsize='small')
plt.grid(True)
plt.title('Simulation of the precession of Mercury')
plt.xlabel('x(AU)')
plt.ylabel('y(AU)')
ax = plt.gca()
ax.set_aspect(1)
#plt.xlim(-1,1)
#plt.ylim(-0.6,0.6)
#plt.scatter(0,0)
A = Precession()
A.run()
A.show_results()
网友评论