def reflections(depth, e_zH, Gam, lrec, lsrc, use_ne_eval):
    r
    for plus in [True, False]:
        if plus:
            pm = 1
            layer_count = np.arange(depth.size-2, min(lrec, lsrc)-1, -1)
            izout = abs(lsrc-lrec)
            minmax = max(lrec, lsrc)
        else:
            pm = -1
            layer_count = np.arange(1, max(lrec, lsrc)+1, 1)
            izout = 0
            minmax = -min(lrec, lsrc)
        shiftplus = lrec < lsrc and lrec == 0 and not plus
        shiftminus = lrec > lsrc and lrec == depth.size-1 and plus
        if shiftplus or shiftminus:
            izout -= pm
        Ref = np.zeros((Gam.shape[0], Gam.shape[1], abs(lsrc-lrec)+1,
                        Gam.shape[3]), dtype=complex)
        for iz in layer_count:
            e_zHa = e_zH[:, None, iz+pm, None]
            Gama = Gam[:, :, iz, :]
            e_zHb = e_zH[:, None, iz, None]
            Gamb = Gam[:, :, iz+pm, :]
            if use_ne_eval:
                rlocstr = "(e_zHa*Gama - e_zHb*Gamb)/(e_zHa*Gama + e_zHb*Gamb)"
                rloc = use_ne_eval(rlocstr)
            else:
                rloca = e_zHa*Gama
                rlocb = e_zHb*Gamb
                rloc = (rloca - rlocb)/(rloca + rlocb)
            if iz == layer_count[0]:
                tRef = rloc.copy()
            else:
                ddepth = depth[iz+1+pm]-depth[iz+pm]
                if use_ne_eval:
                    term = use_ne_eval("tRef*exp(-2*Gamb*ddepth)")
                    tRef = use_ne_eval("(rloc + term)/(1 + rloc*term)")
                else:
                    term = tRef*np.exp(-2*Gamb*ddepth)
                    tRef = (rloc + term)/(1 + rloc*term)
            if lrec != lsrc and pm*iz <= minmax:
                Ref[:, :, izout, :] = tRef[:]
                izout -= pm
        if lsrc == lrec and layer_count.size > 0:
            Ref = tRef
        if plus:
            Rm = Ref
        else:
            Rp = Ref
    return Rm, Rp