2013-04-22 116 views
-1

原来的二阶微分方程是整合轨道轨迹2

x'' - 2 * omega * y' - omega ** 2 * x = - mue * (x + pi2 * r12)/np.sqrt((x + pi2 * r12) ** 2 + y ** 2) ** 3 - mum * (x - pi1 * r12)/np.sqrt((x - pi1 * r12) ** 2 + y ** 2) 
y'' + 2 * omega * x' - omega **2 * y = - mue * y/np.sqrt((x + pi2 * r12) ** 2 + y ** 2) ** 3 - mum * y/np.sqrt((x - pi1 * r12) ** 2 + y ** 2) 
z'' = 0 

因此,这里是我用来解决ODE但首先我把它分成2个第一批订单的代码。

我收到第61行模块不可调用的错误。

线61 u = odeint(deriv, u0, dt)

#!/usr/bin/env python                

import numpy as np 
import scipy.integrate as odeint 
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D 

me = 5.974 * 10 ** (24) # mass of the earth          
mm = 7.348 * 10 ** (22) # mass of the moon          
G = 6.67259 * 10 ** (-20) # gravitational parameter        
re = 6378.0 # radius of the earth in km           
rm = 1737.0 # radius of the moon in km           
r12 = 384400.0 # distance between the CoM of the earth and moon     
M = me + mm 

pi1 = me/M 
pi2 = mm/M 
mue = 398600.0 # gravitational parameter of earth km^3/sec^2      
mum = G * mm # grav param of the moon           
mu = mue + mum 
omega = np.sqrt(mu/r12 ** 3) 
nu = 0.0 # flight path angle              

x = 327156.0 # x location where the moon's SOI effects the spacecraft   
y = 33050.0 # y location              

vbo = 10.85 # velocity at burnout            

gamma = -141.868 * np.pi/180 # angle in radians of true anomaly    

vx = vbo * (np.sin(gamma) * np.cos(nu) - np.cos(gamma) * np.sin(nu)) 
# velocity of the bo in the x direction           
vy = vbo * (np.sin(gamma) * np.sin(nu) + np.cos(gamma) * np.cos(nu)) 
# velocity of the bo in the y direction           

xrel = (re + 300.0) * np.cos(gamma) 
# spacecraft x location relative to the earth          
yrel = (re + 300.0) * np.sin(gamma) 

# r0 = [xrel, yrel, 0]               
# v0 = [vx, vy, 0]    
u0 = [xrel, yrel, 0, vx, vy, 0] 


def deriv(u, dt): 
    n1 = -((mue * (u[0] + pi2 * r12)/np.sqrt((u[0] + pi2 * r12) ** 2 
               + u[1] ** 2) ** 3) 
     - (mum * (u[0] - pi1 * r12)/np.sqrt((u[0] - pi1 * r12) ** 2 
               + u[1] ** 2) ** 3)) 
    n2 = -((mue * u[1]/np.sqrt((u[0] + pi2 * r12) ** 2 + u[1] ** 2) ** 3) 
     - (mum * u[1]/np.sqrt((u[0] - pi1 * r12) ** 2 + u[1] ** 2) ** 3)) 
    return [u[3], # dotu[0] = u[3]            
      u[4], # dotu[1] = u[4]            
      u[5], # dotu[2] = u[5]            
      2 * omega * u[5] + omega ** 2 * u[0] + n1, # dotu[3] = that   
      omega ** 2 * u[1] - 2 * omega * u[4] + n2, # dotu[4] = that   
      0] # dotu[5] = 0              


dt = np.arange(0.0, 250000.0, .1) 
u = odeint(deriv, u0, dt) 
x, y, z, x2, y2, z2 = u.T 

fig = plt.figure() 
ax = fig.add_subplot(111, projection='3d') 
ax.plot(x, y, z) 
plt.show() 

回答

2

假设你的意思是这个错误:因为你想要的功能scipy.integrate命名为odeint

~/coding$ python orbit1.py 
Traceback (most recent call last): 
    File "orbit1.py", line 61, in <module> 
    u = odeint(deriv, u0, dt) 
TypeError: 'module' object is not callable 

这是。您的线路

import scipy.integrate as odeint 

导入整个模块并给它起个名字odeint。尝试

from scipy.integrate import odeint 

代替,或

import scipy.integrate 
[...] 

u = scipy.integrate.odeint(deriv, u0, dt) 

这应该给你 enter image description here