Численное двойное интегрирование, взвешенное функцией распределения вероятностей (PDF) в Python

У меня проблема с выполнением двойного определенного интеграла в функции, которая зависит от двух переменных (q,r) и имеет в ней один дополнительный интеграл. Функция, которую я хочу взвесить с помощью гауссовской функции:

F(q,r)=f(q,r)+int_{0,r}(h(q,r')dr')

И должен быть снова интегрирован, чтобы быть взвешенным с гауссовым:

I(q)=int_{0,inf}(F(q,r)^2*g(r)dr)

Гауссов g(r) центрируется по координате R.

Как видите, основная проблема заключается в том, что я смешиваю массивы со скалярами. Использование того же метода, который используется для гауссиана (np.ogrid и сумма по оси), может быть решением, но я не знаю, как его реализовать.

import numpy as np
from scipy.integrate import quad
import math as m

R=53.
R0=40.
delta=50.
c=2.
qm, rm = np.ogrid[0.0005:2.0:0.0005, 20:100:500j]

#normalized gauss function
#g(r)
def gauss_grid(r,Rmin,pd):
    def gauss(r,Rmin,pd):
        sigma=1.5
        return (1/sigma)*np.exp(-((r-Rmin)**2)/(2*sigma**2))
    gauss_grid = gauss(r,Rmin,pd)
    #normalization of gaussian
    gauss_grid /= np.sum(gauss_grid)
    return gauss_grid

#spherical function 
#f(q,r) 
def form(q,R):
    return (4/3)*m.pi*3*(np.sin(q*R)-q*R*np.cos(q*R))/(q**3)

#FINAL function
#I(q)
def helfand():
    def F(q,R):
        #integral (0,R) of h(q,r)
        def integral(q,Rmax):
            #h(q,r)
            def integrand(r,q):
                return np.sin(q*r)*(r**2)/(q*r*(1+np.exp(c*(R0-r))))
            return quad(integrand, 0, Rmax, args=(q))[0]
        return (form(q,R)+delta*integral(q,R))**2

    FF_hel=F(qm,rm)
    FF_hel *= gauss_grid(rm,R,pd)
    I=FF_hel.sum(axis=1)
    return I,qm.ravel()

helfand()

* ОБНОВЛЕНИЕ * * **

Я пробовал использовать библиотеку scipy.integrate (с quad), но не могу. Это похоже на то, что он не передает правильный аргумент (q) следующей функции. Вот очень упрощенная версия того, что я пытаюсь:

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt

R=53.
R0=41.
pd=15.
sigma=1.5

def I(q):
    #(function with integral inside) squared
    def FF(q,r):
        def integral_f(q,r):
            def f(r1,q):
                return np.sin(q*r1)
            return quad(f,0,r,args=(q))[0]

        def h(q,r):
            return (r*np.cos(q*r))
        return (h(q,r)+integral_f(q,r))**2

    #gaussian function normalized
    def g(r,R0):
        def gauss(r,R0):
            return (1/sigma)*np.exp(-((r-R0)**2)/(2*sigma**2))
        return gauss(r,R0)/(quad(gauss,0,np.inf,args=(R0))[0])

    #main function to be integrated with gaussian
    def function(r,q):
        return FF(q,r)*g(r,R)

    return quad(function,0,np.inf,args=(q))[0]

q=np.arange(0.001,1.,0.001)
plt.plot(q,I(q))

Ошибка говорит:

Поставляемая функция не возвращает допустимое число с плавающей запятой.