code: def nltides_gw_phase_diff_isco(f_low, f0, amplitude, n, m1, m2):
      
      f0, amplitude, n, m1, m2, input_is_array = ensurearray(
          f0, amplitude, n, m1, m2)
  
      f_low = numpy.zeros(m1.shape) + f_low
  
      phi_l = nltides_gw_phase_difference(
                  f_low, f0, amplitude, n, m1, m2)
  
      f_isco = f_schwarzchild_isco(m1+m2)
  
      phi_i = nltides_gw_phase_difference(
                  f_isco, f0, amplitude, n, m1, m2)
  
      return formatreturn(phi_i - phi_l, input_is_array)
docs: Calculate the gravitational-wave phase shift bwtween
      f_low and f_isco due to non-linear tides.
  
      Parameters
      ----------
      f_low: float
          Frequency from which to compute phase. If the other
          arguments are passed as numpy arrays then the value
          of f_low is duplicated for all elements in the array
      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