附二(托勒密与开普勒轨道的比较代码,其中 RK4代码源自 https://rosettacode.org/wiki/Runge-Kutta_method )
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 22 2018
@author: YDX
"""
import matplotlib as mpl
#mpl.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
from matplotlib import animation
import numpy as np
plt.clf()
mpl.rcParams['animation.ffmpeg_path'] = "C:/Program Files/ImageMagick-7.0.8-Q16/ffmpeg.exe"
mpl.rcParams['animation.convert_path'] = "C:/Program Files/ImageMagick-7.0.8-Q16/magick.exe"
mpl.rc('font', family='Sans-Serif')
bodies ={
'Mars': {'r': 141/93, 'period':687./365, 'phi0':0, 'e':0.093,'plot':{'color':'r'}},
#'Venus': {'r': 67/93, 'o': 365/225, 'phi0':0, 'e':0.007,'plot':{'color':'m'}},
}
def RK4(f):
return lambda t, y, dt: (
lambda dy1: (
lambda dy2: (
lambda dy3: (
lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6
)( dt * f( t + dt , y + dy3 ) )
)( dt * f( t + dt/2, y + dy2/2 ) )
)( dt * f( t + dt/2, y + dy1/2 ) )
)( dt * f( t , y ) )
years =1.1
# d\theta/dt = K*[(1+ e cos(theta)) /(1-e^2)]^2
def kepler_v(period,e):
return lambda t, theta: 2*np.pi/period*np.sqrt(1-e*e)*((1+e*np.cos(theta))/(1-e*e))**2
def kepler_point(r,e,theta,xoff=0,yoff=0):
c = np.cos(theta)
rs = r*(1-e*e)/(1+e*c)
x = rs*c+xoff
y= rs*np.sin(theta)+yoff
return (x, y)
#theta is the angle at the equant (-e,0)
def ptolemy_point(r,e, theta, xoff=0, yoff=0):
c = np.cos(theta)
s = np.sin(theta)
sx = e * s
cx = np.sqrt(1-sx*sx)
rr = r*(e*c + cx)
x = rr*c-2*r*e+xoff
y= rr *s+yoff
return (x,y)
def kepler_seq (period, r, e, revolutions, tstep):
kv = kepler_v(period, e)
kvr4 = RK4(kv)
n = int(revolutions*2*np.pi/tstep)+1
th=0
i=0
t=0
xarr=np.empty([n,])
yarr=np.empty([n,])
xarr[0]=r*(1-e)
yarr[0]=0
while i < n-1:
i = i+1
t=t+tstep
th = th + kvr4(t, th,tstep)
x,y = kepler_point(r,e,th)
xarr[i]=x
yarr[i]=y
return (xarr, yarr)
def ptolemy_seq (period, r, e, revolutions, tstep):
omega = 2*np.pi/period
angles =np.arange(0,2*np.pi*revolutions,tstep) *omega
return ptolemy_point(r, e, angles, 0)
plt.gcf().set_size_inches(8,8)
frames =600
tstep =0.005
posk ={}
posp={}
linek = {}
linep={}
linea={}
lineb={}
linec={}
txt = None
n = int(years*2*np.pi/tstep)
ax = plt.subplot(111)
ax.set_aspect('equal')
fig = plt.gcf()
for b in bodies.keys():
r, period, e = [bodies[b][x] for x in ['r','period','e']]
ax.set_xlim(-r*(1+e)-0.2,r*(1-e)+0.2)
ax.set_ylim(-r-0.1,r+0.1)
posk[b] = kepler_seq(period, r, e, years, tstep)
posp[b]= ptolemy_seq(period, r, e, years, tstep)
xoff = 2*e*r
angle_k = np.arctan2(posk[b][1], posk[b][0] )
angle_p = np.arctan2(posp[b][1], posp[b][0] +xoff)
angle_kp = np.arctan2(posk[b][1], posk[b][0]+xoff)
angle_delta = angle_kp - angle_p
ax.plot([-r*e], [0], marker='o', color='b')
ax.plot([-2*r*e], [0], marker='o', color='g')
ax.annotate('EQT', xy=(-2*r*e,0), xytext=(-2*e*r-0.1, -0.1))
linek[b], =ax.plot(posk[b][0][0],posk[b][1][0], **bodies[b]['plot'], ls='-', markersize=3)
linea[b] = mlines.Line2D([0,posk[b][0][0]],[0, posk[b][1][0]], color='r')
ax.add_line(linea[b])
linep[b], =ax.plot(posp[b][0][0],posp[b][1][0], color="b", ls='-', markersize=3)
lineb[b] = mlines.Line2D([0,posp[b][0][0]],[0, posp[b][1][0]], color='b')
ax.add_line(lineb[b])
linec[b] = mlines.Line2D([-2*r*e,posp[b][0][0]],[0, posp[b][1][0]], color='b',linestyle='dashed',lw=1)
ax.add_line(linec[b])
txt = ax.text(0, -1, "")
#fig2, ax2= plt.subplots()
#ax2.plot(angle_k, angle_delta)
ax.plot([0], [0], marker='o', color='y',markersize=13)
ax.set_title("Ptolemy v. Kepler (e="+str(e)+")")
ax.relim()
ax.autoscale_view(True,True,True)
def init():
return list( [lineX[x] for lineX in [linea, lineb, linec, linek, linep] for x in linek.keys()])+ [txt]
def update(i):
ax.autoscale_view(True,True,True)
for b in bodies.keys():
r, period, e = [bodies[b][x] for x in ['r','period','e']]
n2 = int(n*i/frames)
n1 = 0
if(n2>1200):
n1= n2-1200
#print (i, frames, n, n2)
linea[b].set_data([0,posk[b][0][n2]],[0, posk[b][1][n2]])
lineb[b].set_data([0,posp[b][0][n2]],[0, posp[b][1][n2]])
angles = np.arctan2([posk[b][1][n2], posp[b][1][n2]], [posk[b][0][n2], posp[b][0][n2]])
xoff = 2*e*r
angles2 = np.arctan2([posk[b][1][n2], posp[b][1][n2]], [posk[b][0][n2]+xoff, posp[b][0][n2]+xoff])
delta = (angles[1]-angles[0])*180/np.pi
delta2 = (angles2[1]-angles2[0])*180/np.pi
print ("angle=", angles, delta, angles2, delta2)
txtstr = '$\Delta=%.2f^\circ$' % (delta)
txt.set_text(txtstr)
linec[b].set_data([-2*e*r,posp[b][0][n2]],[0, posp[b][1][n2]])
linek[b].set_data(posk[b][0][n1:n2],posk[b][1][n1:n2])
linep[b].set_data(posp[b][0][n1:n2],posp[b][1][n1:n2])
return list( [lineX[x] for lineX in [linea, lineb, linec, linek, linep] for x in linek.keys()])+ [txt]
animate = animation.FuncAnimation(plt.gcf(), update, init_func=init, frames=frames, interval=100, blit=True)
#animate.save("c:/tmp/kepler-ptolemy-compare"+str(e)+".gif",writer='imagemagick', savefig_kwargs =dict(pad_inches=0))
#animate.save("c:/tmp/kepler-ptolemy-compare"+str(e)+".mp4",writer='ffmpeg', savefig_kwargs =dict(pad_inches=0))
plt.show()