homework_11 计算程序

计算程序1

import math
G=4*math.pi**2
m1=1
m2=100
dt=0.001
end_t=20
t=[0]
x_1=[10]
v_x_1=[0]
y_1=[0]
v_y_1=[(4.0*math.pi**2*10/(10.1**2))**0.5]
r_1=[]
theta_1=[]
x_2=[-0.1]
v_x_2=[0]
y_2=[0]
v_y_2=[-(0.1*4*math.pi**2/(10.1**2))**0.5]
r_2=[]
theta_2=[]
r=[]
n_i=[]
n_j=[]
for i in range(int(end_t/dt)):
      r_1.append((x_1[i]**2+y_1[i]**2)**0.5)
      if y_1[i]<0:
         theta_1.append(2*math.pi-math.acos(x_1[i]/(x_1[i]**2+y_1[i]**2)**0.5))
      else:
         theta_1.append(math.acos(x_1[i]/(x_1[i]**2+y_1[i]**2)**0.5))
      r_2.append((x_2[i]**2+y_2[i]**2)**0.5)
      if y_2[i]<0:
         theta_2.append(2*math.pi-math.acos(x_2[i]/(x_2[i]**2+y_2[i]**2)**0.5))
      else:
         theta_2.append(math.acos(x_2[i]/(x_2[i]**2+y_2[i]**2)**0.5))
      #print theta_1[i]
      #print theta_2[i]
      r.append((r_1[i]**2+r_2[i]**2-2*math.cos(theta_1[i]-theta_2[i])*r_1[i]*r_2[i])**0.5)
      #print r[-1]
      n_i.append((r_2[i]*math.cos(theta_2[i])-r_1[i]*math.cos(theta_1[i]))/((r_2[i]*math.cos(theta_2[i])-r_1[i]*math.cos(theta_1[i]))**2+(r_2[i]*math.sin(theta_2[i])-r_1[i]*math.sin(theta_1[i]))**2)**0.5)
      n_j.append((r_2[i]*math.sin(theta_2[i])-r_1[i]*math.sin(theta_1[i]))/((r_2[i]*math.cos(theta_2[i])-r_1[i]*math.cos(theta_1[i]))**2+(r_2[i]*math.sin(theta_2[i])-r_1[i]*math.sin(theta_1[i]))**2)**0.5)
      print n_j[i]
      v_x_1.append(v_x_1[i]+(G*m2/r[i]**2)*n_i[i]*dt)
      x_1.append(x_1[i]+v_x_1[i+1]*dt)
      v_y_1.append(v_y_1[i]+(G*m2/r[i]**2)*n_j[i]*dt)
      y_1.append(y_1[i]+v_y_1[i+1]*dt)
      v_x_2.append(v_x_2[i]-(G*m1/r[i]**2)*n_i[i]*dt)
      x_2.append(x_2[i]+v_x_2[i+1]*dt)
      v_y_2.append(v_y_2[i]-(G*m1/r[i]**2)*n_j[i]*dt)
      y_2.append(y_2[i]+v_y_2[i+1]*dt)
      t.append(t[i]+dt)
import matplotlib.pyplot as plt
import numpy as np
plt.plot(x_1,y_1,'--' ,label='star_1,$r(2,0),v(0,4/3 \pi-0.1,m_1=1$')
plt.plot(x_2,y_2,'--',label='star_2,$r(-1,0),v(0,-2/3 \pi+0.1,m_2=2$')
plt.title(u'the orbit of star',fontsize=14)
plt.xlabel(u'x/AU',fontsize=14)
plt.ylabel(u'y/AU',fontsize=14)
plt.legend(fontsize=12,loc='best')
plt.show()

程序2

import math
G=4*math.pi**2
m1=1
m2=1
m3=1
dt=0.001
end_t=12
t=[0]
x_1=[0]
v_x_1=[-((4*3**0.5*math.pi**2)/3.0)**0.5]
y_1=[1]
v_y_1=[0]
r_1=[]
theta_1=[]
x_2=[-3**0.5/2.0]
v_x_2=[0.5*((4*3**0.5*math.pi**2)/3.0)**0.5]
y_2=[-1/2.0]
v_y_2=[-(3**0.5/2.0)*((4*3**0.5*math.pi**2)/3.0)**0.5]
r_2=[]
theta_2=[]
x_3=[3.0**0.5/2]
v_x_3=[0.5*((4*3**0.5*math.pi**2)/3.0)**0.5]
y_3=[-1.0/2]
v_y_3=[3**0.5/2*((4*3**0.5*math.pi**2)/3.0)**0.5]
r_3=[]
theta_3=[]
r12=[]
r13=[]
r23=[]
n12_i=[]
n13_i=[]
n23_i=[]
n12_j=[]
n13_j=[]
n23_j=[]
for i in range(int(end_t/dt)):
      ''' 
      r_1.append((x_1[i]**2+y_1[i]**2)**0.5)
      if y_1[i]<0:
         theta_1.append(2*math.pi-math.acos(x_1[i]/(x_1[i]**2+y_1[i]**2)**0.5))
      else:
         theta_1.append(math.acos(x_1[i]/(x_1[i]**2+y_1[i]**2)**0.5))
      r_2.append((x_2[i]**2+y_2[i]**2)**0.5)
      if y_2[i]<0:
         theta_2.append(2*math.pi-math.acos(x_2[i]/(x_2[i]**2+y_2[i]**2)**0.5))
      else:
         theta_2.append(math.acos(x_2[i]/(x_2[i]**2+y_2[i]**2)**0.5))
      r_3.append((x_3[i]**2+y_3[i]**2)**0.5)
      if y_1[i]<0:
         theta_3.append(2*math.pi-math.acos(x_3[i]/(x_3[i]**2+y_3[i]**2)**0.5))
      else:
         theta_3.append(math.acos(x_3[i]/(x_3[i]**2+y_3[i]**2)**0.5))
      #print theta_3[i]
      #print theta_2[i]
      r12.append((r_1[i]**2+r_2[i]**2-2*math.cos(theta_1[i]-theta_2[i])*r_1[i]*r_2[i])**0.5)
      r13.append((r_1[i]**2+r_3[i]**2-2*math.cos(theta_1[i]-theta_3[i])*r_1[i]*r_3[i])**0.5)
      r23.append((r_3[i]**2+r_2[i]**2-2*math.cos(theta_3[i]-theta_2[i])*r_3[i]*r_2[i])**0.5)
      print r_1[-1]
      n12_i.append((r_2[i]*math.cos(theta_2[i])-r_1[i]*math.cos(theta_1[i]))/((r_2[i]*math.cos(theta_2[i])-r_1[i]*math.cos(theta_1[i]))**2+(r_2[i]*math.sin(theta_2[i])-r_1[i]*math.sin(theta_1[i]))**2)**0.5)
      n13_i.append((r_3[i]*math.cos(theta_3[i])-r_1[i]*math.cos(theta_1[i]))/((r_3[i]*math.cos(theta_3[i])-r_1[i]*math.cos(theta_1[i]))**2+(r_3[i]*math.sin(theta_3[i])-r_1[i]*math.sin(theta_1[i]))**2)**0.5)      
      n23_i.append((r_2[i]*math.cos(theta_2[i])-r_3[i]*math.cos(theta_3[i]))/((r_2[i]*math.cos(theta_2[i])-r_3[i]*math.cos(theta_3[i]))**2+(r_2[i]*math.sin(theta_2[i])-r_3[i]*math.sin(theta_3[i]))**2)**0.5)
      n12_j.append((r_2[i]*math.sin(theta_2[i])-r_1[i]*math.sin(theta_1[i]))/((r_2[i]*math.cos(theta_2[i])-r_1[i]*math.cos(theta_1[i]))**2+(r_2[i]*math.sin(theta_2[i])-r_1[i]*math.sin(theta_1[i]))**2)**0.5)
      n13_j.append((r_3[i]*math.sin(theta_3[i])-r_1[i]*math.sin(theta_1[i]))/((r_3[i]*math.cos(theta_3[i])-r_1[i]*math.cos(theta_1[i]))**2+(r_3[i]*math.sin(theta_3[i])-r_1[i]*math.sin(theta_1[i]))**2)**0.5)  
      n23_j.append((r_2[i]*math.sin(theta_2[i])-r_3[i]*math.sin(theta_3[i]))/((r_2[i]*math.cos(theta_2[i])-r_3[i]*math.cos(theta_3[i]))**2+(r_2[i]*math.sin(theta_2[i])-r_3[i]*math.sin(theta_3[i]))**2)**0.5)    
      #print n_j[i]
      v_x_1.append(v_x_1[i]+(G*m2/r12[i]**2)*n12_i[i]*dt+(G*m3/r13[i]**2)*n13_i[i]*dt)
      x_1.append(x_1[i]+v_x_1[i+1]*dt)
      v_y_1.append(v_y_1[i]+(G*m2/r12[i]**2)*n12_j[i]*dt+(G*m3/r13[i]**2)*n13_j[i]*dt)
      y_1.append(y_1[i]+v_y_1[i+1]*dt)
      v_x_2.append(v_x_2[i]-(G*m1/r12[i]**2)*n12_i[i]*dt-(G*m3/r23[i]**2)*n23_i[i]*dt)
      x_2.append(x_2[i]+v_x_2[i+1]*dt)
      v_y_2.append(v_y_2[i]-(G*m1/r12[i]**2)*n12_j[i]*dt-(G*m3/r23[i]**2)*n23_j[i]*dt)
      y_2.append(y_2[i]+v_y_2[i+1]*dt)
      v_x_3.append(v_x_3[i]-(G*m1/r13[i]**2)*n13_i[i]*dt+(G*m2/r23[i]**2)*n23_i[i]*dt)
      x_3.append(x_3[i]+v_x_3[i+1]*dt)
      v_y_3.append(v_y_3[i]-(G*m1/r13[i]**2)*n13_j[i]*dt+(G*m2/r23[i]**2)*n23_j[i]*dt)
      y_3.append(y_3[i]+v_y_3[i+1]*dt)
      t.append(t[i]+dt)
      '''
      v_x_1.append(v_x_1[i]+(G*m2*(x_2[i]-x_1[i])/((x_1[i]-x_2[i])**2+(y_1[i]-y_2[i])**2)**(1.5)+G*m3*(x_3[i]-x_1[i])/((x_1[i]-x_3[i])**2+(y_1[i]-y_3[i])**2)**(1.5))*dt)
      v_x_2.append(v_x_2[i]+(G*m1*(x_1[i]-x_2[i])/((x_1[i]-x_2[i])**2+(y_1[i]-y_2[i])**2)**(1.5)+G*m3*(x_3[i]-x_2[i])/((x_2[i]-x_3[i])**2+(y_2[i]-y_3[i])**2)**(1.5))*dt)
      v_x_3.append(v_x_3[i]+(G*m1*(x_1[i]-x_3[i])/((x_1[i]-x_3[i])**2+(y_1[i]-y_3[i])**2)**(1.5)+G*m2*(x_2[i]-x_3[i])/((x_2[i]-x_3[i])**2+(y_2[i]-y_3[i])**2)**(1.5))*dt)
      v_y_1.append(v_y_1[i]+(G*m2*(y_2[i]-y_1[i])/((x_1[i]-x_2[i])**2+(y_1[i]-y_2[i])**2)**(1.5)+G*m3*(y_3[i]-y_1[i])/((x_1[i]-x_3[i])**2+(y_1[i]-y_3[i])**2)**(1.5))*dt)
      v_y_2.append(v_y_2[i]+(G*m1*(y_1[i]-y_2[i])/((x_1[i]-x_2[i])**2+(y_1[i]-y_2[i])**2)**(1.5)+G*m3*(y_3[i]-y_2[i])/((x_2[i]-x_3[i])**2+(y_2[i]-y_3[i])**2)**(1.5))*dt)
      v_y_3.append(v_y_3[i]+(G*m1*(y_1[i]-y_3[i])/((x_1[i]-x_3[i])**2+(y_1[i]-y_3[i])**2)**(1.5)+G*m2*(y_2[i]-y_3[i])/((x_2[i]-x_3[i])**2+(y_2[i]-y_3[i])**2)**(1.5))*dt)
      x_1.append(x_1[i]+v_x_1[i+1]*dt)
      x_2.append(x_2[i]+v_x_2[i+1]*dt)
      x_3.append(x_3[i]+v_x_3[i+1]*dt)
      y_1.append(y_1[i]+v_y_1[i+1]*dt)
      y_2.append(y_2[i]+v_y_2[i+1]*dt)
      y_3.append(y_3[i]+v_y_3[i+1]*dt)
import matplotlib.pyplot as plt
import numpy as np
plt.plot(x_1,y_1,'--' ,label='star_1,$r_1=1,m_1=1$')
plt.plot(x_2,y_2,'--',label='star_2,$r_2=1,m_2=1$')
plt.plot(x_3,y_3,'--',label='star_2,$r_3=1,m_3=1$')
plt.title(u'the orbit of star',fontsize=14)
plt.xlabel(u'x/AU',fontsize=14)
plt.ylabel(u'y/AU',fontsize=14)
plt.legend(fontsize=12,loc='best')
plt.show()