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. / 12. - 19. / 3. * eta) * chiS
    beta += 113. / 12. * delta * chiA
    sigma = eta / 48. * (474 * spinspin)
    sigma += (1 - 2 * eta) * (81. / 16. * (chiS * chiS + chiA * chiA))
    sigma += delta * (81. / 8. * (chiS * chiA))
    gamma = (732985. / 2268. - 24260. / 81. * eta - \
            340. / 9. * eta * eta) * chiS
    gamma += (732985. / 2268. + 140. / 9. * eta) * delta * chiA
    return beta, sigma, gamma
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.