def Gp_Gm_vectors(ssalb, kestar, g, mu):
    """return Gp and Gm vectors at one wavelength

    :param ssalb: single scattering albedo of each layer (no unit)
    :type ssalb: array
    :param kestar: delta-eddington AFEC (no unit)
    :type kestar: array
    :param g: asymmetry factor  of each layer (no unit)
    :type g: array
    :param mu: cosine of the incident angle

    :returns: Gp and Gm
"""
    g_star = g / (1 + g)
    ssalb_star = ssalb * (1 - g ** 2) / (1 - g ** 2 * ssalb)
    gamma1 = 0.25 * (7 - ssalb_star * (4 + 3 * g_star))
    gamma2 = -0.25 * (1 - ssalb_star * (4 - 3 * g_star))
    gamma3 = 0.25 * (2 - 3 * g_star * mu)
    gamma4 = 0.25 * (2 + 3 * g_star * mu)
    G = mu ** 2 * ssalb_star / ((kestar * mu) ** 2 - 1)
    Gp = G * ((gamma1 - 1 / mu) * gamma3 + gamma2 * gamma4)
    Gm = G * ((gamma1 + 1 / mu) * gamma4 + gamma2 * gamma3)
    return Gp, Gm
