求水星近日点轨道进动模拟轨道代码，关键是用python编写的

尼奇怪

RT
这是咱本周计算物理的作业，本来上周的FFT（快速傅立叶变换）就够痛苦了……
以下是我用kepler轨道改编的，不过怎么画出来都觉得不像进动轨道……

1. from pylab import *
   
2. import matplotlib as mpl
   
3. from math import sqrt, pi
   
4. from mpl_toolkits.mplot3d import axes3d, Axes3D
   
5. import numpy as np
   
6. import matplotlib.pyplot as plt
   
7. 
   
8. 
   
9. class Mecury():
   
10. 
    
11.     def __init__(self, initx=0.47, inity=0, initvx=0, initvy=8.2, dt=0.001, powerlaw=2., alpha=0.008, outputfile=''):
    
12.         self.dt=dt
    
13.         self.beta=powerlaw
    
14.         self.alpha=alpha
    
15.         self.filename=outputfile
    
16.         self.calculate(sqrt(initx*initx+inity*inity), initx, inity, initvx, initvy, 0.)
    
17.         self.dth=0.
    
18. 
    
19.     def stepOne(self, r, x, y, vx, vy, t):
    
20.         r=sqrt(x*x+y*y)
    
21.         vx=vx-(4*pi*pi*x*self.dt)/r**(self.beta+1)/(1+self.alpha/r**(self.beta))
    
22.         vy=vy-(4*pi*pi*y*self.dt)/r**(self.beta+1)/(1+self.alpha/r**(self.beta))
    
23.         x=x+vx*self.dt
    
24.         y=y+vy*self.dt
    
25.         t=t+self.dt
    
26.         return r, x, y, vx, vy, t
    
27. 
    
28.     def calculate(self, r, x ,y ,vx, vy, t):
    
29.         self.data=[(r, x,y,vx,vy,t)]
    
30.         while t<100:
    
31.             r, x, y, vx, vy, t=self.stepOne(r, x ,y ,vx, vy, t)
    
32.             self.data.append((r, x,y,vx,vy,t))
    
33.         if self.filename!='':
    
34.             self.store()
    
35. 
    
36. 
    
37. 
    
38.     def store(self):
    
39.         pass
    
40. 
    
41.     def plotMecury(self, figType=1):
    
42.         x=[]; y=[]; z=[0.]; vx=[]; vy=[]; t=[]; r=[]; xp=[]; yp=[]; k=[]
    
43.         for i in range(5000):
    
44.             (r1, x1, y1, vx1, vy1, t1)=self.data[i]
    
45.             x.append(x1); y.append(y1); vx.append(vx1); vy.append(vy1)
    
46.             t.append(t1); z.append(0.); r.append(r1)
    
47.             if i>2:
    
48.                 dr1=r[i-1]-r[i-2]
    
49.                 dr2=r[i]-r[i-1]
    
50.                 if dr1>0 and dr2<0:
    
51.                     xp.append(x[i-1])
    
52.                     yp.append(y[i-1])
    
53.             m=len(self.data)
    
54.         if figType==1:
    
55.             plot(x,y)
    
56.             plot(0,0,'r*')
    
57.             for j in range(len(xp)):
    
58.                 plot([0,xp[j]], [0,yp[j]],'k-')
    
59.             xlim(-0.55,0.55)
    
60.             ylim(-0.55,0.55)
    
61.             show()
    
62.         else:
    
63.             fig=plt.figure()
    
64.             ax=Axes3D(fig)
    
65.             ax.plot(x[0:m], y[0:m], z[0:m], label='Kepler')
    
66.             ax.legend()
    
67.             plt.show()
    
68.             
    
69.     def plotPrecessionrate(self, figType=1):
    
70.         x=[]; y=[]; z=[0.]; vx=[]; vy=[]; t=[]; r=[]; xp=[]; yp=[]; tp=[]; k=[];
    
71.         theta=[]; th=0.; dth=0.#dth=[]
    
72.         for i in range(3000):
    
73.             (r1, x1, y1, vx1, vy1, t1)=self.data[i]
    
74.             x.append(x1); y.append(y1); vx.append(vx1); vy.append(vy1)
    
75.             t.append(t1); z.append(0.); r.append(r1)
    
76.             if i>2:
    
77.                 dr1=r[i-1]-r[i-2]
    
78.                 dr2=r[i]-r[i-1]
    
79.                 if dr1>0 and dr2<0:
    
80.                     xp.append(x[i-1])
    
81.                     yp.append(y[i-1])
    
82.                     tp.append(t[i-1]-0.24500000000000019)
    
83.             m=len(self.data)
    
84.             n=len(xp)
    
85.         if figType==1:
    
86.             for j in range(n):
    
87.                 k.append(yp[j]/xp[j])
    
88.                 theta.append((180/pi)*(arctan((k[0]-k[j])/(1+k[0]*k[j]))))
    
89.             for l in range(n-1):
    
90.                 th=th+theta[l+1]/tp[l+1]
    
91.             dth=th/n
    
92.             print theta
    
93.             print tp
    
94.             print dth
    
95. 
    
96.             
    
97.                
    
98.             #print th
    
99.             #print theta[0]/tp[0]
    
100.             plot(tp,theta)
     
101.             #plot(0,0,'r*')
     
102.             xlim(0,3)
     
103.             ylim(0,20)
     
104.             show()
     
105.             
     
106.         else:
     
107.             fig=plt.figure()
     
108.             ax=Axes3D(fig)
     
109.             ax.plot(x[0:m], y[0:m], z[0:m], label='Kepler')
     
110.             ax.legend()
     
111.             plt.show()
     
112.         #return dth
     
113.    
     
114. 
     
115. 
     
116.     def calculatedth(self, figType=1):
     
117.         x=[]; y=[]; z=[0.]; vx=[]; vy=[]; t=[]; r=[]; xp=[]; yp=[]; tp=[]; k=[];
     
118.         theta=[]; th=0.; dth=0.#dth=[]
     
119.         for i in range(3000):
     
120.             (r1, x1, y1, vx1, vy1, t1)=self.data[i]
     
121.             x.append(x1); y.append(y1); vx.append(vx1); vy.append(vy1)
     
122.             t.append(t1); z.append(0.); r.append(r1)
     
123.             if i>2:
     
124.                 dr1=r[i-1]-r[i-2]
     
125.                 dr2=r[i]-r[i-1]
     
126.                 if dr1>0 and dr2<0:
     
127.                     xp.append(x[i-1])
     
128.                     yp.append(y[i-1])
     
129.                     tp.append(t[i-1]-0.24500000000000019)
     
130.             m=len(self.data)
     
131.             n=len(xp)
     
132.         if figType==1:
     
133.             for j in range(n):
     
134.                 k.append(yp[j]/xp[j])
     
135.                 theta.append((180/pi)*(arctan((k[0]-k[j])/(1+k[0]*k[j]))))
     
136.             for l in range(n-1):
     
137.                 th=th+theta[l+1]/tp[l+1]
     
138.             self.dth=th/n
     
139. 
     
140. def calculaterate():
     
141.     listalpha=[]
     
142.     listdth=[]
     
143.     list1=[0,1,2,4,8,16]
     
144.     rate=0.
     
145.     for k in range(6):
     
146.         list1[k]=0.0002*list1[k]
     
147.         c=Mecury(alpha=list1[k])
     
148.         c.calculatedth()
     
149.         listalpha.append(c.alpha)
     
150.         #print c.alpha
     
151.         listdth.append(c.dth)
     
152.         #print c.dth
     
153.     print listalpha,listdth   
     
154.     plot(listalpha, listdth, 'ro')
     
155.     for m in range(5):
     
156.         rate+=listdth[m+1]/listalpha[m+1]
     
157.     drate=rate/6
     
158.     print drate
     
159.     show()
     
160. 
     
161. a=Mecury(alpha=0.01)
     
162. a.plotMecury()
     
163. #b=Mecury(alpha=0.0008)
     
164. #b.plotPrecessionrate()
     
165. 
     
166. #calculaterate()
     
167. 
     

复制代码

最终成品：

图中黑线为太阳与远日点连线
看起来好圆，好假，而且运算出来的进动角速度和理论值差了25%……

bxfsquare

[S::funk:]没仔细看代码,好多4*pi*pi  
看了半天说明,我也没看懂,下面就当我胡说把
R要是有公式可以得N步以后的R
也许你可以把每圈的r分成2区间,然后各取一个最有代表性的R再按条件比较,得R,如此反复迭代得R(x,y)
就是不断划分R,去逼近远日点
另外水星轨道每一圈的世界坐标系旋转后好象都是一样的?也许计算机图形学里面的坐标旋转Matrix可以用在这里减少计算量

大闸蟹

嗯，各种看不懂

4qwer

你们再说些什么
完全没有看懂啊

lipancs

我连基本的经典物理学直线运动相关方程都想不起来了

黑暗枪骑兵

尼奇怪 发表于 2011-4-26 23:10
回复 黑暗枪骑兵 的帖子

我需要的效果不是那么简单的

我擦。。。我完全理解错了= =

尼奇怪

2qwer 发表于 2011-4-27 18:34
那你為何不用大小姐的方法來取r的極大值後才代入公式中計算呢？

喵~

是这个样子的，x,y为自变量，r为x，y所决定的一个因变量（在这里就是r为x,y的平方和），所以不是简单的说知道了r的极值后代入算式就能获得对应的x，y值（满足条件的x，y可有无穷组）

呃……其实我已经做完了，不过现在觉得原来解决这个问题的办法十分蛋疼（x,y,r各使用了一个长度为n的数组保存所有数据，然后把r极值处的元素序号对应的x，y抽出来）所以想知道有没有更简洁的思路

2qwer

尼奇怪 发表于 2011-4-27 12:10
回复 黑暗枪骑兵 的帖子

我需要的效果不是那么简单的

那你為何不用大小姐的方法來取r的極大值後才代入公式中計算呢？

喵~
p.s.其實一點也不懂的路過的說~~~

尼奇怪

回复 黑暗枪骑兵 的帖子

我需要的效果不是那么简单的
1.我要找的是r的极大值，而不是最大值，就是说在一定区间内的最大值，而非全区间内最大值
2.我需要的不是r的极值，而是对应于r出现极大值时的x，y的值

黑暗枪骑兵

尼奇怪 发表于 2011-4-26 08:07
回复 黑暗枪骑兵 的帖子

询问一下：

你把那些数值放入一个Array里，然后用max(array名）来得到最大值