code: 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