# Numerical Double Integration using numba and scipy

import numpy as np
from scipy import integrate
from scipy.special import j1
from scipy.special import j0
from numpy import sin
from numpy import cos
from numpy import pi
import time
import numba as nb
from numba import jit
from numba import cfunc
from numba.types import intc, CPointer, float64
from scipy import LowLevelCallable
from matplotlib import pyplot as plt

q = np.linspace(0.03, 1.0, 1000)
start = time.time()
length = 20000
pd = 0.1
start = time.time()

def jit_integrand_function(integrand_function):
jitted_function = jit(integrand_function, nopython=True)
@cfunc(float64(intc, CPointer(float64)),error_model="numpy",fastmath=True)
def wrapped(n, xx):
ar = nb.carray(xx, n)
return jitted_function(ar[0], ar[1], ar[2])
return LowLevelCallable(wrapped.ctypes)

@jit_integrand_function
def f(theta,z,q):
def ff(theta, z, q):
l = length / 2
qz = q * sin(theta)
qr = q * cos(theta)
return (4 * pi * sin(qz * l) * z * j1(qr * z) / (qr * qz)) ** 2

s = pd * c
return (1 / s)) * np.exp(-0.5 * (((z - radius) / s) ** 2))

def lower_inner(z):
return 0.

def upper_inner(z):
return 120.

y = np.empty(len(q))
for n in range(len(q)):
y[n] = integrate.dblquad(f, 0, 1*pi/180, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()
print(end - start)
plt.loglog(q,y)
plt.show()


Also, I was wondering if I can use GPU for the above code? How much gain I would have for the trouble(I have zero knowledge of programming for GPU use)?