def empy_hankel(ftype, zsrc, zrec, res, freqtime, depth=None, aniso=None,
                epermH=None, epermV=None, mpermH=None, mpermV=None,
                htarg=None, verblhs=0, verbrhs=0):
    r"""Numerical transform pair with empymod.

    All parameters except ``ftype``, ``verblhs``, and ``verbrhs`` correspond to
    the input parameters to ``empymod.dipole``. See there for more information.

    Note that if depth=None or [], the analytical full-space solutions will be
    used (much faster).

    Parameters
    ----------
    ftype : str or list of strings
        Either of: {'j0', 'j1', 'j2', ['j0', 'j1']}

        - 'j0': Analyze J0-term with ab=11, angle=45°
        - 'j1': Analyze J1-term with ab=31, angle=0°
        - 'j2': Analyze J0- and J1-terms jointly with ab=12, angle=45°
        - ['j0', 'j1']: Same as calling empy_hankel twice, once with 'j0' and
                        one with 'j1'; can be provided like this to
                        fdesign.design.

    verblhs, verbrhs: int
        verb-values provided to empymod for lhs and rhs.

    Note that ftype='j2' only works for fC, not for fI.

    """

    # Loop over ftypes, if there are several
    if isinstance(ftype, list):
        out = []
        for f in ftype:
            out.append(empy_hankel(f, zsrc, zrec, res, freqtime, depth, aniso,
                                   epermH, epermV, mpermH, mpermV, htarg,
                                   verblhs, verbrhs))
        return out

    # Collect model
    model = {'src': [0, 0, zsrc],
             'depth': depth,
             'res': res,
             'aniso': aniso,
             'epermH': epermH,
             'epermV': epermV,
             'mpermH': mpermH,
             'mpermV': mpermV}

    # Finalize model depending on ftype
    if ftype == 'j0':  # J0: 11, 45°
        model['ab'] = 11
        x = 1/np.sqrt(2)
        y = 1/np.sqrt(2)

    elif ftype == 'j1':  # J1: 31, 0°
        model['ab'] = 31
        x = 1
        y = 0

    elif ftype == 'j2':  # J2: 12, 45°
        model['ab'] = 12
        x = 1/np.sqrt(2)
        y = 1/np.sqrt(2)

    # rhs: empymod.model.dipole
    # If depth=[], the analytical full-space solution will be used internally
    def rhs(r):
        out = dipole(rec=[r*x, r*y, zrec], ht='qwe', xdirect=True,
                     verb=verbrhs, htarg=htarg, freqtime=freqtime, **model)
        return out

    # lhs: empymod.model.dipole_k
    def lhs(k):
        lhs0, lhs1 = dipole_k(rec=[x, y, zrec], wavenumber=k, verb=verblhs,
                              freq=freqtime, **model)
        if ftype == 'j0':
            return lhs0
        elif ftype == 'j1':
            return lhs1
        elif ftype == 'j2':
            return (lhs0, lhs1)

    return Ghosh(ftype, lhs, rhs)
r"""Numerical transform pair with empymod.

    All parameters except ``ftype``, ``verblhs``, and ``verbrhs`` correspond to
    the input parameters to ``empymod.dipole``. See there for more information.

    Note that if depth=None or [], the analytical full-space solutions will be
    used (much faster).

    Parameters
    ----------
    ftype : str or list of strings
        Either of: {'j0', 'j1', 'j2', ['j0', 'j1']}

        - 'j0': Analyze J0-term with ab=11, angle=45°
        - 'j1': Analyze J1-term with ab=31, angle=0°
        - 'j2': Analyze J0- and J1-terms jointly with ab=12, angle=45°
        - ['j0', 'j1']: Same as calling empy_hankel twice, once with 'j0' and
                        one with 'j1'; can be provided like this to
                        fdesign.design.

    verblhs, verbrhs: int
        verb-values provided to empymod for lhs and rhs.

    Note that ftype='j2' only works for fC, not for fI.