首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >“numpy.float64”对象不能解释为整数,但我无法确定它在哪里请求整数

“numpy.float64”对象不能解释为整数,但我无法确定它在哪里请求整数
EN

Stack Overflow用户
提问于 2018-05-30 10:27:17
回答 1查看 8.6K关注 0票数 0

在运行以下代码时,我遇到此错误消息'numpy.float64' object cannot be interpreted as an integer

代码语言:javascript
复制
from scipy.integrate import nquad
from numpy import sqrt, sin, cos, pi, arcsin, maximum, log, log10, exp
from matplotlib.pyplot import xscale, yscale, plot, show, savefig

#define constants
msun=1500
m1=30*msun
m2=30*msun
M=m1+m2
eta=m1*m2/M**2
mSMBH=4.2*10**6*msun
v=5/10**4
pc=3.09*10**16
r0=1.005*10**21              #nfw parameter#
rho0=1.27/10**49             #nfw parameter#
pmtpy=3*10**8*365*24*3600    #rate from per-meter to per-yr
rhoc=1.35*10**(-53)
cm=2.16

#define functions
def g(x):
    return log(1+x)-x/(1+x)
def y(x):
    return 1.43*10**12*1500/x
def sig(x):
    return 13.216*y(x)**0.41/(1+1.102*y(x)**0.2+6.22*y(x)**0.333)
def c(x):
    return 2.881*(1+(sig(x)/1.257)**1.022)*exp(0.06/sig(x)**2)
def r0(x):
    return (3*x/(800*c(x)**3*rhoc*pi))**(1/3)
def rho0(x):
    return 200*rhoc*c(x)**3/3/g(c(x))
def sigma(x):
    return sqrt(4*pi*rho0(x)*r0(x)**2*g(cm)/cm)
def mSMBH(x):
    return 1500*10**(8.15+4.38*log10(1500*sigma(x)))
def rcap(x):
    return 8*mSMBH(x)

def Rsp(x):
    return 0.122*sqrt(mSMBH(x)/rho0(x)/r0(x))
def rhosp(r,x):
    return 3*rho0(x)*r0(x)*Rsp(x)**(4/3)*(1-rcap(x)/r)**3/r**(7/3)

#x,y,z are to be integrated, and m is the x-axis parameter to be searched through using while loop
def v(z,m):
    return sqrt(mSMBH(m)/z)
def vrel(y,z,m):
    return 2*v(z,m)*sin(y/2)
def vcm(y,z,m):
    return 2*v(z,m)*cos(y/2)
def bmin(y,z,m):
    return 2*M/vrel(y,z,m)
def bmax(y,z,m):
    return maximum(bmin(y,z,m), M*(340*eta/(3*cos(1)*vrel(y,z,m)**9))**(1/7))
def a(x,y,z,m):
    if bmax(y,z,m)>bmin(y,z,m):
        return M/(vrel(y,z,m)**2*((bmax(y,z,m)/x)**7-1))
    else:
        return 0
def rt(x,y,z,m):
    return a(x,y,z,m)*(mSMBH(m)/M)**(1/3)
def rd(x,y,z,m):
    return maximum(rcap(m),rt(x,y,z,m))
def rz(x,y,z,m):
    return rd(x,y,z,m)*(sqrt(1+8*mSMBH(m)/(rd(x,y,z,m)*vcm(y,z,m)**2))-1)/2
def thetal(x,y,z,m):
    if rd(x,y,z,m)<z and rz(x,y,z,m)<z:
        return arcsin(rd(x,y,z,m)*sqrt(2*(vcm(y,z,m)**2/2-mSMBH(m)/z+mSMBH(m)/rd(x,y,z,m)))/(vcm(y,z,m)*z))
    else:
        return pi/2

#prepare looping tools    
m=10**9*msun
M=[]           #x value array
Rate=[]        #result estimation array
RateUp=[]      #result upper bound array
RateDown=[]    #result lower bound array
options={'limit':200}

while m<10**14*msun:
    def zbound():
        return [rcap(m),Rsp(m)]
    def ybound(z_foo):
        return [10**(-10),pi-10**(-10)]
    def xbound(y,z):
        return [bmin(y,z,m)-10**(-6),bmax(y,z,m)-10**(-6)]
    def f(x,y,z):
        return pmtpy*8*pi**2*v(z,m)/m1**2*rhosp(z,m)**2*z**2*sin(y/2)**2*cos(y/2)*x*cos(thetal(x,y,z,m))
    ans=nquad(f, [xbound,ybound,zbound],opts=[options,options,options])   
    M.append(m)
    Rate.append(ans[0])
    RateUp.append(ans[0]+ans[1])
    RateDown.append(ans[0]-ans[1])
    m+=10**15

xscale('log')
yscale('log')
plot(M,Rate,M,RateUp,M,RateDown)
show() 

基本上,我试图求一个依赖于参数m的三重积分。

代码语言:javascript
复制
    File "C:/Users/user/.spyder-py3/temp.py", line 83, in xbound
        return [bmin(y,z,m)-10**(-6),bmax(y,z,m)-10**(-6)]
    File "C:/Users/user/.spyder-py3/temp.py", line 52, in bmax
        return maximum(bmin(y,z,m), M*(340*eta/(3*cos(1)*vrel(y,z,m)**9))**(1/7))

TypeError: 'numpy.float64' object cannot be interpreted as an integer

但是我真的不知道为什么在我的代码中“bmax”中的任何地方都需要一个整数。知道有什么问题吗?真的很感谢你的帮助!

EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2018-05-30 10:56:38

问题的出现是因为您有M=[],这使它成为一个list。这完全超越了您在顶部所做的M = m1 + m2。然后在第55行,有M * BLAH,其中M不是一个numpy数组(它是一个列表),因此*操作符被解释为“重复列表的时间”。这就是问题所在。

PS:将其更改为np.array(M) * (340*eta/(3*cos(1.)*vrel(y,z,m)**9))**(1/7)解决了这个问题,但创建了其他问题。例如,如果您使用的是python2,那么(1/7) = 0,您可能需要(1./7.),等等。您可以预先了解大量信息,因此使用np.arange创建完整的M数组并进行处理可能会更快。

PS:既然你有天文学相关的计算,请查看astropy软件包,它会做一些你正在寻找的事情。

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/50602337

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档