<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def wavenumber(zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH, zetaV, lambd, ab, xdirect, msrc, mrec, use_ne_eval):
 r"""Calculate wavenumber domain solution. Return the wavenumber domain solutions ``PJ0``, ``PJ1``, and ``PJ0b``, which have to be transformed with a Hankel transform to the frequency domain. ``PJ0``/``PJ0b`` and ``PJ1`` have to be transformed with Bessel functions of order 0 (:math:`J_0`) and 1 (:math:`J_1`), respectively. This function corresponds loosely to equations 105--107, 111--116, 119--121, and 123--128 in [HuTS15]_, and equally loosely to the file ``kxwmod.c``. [HuTS15]_ uses Bessel functions of orders 0, 1, and 2 (:math:`J_0, J_1, J_2`). The implementations of the *Fast Hankel Transform* and the *Quadrature-with-Extrapolation* in ``transform`` are set-up with Bessel functions of order 0 and 1 only. This is achieved by applying the recurrence formula .. math:: J_2(kr) = \frac{2}{kr} J_1(kr) - J_0(kr) \ . .. note:: ``PJ0`` and ``PJ0b`` could theoretically be added here into one, and then be transformed in one go. However, ``PJ0b`` has to be multiplied by ``factAng`` later. This has to be done after the Hankel transform for methods which make use of spline interpolation, in order to work for offsets that are not in line with each other. This function is called from one of the Hankel functions in :mod:`transform`. Consult the modelling routines in :mod:`model` for a description of the input and output parameters. If you are solely interested in the wavenumber-domain solution you can call this function directly. However, you have to make sure all input arguments are correct, as no checks are carried out here. """


    # ** CALCULATE GREEN'S FUNCTIONS
    # Shape of PTM, PTE: (nfreq, noffs, nfilt)
    PTM, PTE = greenfct(zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH,
                        zetaV, lambd, ab, xdirect, msrc, mrec, use_ne_eval)

    # ** AB-SPECIFIC COLLECTION OF PJ0, PJ1, AND PJ0b

    # Pre-allocate output
    PJ0 = None
    PJ1 = None
    PJ0b = None

    # Calculate Ptot which is used in all cases
    Ptot = (PTM + PTE)/(4*np.pi)

    # If rec is magnetic switch sign (reciprocity MM/ME => EE/EM).
    if mrec:
        sign = -1
    else:
        sign = 1

    # Group into PJ0 and PJ1 for J0/J1 Hankel Transform
    if ab in [11, 12, 21, 22, 14, 24, 15, 25]:    # Eqs 105, 106, 111, 112,
        # J2(kr) = 2/(kr)*J1(kr) - J0(kr)         #     119, 120, 123, 124
        if ab in [14, 22]:
            sign *= -1
        PJ0b = sign/2*Ptot*lambd
        PJ1 = -sign*Ptot
        if ab in [11, 22, 24, 15]:
            if ab in [22, 24]:
                sign *= -1
            PJ0 = sign*(PTM - PTE)/(8*np.pi)*lambd

    elif ab in [13, 23, 31, 32, 34, 35, 16, 26]:  # Eqs 107, 113, 114, 115,
        PJ1 = sign*Ptot*lambd*lambd               # .   121, 125, 126, 127
        if ab in [34, 26]:
            PJ1 *= -1

    elif ab in [33, ]:                            # Eq 116
        PJ0 = sign*Ptot*lambd*lambd*lambd

    # Return PJ0, PJ1, PJ0b
    return PJ0, PJ1, PJ0b