code: def bh_effective_spin(chi,incl):
      
      if incl == 0:
          chi_eff = chi
      else:
          # ISSO radius for the given BH spin magnitude and inclination
          rISSO = PG_ISSO_solver(chi,incl)
          # Angle at which the branch of positive solutions has its minumum
          incl_flip = scipy.optimize.fmin(pos_branch, math.pi/4, args=tuple([chi]), full_output=False, disp=False)[-1]
          # Use incl_flip to determine the initial guess: the sign difference
          # in the initial_guess ensures that chi_eff has the correct sign
          if incl>incl_flip:
              initial_guess = -1.1
          else:
              initial_guess = 1.0
          chi_eff = scipy.optimize.fsolve(ISCO_eq_chi_first, initial_guess, args=(rISSO))
  
      return chi_eff
docs: Determines the effective [as defined in Stone, Loeb,
      Berger, PRD 87, 084053 (2013)] aligned dimensionless
      spin parameter of a NS-BH binary with tilted BH spin.
      This means finding the root chi_eff of
      ISCO_eq_chi_first(chi_eff, PG_ISSO_solver(chi,incl))
      with the correct sign.
  
      Parameters
      -----------
      chi: float
          the BH dimensionless spin parameter
      incl: float
          the inclination angle between the BH spin and the orbital
          angular momentum in radians
  
      Returns
      ----------
      chi_eff: float
          the effective dimensionless spin parameter solution