code: def f_BKLISCO(m1, m2):
      
      # q is defined to be in [0,1] for this formula
      q = numpy.minimum(m1/m2, m2/m1)
      return f_SchwarzISCO(m1+m2) * ( 1 + 2.8*q - 2.6*q*q + 0.8*q*q*q )
docs: Mass ratio dependent ISCO derived from estimates of the final spin
      of a merged black hole in a paper by Buonanno, Kidder, Lehner
      (arXiv:0709.3839).  See also arxiv:0801.4297v2 eq.(5)
  
      Parameters
      ----------
      m1 : float or numpy.array
          First component mass in solar mass units
      m2 : float or numpy.array
          Second component mass in solar mass units
  
      Returns
      -------
      f : float or numpy.array
          Frequency in Hz