<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def hquad(zsrc, zrec, lsrc, lrec, off, factAng, depth, ab, etaH, etaV, zetaH, zetaV, xdirect, quadargs, use_ne_eval, msrc, mrec):
 r"""Hankel Transform using the ``QUADPACK`` library. This routine uses the ``scipy.integrate.quad`` module, which in turn makes use of the Fortran library ``QUADPACK`` (``qagse``). It is massively (orders of magnitudes) slower than either ``fht`` or ``hqwe``, and is mainly here for completeness and comparison purposes. It always uses interpolation in the wavenumber domain, hence it generally will not be as precise as the other methods. However, it might work in some areas where the others fail. The function is called from one of the modelling routines in :mod:`model`. Consult these modelling routines for a description of the input and output parameters. Returns ------- fEM : array Returns frequency-domain EM response. kcount : int Kernel count. For HQUAD, this is 1. conv : bool If true, QUAD converged. If not, <htarg> might have to be adjusted. """


    # Get quadargs
    rtol, atol, limit, a, b, pts_per_dec = quadargs

    # Get required lambdas
    la = np.log10(a)
    lb = np.log10(b)
    ilambd = np.logspace(la, lb, (lb-la)*pts_per_dec + 1)

    # Call the kernel
    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)

    # Interpolation in wavenumber domain: Has to be done separately on each PJ,
    # in order to work with multiple offsets which have different angles.
    # We check if the kernels are zero, to avoid unnecessary calculations.
    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

    # Pre-allocate output array
    fEM = np.zeros(off.size, dtype=complex)
    conv = True

    # Input-dictionary for quad
    iinp = {'a': a, 'b': b, 'epsabs': atol, 'epsrel': rtol, 'limit': limit}

    # Loop over offsets
    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 the electromagnetic field
    # Second argument (1) is the kernel count, last argument is only for QWE.
    return fEM, 1, conv