def greenfct(zsrc, zrec, lsrc, lrec, depth, etaH, etaV, zetaH, zetaV, lambd):
    r"""Calculate Green's function for TM and TE.

    This is a modified version of empymod.kernel.greenfct(). See the original
    version for more information.

    """
    # GTM/GTE have shape (frequency, offset, lambda).
    # gamTM/gamTE have shape (frequency, offset, layer, lambda):

    for TM in [True, False]:

        # Define eta/zeta depending if TM or TE
        if TM:
            e_zH, e_zV, z_eH = etaH, etaV, zetaH   # TM: zetaV not used
        else:
            e_zH, e_zV, z_eH = zetaH, zetaV, etaH  # TE: etaV not used

        # Uppercase gamma
        Gam = np.sqrt((e_zH/e_zV)[:, None, :, None] *
                      (lambd*lambd)[None, :, None, :] +
                      (z_eH*e_zH)[:, None, :, None])

        # Gamma in receiver layer
        lrecGam = Gam[:, :, lrec, :]

        # Reflection (coming from below (Rp) and above (Rm) rec)
        Rp, Rm = reflections(depth, e_zH, Gam, lrec, lsrc, False)

        # Field propagators
        # (Up- (Wu) and downgoing (Wd), in rec layer); Eq 74
        if lrec != depth.size-1:  # No upgoing field prop. if rec in last
            ddepth = depth[lrec + 1] - zrec
            Wu = np.exp(-lrecGam*ddepth)
        else:
            Wu = np.full_like(lrecGam, 0+0j)
        if lrec != 0:     # No downgoing field propagator if rec in first
            ddepth = zrec - depth[lrec]
            Wd = np.exp(-lrecGam*ddepth)
        else:
            Wd = np.full_like(lrecGam, 0+0j)

        # Field at rec level (coming from below (Pu) and above (Pd) rec)
        Puu, Pud, Pdu, Pdd = fields(depth, Rp, Rm, Gam, lrec, lsrc, zsrc, TM)

        # Store in corresponding variable PT* = [T*uu, T*ud, T*du, T*dd]
        df = np.exp(-lrecGam*abs(zsrc - zrec))  # direct field
        fTM = Gam[:, :, lrec, :]/etaH[:, None, lrec, None]
        fTE = zetaH[:, None, lsrc, None]/Gam[:, :, lsrc, :]
        if TM:
            PTM = [Puu*Wu*fTM, Pud*Wu*fTM, Pdu*Wd*fTM, Pdd*Wd*fTM, -df*fTM]
        else:
            PTE = [Puu*Wu*fTE, Pud*Wu*fTE, Pdu*Wd*fTE, Pdd*Wd*fTE, df*fTE]

    # Return Green's functions
    return PTM, PTE