Numerical integration of highly oscillating 1-D integrand (containing Bessel functions) in python

I am attempting to numerically evaluate a real-valued integrand composed of several Bessel functions (of the first and second kind). The integrand is oscillating and decaying, and needs to be evaluated between 0 and +∞. So far, my attempts using the scipy.integrate subpackage (quad and fixed_quad) have been unsuccessful. The evaluated value jumps around when in reality it should be smooth. For certain sets of parameter values, I also receive warnings: "IntegrationWarning: The integral is probably divergent, or slowly convergent." (it is known to be convergent) or "IntegrationWarning: The maximum number of subdivisions (50) has been achieved."

The equation is from: http://dx.doi.org/10.1029/WR003i001p00241

It is also available here: http://www.aqtesolv.com/papadopu.htm

Thank you for any assistance you have in numerical integration of fussy functions in python...

CODE EXAMPLE

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import special as sps
import scipy.integrate as integrate

# define constants and variables (SI mks units):
r_w = 0.15
r_c = 0.16
b = 10
S_s = 1E-6
Q = 0.001
S = S_s*b
K=1E-8
T=K*b
alpha = (r_w**2)*S/r_c**2
def r_D(r):
    return r/r_w
def u(r,t):
    return r**2*S/(4*T*t)
def B(beta):
    return beta*sps.jv(0,beta) - 2*alpha*sps.jv(1,beta)
def A(beta):
    return beta*sps.yn(0,beta) - 2*alpha*sps.yn(1,beta)
def intd_f(beta,r,t): 
    TOP = (1-np.exp(-beta**2*r_D(r)**2)/(4*u(r,t)))*( sps.jv(0,beta*r_D(r))*A(beta) - sps.yn(0,beta*r_D(r))*B(beta) )
    BOT = (A(beta)**2 + B(beta)**2)*beta**2
    return TOP/BOT
def s(r,t):
    banana = (2*alpha*Q)/(np.pi**2*K*b)
    apple = integrate.quad(intd_f, 0, np.inf, args=(r,t))[0]
    #apple = integrate.fixed_quad(intd_f, 0, 1E100, args=(r,t))[0] # another option I have tried
    return banana*apple

#%% simple example usage...
r=np.arange(1,10,.1)
t=60*60*24*pd.Series([1/24,1,365,3650])
plt.figure()
for tt in t:
    print('time=',tt)
    snow=[]
    for rr in r:
        print('r=',rr)
        snow.append(s(rr,tt))
    plt.subplot(2,1,1)
    plt.plot(r,snow,label='t='+str(tt/(60*60*24))+'d')
    plt.subplot(2,1,2)
    plt.semilogy(r,np.abs(snow))    
plt.subplot(2,1,1)
plt.legend()

I have bad news for you: Oscillatory integrands must be dealt with by special methods, for example Filon's method. Reading through the scipy docs, I see no way to do this. However, it appears MPMath implements precisely what you are looking for, and even gives an example using Bessel functions.