def hquad(zsrc, zrec, lsrc, lrec, off, factAng, depth, ab, etaH, etaV, zetaH,
          zetaV, xdirect, quadargs, use_ne_eval, msrc, mrec):
    r
    rtol, atol, limit, a, b, pts_per_dec = quadargs
    la = np.log10(a)
    lb = np.log10(b)
    ilambd = np.logspace(la, lb, (lb-la)*pts_per_dec + 1)
    PJ0, PJ1, PJ0b = kernel.wavenumber(zsrc, zrec, lsrc, lrec, depth, etaH,
                                       etaV, zetaH, zetaV,
                                       np.atleast_2d(ilambd), ab, xdirect,
                                       msrc, mrec, use_ne_eval)
    if PJ0 is not None:
        sPJ0r = iuSpline(np.log(ilambd), PJ0.real)
        sPJ0i = iuSpline(np.log(ilambd), PJ0.imag)
    else:
        sPJ0r = None
        sPJ0i = None
    if PJ1 is not None:
        sPJ1r = iuSpline(np.log(ilambd), PJ1.real)
        sPJ1i = iuSpline(np.log(ilambd), PJ1.imag)
    else:
        sPJ1r = None
        sPJ1i = None
    if PJ0b is not None:
        sPJ0br = iuSpline(np.log(ilambd), PJ0b.real)
        sPJ0bi = iuSpline(np.log(ilambd), PJ0b.imag)
    else:
        sPJ0br = None
        sPJ0bi = None
    fEM = np.zeros(off.size, dtype=complex)
    conv = True
    iinp = {'a': a, 'b': b, 'epsabs': atol, 'epsrel': rtol, 'limit': limit}
    for i in range(off.size):
        fEM[i], tc = quad(sPJ0r, sPJ0i, sPJ1r, sPJ1i, sPJ0br, sPJ0bi, ab,
                          off[i], factAng[i], iinp)
        conv *= tc
    return fEM, 1, conv