def merger_final_spin(q, a_1, a_2):
    """Compute final BH spin after a GW merger between two BHs. Expression from
    Tichy & Marronetti 2008.

    Parameters
    ----------
    q : `float/array`
        Mass ratio of binary
    a_1 : `float/array`
        Spin of primary BH
    a_2 : `float/array`
        Spin of secondary BH

    Returns
    -------
    a_f : `float/array`
        Final spin of the BH after merger
    """
    v = q / (1 + q)**2
    a_f = 0.686 * (5.04 * v - 4.16 * v**2)\
        + 0.4 * ((a_1 / (0.632 + 1 / q)**2) + (a_2 / (0.632 + q)**2))
    return a_f