I am trying to run this code but I get the error in the title. I looked for the documentation and examples about tplquad, but I couldn't understand my problem. Thank you very much in advance!
here my code:
from numpy import *
from pylab import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from math import *
from scipy.integrate import quad,dblquad,tplquad
from scipy.integrate import nquad
fig_width = 6.
fig_height = fig_width*0.75
fig_size = [fig_width,fig_height]
params = {'backend': 'TkAgg',
'axes.labelsize': 30,
'text.fontsize': 20,
'title.fontsize': 20,
'legend.fontsize': 20,
'xtick.labelsize': 20,
'ytick.labelsize': 20,
'text.usetex': False,
'font.family': 'sans-serif',
'figure.figsize': fig_size}
rcParams.update(params)
pi=3.14
pt_T=3.
#T=0.47
thetaP= -pi
precision=5
y=0
M_T=linspace(1.,7.,precision)
integral1d=[0]*precision
#chi now is def with a plus instead of the minus in the article
def chi(thetap1,p1,thetaP,T,M_T):
return abs((2*p1*T*sqrt(pt_T**2+(M_T**2+pt_T**2)*sinh(y)**2)*sin(thetaP)*sin(thetap1))**2 - (2*p1*T*(sqrt(M_T**2+pt_T**2)*cosh(y)- sqrt(pt_T**2+(M_T**2+pt_T**2)*sinh(y)**2)*cos(thetaP) *cos(thetap1) )-(T**2)* M_T**2)**2)+1
def p1max(thetaP, thetap1,T,M_T):
return (M_T**2)*T/(2*(sqrt(M_T**2+pt_T**2)*cosh(y)- sqrt(pt_T**2+(M_T**2+pt_T**2))*sinh(y)**2*cos(thetaP-thetap1)))-0.1
def p1min(thetaP, thetap1,T,M_T):
#NOT SURE ABOUT THE T AT DENOMINATOR
return (M_T**2)*T/(2*(sqrt(M_T**2+pt_T**2)*cosh(y)- sqrt(pt_T**2+(M_T**2+pt_T**2))*sinh(y)**2*cos(thetaP+thetap1))) +0.1
def integral(thetaP,T,M_T):
area =dblquad(lambda p1, thetap1: 5*(1/(18*pi**5))*sin(thetap1)*(p1/(sqrt(chi(thetap1,p1,thetaP,T,M_T))))*(1/(exp(p1/T) + 1))*(1/(exp((sqrt(M_T**2 + pt_T**2)*cosh(y) - p1/T) +1))) , -pi+0.1, -0.1, lambda p1: p1min(thetaP, p1,T,M_T), lambda p1: p1max(thetaP,p1,T,M_T)) #CHANGE 1., lambda p1:10.)
return area[0]
def integrand(M_T, p1,thetap1,T):
return pt_T*T*2*pi*5*(1/(18*pi**5))*sin(thetap1)*(p1/(sqrt(chi(thetap1,p1,thetaP,T,M_T))))*(1/(exp(p1/T) + 1))*(1/(exp((sqrt(M_T**2 + pt_T**2)*cosh(y) - p1/T) +1)))
def formula151(M_T):
area =tplquad(lambda p1, thetap1,T: pt_T*T*2*pi*5*(1/(18*pi**5))*sin(thetap1)*(p1/(sqrt(chi(thetap1,p1,thetaP,T,M_T))))*(1/(exp(p1/T) + 1))*(1/(exp((sqrt(M_T**2 + pt_T**2)*cosh(y) - p1/T) +1))) ,0.333, 20./3,lambda thetap1: -pi+0.1, -0.1, lambda thetap1, p1: p1min(thetaP, p1,T,M_T),lambda thetap1,p1: p1max(thetaP,p1,T,M_T) )
return area[0]
#solving the integral
for ind in range(0, precision):
integral1d[ind]=formula151( M_T[ind])
print integral1d[ind]
plot(M_T,integral1d)
xlabel('M/T')
ylabel('prod rate')
title('thetaP =-3.12')
plt.yscale('log')
#plt.xscale('log')
show()
the error comes from line 57, where tplquad is used, and the full traceback is
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
/usr/lib/python2.7/dist-packages/IPython/utils/py3compat.pyc in execfile(fname, *where)
202 else:
203 filename = fname
--> 204 __builtin__.execfile(filename, *where)
/home/chiara/Scrivania/formula15a.py in <module>()
61 #solving the integral
62 for ind in range(0, precision):
---> 63 integral1d[ind]=formula151( M_T[ind])
64 print integral1d[ind]
65
/home/chiara/Scrivania/formula15a.py in formula151(M_T)
55
56 def formula151(M_T):
---> 57 area =tplquad(lambda p1, thetap1,T: pt_T*T*2*pi*5*(1/(18*pi**5))*sin(thetap1)*(p1/(sqrt(chi(thetap1,p1,thetaP,T,M_T))))*(1/(exp(p1/T) + 1))*(1/(exp((sqrt(M_T**2 + pt_T**2)*cosh(y) - p1/T) +1))) ,0.333, 20./3,lambda thetap1: -pi+0.1, -0.1, lambda thetap1, p1: p1min(thetaP, p1,T,M_T),lambda thetap1,p1: p1max(thetaP,p1,T,M_T) )
58 return area[0]
59
/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in tplquad(func, a, b, gfun, hfun, qfun, rfun, args, epsabs, epsrel)
498
499 """
--> 500 return dblquad(_infunc2,a,b,gfun,hfun,(func,qfun,rfun,args),epsabs=epsabs,epsrel=epsrel)
501
502
/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in dblquad(func, a, b, gfun, hfun, args, epsabs, epsrel)
433
434 """
--> 435 return quad(_infunc,a,b,(func,gfun,hfun,args),epsabs=epsabs,epsrel=epsrel)
436
437
/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
252 args = (args,)
253 if (weight is None):
--> 254 retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
255 else:
256 retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts)
/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
317 if points is None:
318 if infbounds == 0:
--> 319 return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
320 else:
321 return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _infunc(x, func, gfun, hfun, more_args)
379 def _infunc(x,func,gfun,hfun,more_args):
380 a = gfun(x)
--> 381 b = hfun(x)
382 myargs = (x,) + more_args
383 return quad(func,a,b,args=myargs)[0]
TypeError: 'float' object is not callable
NOTE: the functions "integral" and "integrand" are defined but finally not used...I just left them there
From the scipy.integrate docs, the signature of tplquad (which computes a numerical approximation to a triple integral) is:
scipy.integrate.tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)
where func is the function of three variables to be integrated, a and b are floating-point limits for the outer integral, gfun and hfun are functions of one variable giving the limits of the middle integral, and qfun and rfun are functions of two variables giving the limits of the innermost integral.
I had a hard time working out what was going on until I reformatted your code to be a bit more readable. Here's your call to tplquad, reformatted to make the line lengths a bit shorter:
area = tplquad(
lambda p1, thetap1, T: (
pt_T*T*2*pi*5*(1/(18*pi**5))*sin(thetap1)*
(p1/(sqrt(chi(thetap1,p1,thetaP,T,M_T))))*
(1/(exp(p1/T) + 1))*
(1/(exp((sqrt(M_T**2 + pt_T**2)*cosh(y) - p1/T) + 1)))
),
0.333, # a
20./3, # b
lambda thetap1: -pi + 0.1, # gfun
-0.1, # hfun
lambda thetap1, p1: p1min(thetaP, p1, T, M_T), # qfun
lambda thetap1, p1: p1max(thetaP, p1, T, M_T), # rfun
)
(As @user2357112 suggested, it would also aid readability to pull these lambda expressions out of the call and define them as separate functions. In particular, if you define a separate function for the integrand you'll be able to perform the computation piece by piece, and won't need to put one enormous expression onto one line.)
After the reformatting, it's much easier to see what the issue is: in your call to tplquad, you're passing the constant -0.1 for hfun. That won't work: in mathematics, one can (ab)use a constant value to represent a constant function, but programming languages (and some mathematicians, come to that) tend to be a bit more picky: you're going to need an actual function here. Replace -0.1 with lambda thetap1: -0.1.
By the way, I'm also a bit suspicious of your variable orders here. The doc page says that while the order of inputs for the integrand should be (z, y, x), gfun and hfun should be functions of just x, and qfun and rfun should be functions of (x, y) (in that order). That doesn't seem to match what you have.