Using nquad for a double integral

Question

Having a problem here. Here's my code so far:

from scipy import integrate
import math
import numpy as np

a = 0.250
s02 = 214.0
a_s = 0.0163

def integrand(r, R, s02, a_s, a):
        return 2.0 * r * (r/a)**(-0.1) * (1.0 + (r**2/a**2))**(-2.45)\\
*(math.sqrt(r**2 - R**2))**(-1.0) * (a_s/(1 + (R-0.0283)**2/a_s**2 ))

def bounds_R(s02, a_s, a):
        return [0, np.inf]
def bounds_r(R, s02, a_s, a):
        return [R, np.inf]

result = integrate.nquad(integrand, [bounds_r(R, s02, a_s, a), bounds_R(s02, a_s, a)])

a, s02 and a_s are constants. I need to perform the first integral over r, and then the second integral over R. The problem I think is the fact that R appears in the limits for the first integral (something called an Abel transform). Tried a few things and every time getting an error that there are either too few arguments in the boundary functions or too little.

Please help!

Answer 1

If you write integrate.nquad(integrand, [bounds_r(R, s02, a_s, a), bounds_R(s02, a_s, a)]), python is expecting you to affect a value to R. But you didn't because integration is carried out over R.

This syntax should work :

result = integrate.nquad(integrand, [bounds_r, bounds_R], args=(s02,a_s,a))

Take a look to the second example in the documentation of integrate.nquad.