code: def dipole_k(src, rec, depth, res, freq, wavenumber, ab=11, aniso=None,
               epermH=None, epermV=None, mpermH=None, mpermV=None, verb=2):
      r
      t0 = printstartfinish(verb)
      modl = check_model(depth, res, aniso, epermH, epermV, mpermH, mpermV,
                         False, verb)
      depth, res, aniso, epermH, epermV, mpermH, mpermV, _ = modl
      f = check_frequency(freq, res, aniso, epermH, epermV, mpermH, mpermV, verb)
      freq, etaH, etaV, zetaH, zetaV = f
      ab_calc, msrc, mrec = check_ab(ab, verb)
      src, nsrc = check_dipole(src, 'src', verb)
      rec, nrec = check_dipole(rec, 'rec', verb)
      off, angle = get_off_ang(src, rec, nsrc, nrec, verb)
      factAng = kernel.angle_factor(angle, ab, msrc, mrec)
      lsrc, zsrc = get_layer_nr(src, depth)
      lrec, zrec = get_layer_nr(rec, depth)
      if off.size == 1 and np.ndim(wavenumber) == 2:
          PJ0 = np.zeros((freq.size, wavenumber.shape[0], wavenumber.shape[1]),
                         dtype=complex)
          PJ1 = np.zeros((freq.size, wavenumber.shape[0], wavenumber.shape[1]),
                         dtype=complex)
      else:
          PJ0 = np.zeros((freq.size, off.size, wavenumber.size), dtype=complex)
          PJ1 = np.zeros((freq.size, off.size, wavenumber.size), dtype=complex)
      if ab_calc not in [36, ]:
          J0, J1, J0b = kernel.wavenumber(zsrc, zrec, lsrc, lrec, depth, etaH,
                                          etaV, zetaH, zetaV,
                                          np.atleast_2d(wavenumber), ab_calc,
                                          False, msrc, mrec, False)
          if J1 is not None:
              PJ1 += factAng[:, np.newaxis]*J1
              if ab in [11, 12, 21, 22, 14, 24, 15, 25]:
                  PJ1 /= off[:, None]
          if J0 is not None:
              PJ0 += J0
          if J0b is not None:
              PJ0 += factAng[:, np.newaxis]*J0b
      printstartfinish(verb, t0, 1)
      return np.squeeze(PJ0), np.squeeze(PJ1)