code: def xi_eq(x, kappa, chi_eff, q):
      
      return x**3*(x**2-3*kappa*x+2*chi_eff*kappa*math.sqrt(kappa*x))-3*q*(x**2-2*kappa*x+(chi_eff*kappa)**2)
docs: The roots of this equation determine the orbital radius
      at the onset of NS tidal disruption in a nonprecessing
      NS-BH binary [(7) in Foucart PRD 86, 124007 (2012)]
  
      Parameters
      -----------
      x: float
          orbital separation in units of the NS radius
      kappa: float
          the BH mass divided by the NS radius
      chi_eff: float
          the BH dimensionless spin parameter
      q: float
          the binary mass ratio (BH mass / NS mass)
  
      Returns
      ----------
      float
          x**3*(x**2-3*kappa*x+2*chi_eff*kappa*sqrt[kappa*x)
          -3*q*(x**2-2*kappa*x+(chi_eff*kappa)**2)