def nltides_gw_phase_difference(f, f0, amplitude, n, m1, m2):
    """Calculate the gravitational-wave phase shift bwtween
    f and f_coalescence = infinity due to non-linear tides.
    To compute the phase shift between e.g. f_low and f_isco,
    call this function twice and compute the difference.

    Parameters
    ----------
    f: float or numpy.array
        Frequency from which to compute phase
    f0: float or numpy.array
        Frequency that NL effects switch on
    amplitude: float or numpy.array
        Amplitude of effect
    n: float or numpy.array
        Growth dependence of effect
    m1: float or numpy.array
        Mass of component 1
    m2: float or numpy.array
        Mass of component 2

    Returns
    -------
    delta_phi: float or numpy.array
        Phase in radians
    """
    f, f0, amplitude, n, m1, m2, input_is_array = ensurearray(
        f, f0, amplitude, n, m1, m2)

    delta_phi = numpy.zeros(m1.shape)

    f_ref, _, phi_of_f_factor = nltides_coefs(amplitude, n, m1, m2)

    mask = f <= f0
    delta_phi[mask] = - phi_of_f_factor[mask] * (f0[mask]/f_ref)**(n[mask]-3.)

    mask = f > f0
    delta_phi[mask] = - phi_of_f_factor[mask] * (f[mask]/f_ref)**(n[mask]-3.)

    return formatreturn(delta_phi, input_is_array)