code: def get_beta_sigma_from_aligned_spins(eta, spin1z, spin2z):
      
      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
docs: 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.