@author: xyh
"""
import math
import matplotlib.pyplot as plt
GM=4(math.pi*2)
class precession:
def init(self,e=0.5,time=9,dt=0.0001,vcoefficient=1,beta=2,alpha=0,intial_w=0,intial_sita=0):
self.intial_sita=intial_sita
self.alpha=alpha
self.beta=beta
self.vcoefficient=vcoefficient
self.e=e
self.a=1/(1+e)
self.x0=self.a*(1+e)
self.y0=0
self.vx0=0
self.vy0=self.vcoefficientmath.sqrt((GM(1-e))/(self.a*(1+e)))
self.X=[self.x0]
self.Y=[self.y0]
self.Vx=[self.vx0]
self.Vy=[self.vy0]
self.w=[intial_w]
self.sita=[intial_sita]
self.T=[0]
self.dt=dt
self.time=time
self.ThetaPrecession=[]
self.TimePrecession=[]
def calculate_reset(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/r2)self.X[-1]/r*(1+self.beta))self.dt
newX=self.X[-1]+newVx*self.dt
newVy=self.Vy[-1]-(GM(1+self.alpha/r2)self.Y[-1]/r*(1+self.beta))self.dt
newY=self.Y[-1]+newVy*self.dt
newW=self.w[-1]-3GM/(r5)(self.X[-1]math.sin(self.sita[-1])-self.Y[-1]math.cos(self.sita[-1]))(self.X[-1]math.cos(self.sita[-1])+self.Y[-1]math.sin(self.sita[-1]))self.dt
newSita=self.dt*self.w[-1]+self.sita[-1]
self.Vx.append(newVx)
self.Vy.append(newVy)
self.X.append(newX)
self.Y.append(newY)
self.w.append(newW)
self.sita.append(newSita)
self.T.append(self.T[-1]+self.dt)
while self.sita[-1]>=1*math.pi:
self.sita[-1]=self.sita[-1]-2*math.pi
while self.sita[-1]<=-1*math.pi:
self.sita[-1]=self.sita[-1]+2*math.pi
def plot(self,slogan=''):
plt.scatter(self.sita,self.w,label=slogan,s=0.1)
plt.legend(loc='upper right',frameon=False,fontsize='small')
plt.grid(True)
plt.title('Hyperion(Elliptical orbit) $\omega$ versus $\theta$')
plt.ylabel('$\omega(radians/yr)$')
plt.xlabel('$\theta(radians)$')
