Объект numpy.float64 не может быть интерпретирован как целое число, но я не могу определить, где он запрашивает целое число

Question

Я сталкиваюсь с этим сообщением об ошибке 'numpy.float64' object cannot be interpreted as an integerпри запуске следующего кода:

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. Ошибка происходит именно на

    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' в моем коде. Есть идеи, в чем может быть проблема? Очень ценю вашу помощь!

Answer 1

Проблема возникает потому, что у вас есть M=[], что делает его list, Это полностью отменяет M = m1 + m2 Вы сделали на вершине. Затем на линии 55 у вас есть M * BLAH, где M НЕ является пустым массивом (это список), следовательно, * Оператор интерпретируется как "Повторите список раз BLAH". Вот в чем проблема.

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 пакет, он делает некоторые вещи, которые вы ищете.