docs: def kbn_conf (N, B, CL):
      """Given a (integer) number of observed Poisson events `N` and a (real)
      expected number of background events `B` and a confidence limit `CL`
      (between 0 and 1), return the confidence interval on the source event
      rate.
  
      Returns: (Smin, Smax)
  
      This interval is calculated using the Bayesian formalism of Kraft, Burrows, &
      Nousek (1991), which assumes no uncertainty in `B` and returns the smallest
      possible interval that satisfies the above properties.
  
      Example: in a certain time interval, 3 events were recorded. Based on
      external knowledge, it is expected that on average 0.5 background events will
      be recorded in the same interval. The 95% confidence interval on the source
      event rate is
  
      >>> kbn_conf.kbn_conf (3, 0.5, 0.95)
      <<< (0.22156, 7.40188)
  
      which agrees with the entry in Table 2 of KBN91.
  
      Reference info: 1991ApJ...374..344K, doi:10.1086/170124
  
      """
  
      tol = 1e-6
  
      origN = N
      try:
          N = int (N)
          assert N == origN
      except Exception:
          raise ValueError ('N must be an integer')
  
      CL = float (CL)
      if CL <= 0. or CL >= 1.:
          raise ValueError ('CL must be between 0 and 1, noninclusive')
  
      B = float (B)
      if B < 0:
          raise ValueError ('B must be nonnegative')
  
      # OK, arg-checking is out of the way. Precompute some things ...
  
      C = _cconst (N, B)
      f = lambda s: _fcnb (C, N, B, s)
  
      # The goal is find Smin and Smax such that the integral of _fnb between
      # Smin and Smax is CL, and _fnb (Smin) = _fnb (Smax). Follow the
      # suggestion in Kraft, Burrows, & Nousek (1991) and start at the
      # maximum-probability value, integrating outwards trying to maintain the
      # constraints. We have to be careful because smin cannot go below zero,
      # and to ignore the enormous typo ("local maximum at S = B + N"!) in the
      # paper!
  
      smin = smax = max (N - B, 0.)
      fmin = f (smin)
      fmax = f (smax)
      conf = 0.
  
      while conf < CL:
          if smin == 0. or fmin < fmax:
              stepsize = max (0.2 * abs (CL - conf) / CL / fmax, tol)
              conf += quad (f, smax, smax + stepsize)[0]
              smax += stepsize
              fmax = f (smax)
          else:
              stepsize = max (min (0.2 * abs (CL - conf) / CL / fmin, 0.1 * smin), tol)
              if smin - stepsize < tol:
                  conf += quad (f, 0, smin)[0]
                  smin = 0.
              else:
                  conf += quad (f, smin - stepsize, smin)[0]
                  smin -= stepsize
              fmin = f (smin)
  
      return smin, smax