Solving a double integral with scipy dblquad

Question

I am trying to solve the double integral of this complicated function. The function contains 3 symbolic variables which are defined with sympy.symbols. My goal is to integrate the function with respect to only two of the variables. sym.integrate run for 2 hours with no results to show. I tried numerical integration with scipy.integrate.dblquad. But I am having troubles which I suspect is due to the third symbolic variable.

Is there a way to do this?

Problem summarized:

sym.symbols('x y z')
My_function(x,y,z)

Integrate My_function with respect to x and y (Both from 0 to inf i.e. definite integral).

h, t, w, r, q = sym.symbols('h t  w r q')  # Define symbols

wn = 2.0
alpha = 1.76
beta = 1.59
a0 = 0.7
a1 = 0.282
a2 = 0.167
b0 = 0.07
b1 = 0.3449
b2 = -0.2073
psi = 0.05

F_H = 1.0 - sym.exp(-(h / alpha) ** beta)
mu_h = a0 + a1 * h ** a2
sig_h = b0 + b1 * sym.exp(b2 * h)
F_TIH = (1 / 2) * (1 + sym.erf((sym.log(t) - mu_h) / (sig_h * sym.sqrt(2))))

f_h = sym.diff(F_H, h)
f_tzIhs = sym.diff(F_TIH, t)
f_S = f_h * f_tzIhs

H = (1.0 - (w / wn) ** 2.0 + 1.0j * 2.0 * psi * w / wn) ** (-1.0)
S_eta_h_t = h ** 2.0 * t / (8.0 * pi ** 2.0) * (w * t / (2.0 * pi)) ** (-5.0) * sym.exp(-1.0 / pi * (w * t / (2.0 * pi)) ** (-4.0))

S_RIS_hu_tu = abs(H) ** 2.0 * S_eta_h_t

m0_s = sym.integrate((w ** 0 * S_RIS_hu_tu), (w, 0, np.inf))
m0_s.doit()

m2_s = sym.integrate((w ** 2 * S_RIS_hu_tu), (w, 0, np.inf))
m2_s.doit()

v_rIs = 1 / (2 * pi) * sym.sqrt(m2_s / m0_s) * sym.exp(-r ** 2 / (2 * m0_s))

fun = v_rIs * f_S

# The integral I am trying to solve is a function of h,t and r. 
integ_ht = sym.integrate(fun,(h,0,np.inf),(t,0,np.inf))