code: def nltides_gw_phase_difference(f, f0, amplitude, n, m1, m2):
      
      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)
docs: 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