def get_beta_sigma_from_aligned_spins(eta, spin1z, spin2z):
    """
    Calculate the various PN spin combinations from the masses and spins.
    See <http://arxiv.org/pdf/0810.5336v3.pdf>.

    Parameters
    -----------
    eta : float or numpy.array
        Symmetric mass ratio of the input system(s)
    spin1z : float or numpy.array
        Spin(s) parallel to the orbit of the heaviest body(ies)
    spin2z : float or numpy.array
        Spin(s) parallel to the orbit of the smallest body(ies)

    Returns
    --------
    beta : float or numpy.array
        The 1.5PN spin combination
    sigma : float or numpy.array
        The 2PN spin combination
    gamma : float or numpy.array
        The 2.5PN spin combination
    chis : float or numpy.array
        (spin1z + spin2z) / 2.
    """
    chiS = 0.5 * (spin1z + spin2z)
    chiA = 0.5 * (spin1z - spin2z)
    delta = (1 - 4 * eta) ** 0.5
    spinspin = spin1z * spin2z
    beta = (113.0 / 12.0 - 19.0 / 3.0 * eta) * chiS
    beta += 113.0 / 12.0 * delta * chiA
    sigma = eta / 48.0 * (474 * spinspin)
    sigma += (1 - 2 * eta) * (81.0 / 16.0 * (chiS * chiS + chiA * chiA))
    sigma += delta * (81.0 / 8.0 * (chiS * chiA))
    gamma = (732985.0 / 2268.0 - 24260.0 / 81.0 * eta - 340.0 / 9.0 * eta * eta
        ) * chiS
    gamma += (732985.0 / 2268.0 + 140.0 / 9.0 * eta) * delta * chiA
    return beta, sigma, gamma
