求水星近日点轨道进动模拟轨道代码，关键是用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%……

0809000

咱上周的作业是用Java Jsp MySQL......
写个选课系统~~~~
蛋疼~~~
我们才讲到JSP就要交作业~~其实课上什么也没讲~~
下次的是用Java写个多人聊天系统~~
现在天朝大学都是这么折腾人的吗~~

输代码要点开那个<>标志的之后输入就行了~

stargazer

tx们都是学计算机的么？可惜我只看过c。

water_long

總覺得用netbean弄個視窗是比用visual studio的幸福

尼奇怪

回复 stargazer 的帖子

Python这东西比C可读性要强，但是格式变态

计算机对于理工科学生是技能要求，而不是知识要求，所以很多非计算机专业也要学习
虽然计算物理也许不是我的研究方向，不过作业总是得做的

其实我倒不觉得这是折腾人，只不过中国学生长期抄过来的，遇到一点要自己独创的东西就容易卡壳，我也一样
-----------------------------------------------------------------
以下是实验结果：

步数次数较少时

步数较多时

三个远日点

更多……这回终于是圆的了
找到远日点，如此一看，似乎还有点像进动轨道了
虽然悲剧的宽屏把正方形压扁了导致轨道看起来是椭圆……

0809000

回复 尼奇怪 的帖子

还真是抄习惯了~~~
不过抄有时也是很麻烦的~~~
你们学着是业余~~

Python很不错的~~用于计算~~格式更变态的JCL也学过~~~Python不限数字格式已经不错了~~

我们可就是直接指一行代码~~这啥意思~~
在就是~~某某功能的代码在那里~~

黑暗枪骑兵

Python看得懂。。但这物理没学过= =

skulldownz

不会看,学,话说 pass 是什么意思

进动轨道什么的也不知道....
不过找了下: 这个会不会有用...虽然只是理论而不是程序...不过我连理论都没听过撒..哈哈哈
http://tieba.baidu.com/f?kz=595300156

skulldownz

喔..又学到一点东西了..哈哈

这么短时间理论是搞不灵的了,但实际的图是类似这样的吗? (网上的图片)


http://zh.wikipedia.org/wiki/广义相对论

恩..不知道网上这些图是否有夸张过...不过如果确实如这些图一般的话,只能看看有没有比例错的地方撒....

stargazer

曾经中学时代对物理尤其是天体物理很感兴趣，成绩也不错，可惜当时数学差，大学后和物理无缘，基本耽误在电脑上了。倒是现在偶尔抽空学习下计算机的东西。
这个轨道即使是和万有引力相关的方程？俺的物理退化到初中水平了