docs: def Ntubes_Phadkeb(DBundle, Do, pitch, Ntp, angle=30):
      r'''Using tabulated values and correction factors for number of passes,
      the highly accurate method of [1]_ is used to obtain the tube count
      of a given tube bundle outer diameter for a given tube size and pitch.
  
      Parameters
      ----------
      DBundle : float
          Outer diameter of tube bundle, [m]
      Do : float
          Tube outer diameter, [m]
      pitch : float
          Pitch; distance between two orthogonal tube centers, [m]
      Ntp : int
          Number of tube passes, [-]
      angle : float, optional
          The angle the tubes are positioned; 30, 45, 60 or 90, [degrees]
  
      Returns
      -------
      Nt : int
          Total number of tubes that fit in the heat exchanger, [-]
  
      Notes
      -----
      For single-pass cases, the result is exact, and no tubes need to be removed
      for any reason. For 4, 6, 8 pass arrangements, a number of tubes must be 
      removed to accomodate pass partition plates. The following assumptions
      are involved with that:
          * The pass partition plate is where a row of tubes would have been. 
            Only one or two rows are assumed affected.
          * The thickness of partition plate is < 70% of the tube outer diameter.
          * The distance between the centerline of the partition plate and the 
            centerline of the nearest row of tubes is equal to the pitch.    
      
      This function will fail when there are more than 100,000 tubes.
      [1]_ tabulated values up to approximately 3,000 tubes derived with 
      number theory. The sequesnces of integers were identified in the
      On-Line Encyclopedia of Integer Sequences (OEIS), and formulas listed in
      it were used to generate more coefficient to allow up to 100,000 tubes.
      The integer sequences are A003136, A038590, A001481, and A057961. The 
      generation of coefficients for A038590 is very slow, but the rest are
      reasonably fast.
      
      The number of tubes that fit generally does not increase one-by-one, but by
      several.
      
      >>> Ntubes_Phadkeb(DBundle=1.007, Do=.028, pitch=.036, Ntp=2, angle=45.)
      558
      >>> Ntubes_Phadkeb(DBundle=1.008, Do=.028, pitch=.036, Ntp=2, angle=45.)
      574
      
      Because a pass partition needs to be installed in multiple tube pass
      shells, more tubes fit in an exchanger the fewer passes are used.
      
      >>> Ntubes_Phadkeb(DBundle=1.008, Do=.028, pitch=.036, Ntp=1, angle=45.)
      593
  
      Examples
      --------
      >>> Ntubes_Phadkeb(DBundle=1.200-.008*2, Do=.028, pitch=.036, Ntp=2, angle=45.)
      782
  
      References
      ----------
      .. [1] Phadke, P. S., Determining tube counts for shell and tube
         exchangers, Chem. Eng., September, 91, 65-68 (1984).
      '''
      if DBundle <= Do*Ntp:
          return 0
      
      if Ntp == 6:
          e = 0.265
      elif Ntp == 8:
          e = 0.404
      else:
          e = 0.
  
      r = 0.5*(DBundle - Do)/pitch
      s = r*r
      Ns, Nr = floor(s), floor(r)
      # If Ns is between two numbers, take the smaller one
      # C1 is the number of tubes for a single pass arrangement.
      if angle == 30 or angle == 60:
          i = np.searchsorted(triangular_Ns, Ns, side='right')
          C1 = int(triangular_C1s[i-1])
      elif angle == 45 or angle == 90:
          i = np.searchsorted(square_Ns, Ns, side='right')
          C1 = int(square_C1s[i-1])
  
      Cx = 2*Nr + 1.
  
      # triangular and rotated triangular
      if (angle == 30 or angle == 60):
          w = 2*r/3**0.5
          Nw = floor(w)
          if Nw % 2 == 0:
              Cy = 3*Nw
          else:
              Cy = 3*Nw + 1
          if Ntp == 2:
              if angle == 30 :
                  C2 = C1 - Cx
              else:
                  C2 = C1 - Cy - 1
          else: # 4 passes, or 8; this value is needed
              C4 = C1 - Cx - Cy
  
      if (angle == 30 or angle == 60) and (Ntp == 6 or Ntp == 8):
          if angle == 30: # triangular
              v = 2*e*r/3**0.5 + 0.5
              Nv = floor(v)
              u = 3**0.5*Nv/2.
              if Nv % 2 == 0:
                  z = (s-u*u)**0.5
              else:
                  z = (s-u*u)**0.5 - 0.5
              Nz = floor(z)
              if Ntp == 6:
                  C6 = C1 - Cy - 4*Nz - 1
              else:
                  C8 = C4 - 4*Nz
          else: # rotated triangular
              v = 2.*e*r
              Nv = floor(v)
              u1 = 0.5*Nv
              z = (s - u1*u1)**0.5
              w1 = 2*z/2**0.5
  #            w1 = 2**2**0.5 # WRONG
              u2 = 0.5*(Nv + 1)
              zs = (s-u2*u2)**0.5
              w2 = 2.*zs/3**0.5
              if Nv%2 == 0:
                  z1 = 0.5*w1
                  z2 = 0.5*(w2+1)
              else:
                  z1 = 0.5*(w1+1)
                  z2 = 0.5*w2
              Nz1 = floor(z1)
              Nz2 = floor(z2)
              if Ntp == 6:
                  C6 = C1 - Cx - 4.*(Nz1 + Nz2)
              else: # 8
                  C8 = C4 - 4.*(Nz1 + Nz2)
  
      if (angle == 45 or angle == 90):
          if angle == 90:
              Cy = Cx - 1.
              # eq 6 or 8 for c2 or c4
              if Ntp == 2:
                  C2 = C1 - Cx
              else: # 4 passes, or 8; this value is needed
                  C4 = C1 - Cx - Cy
          else: # rotated square
              w = r/2**0.5
              Nw = floor(w)
              Cx = 2.*Nw + 1
              Cy = Cx - 1
              if Ntp == 2:
                  C2 = C1 - Cx
              else: # 4 passes, or 8; this value is needed
                  C4 = C1 - Cx - Cy
  
      if (angle == 45 or angle == 90) and (Ntp == 6 or Ntp == 8):
          if angle == 90:
              v = e*r + 0.5
              Nv = floor(v)
              z = (s - Nv*Nv)**0.5
              Nz = floor(z)
              if Ntp == 6:
                  C6 = C1 - Cy - 4*Nz - 1
              else:
                  C8 = C4 - 4*Nz
          else:
              w = r/2**0.5
              Nw = floor(w)
              Cx = 2*Nw + 1
  
              v = 2**0.5*e*r
              Nv = floor(v)
              u1 = Nv/2**0.5
              z = (s-u1*u1)**0.5
              w1 = 2**0.5*z
              u2 = (Nv + 1)/2**0.5
              zs = (s-u2*u2)**0.5
              w2 = 2**0.5*zs
              # if Nv is odd, 21a and 22a. If even, 21b and 22b. Nz1, Nz2
              if Nv %2 == 0:
                  z1 = 0.5*w1
                  z2 = 0.5*(w2 + 1)
              else:
                  z1 = 0.5*(w1 + 1)
                  z2 = 0.5*w2
              Nz1 = floor(z1)
              Nz2 = floor(z2)
              if Ntp == 6:
                  C6 = C1 - Cx - 4*(Nz1 + Nz2)
              else: # 8
                  C8 = C4 - 4*(Nz1 + Nz2)
  
      if Ntp == 1:
          ans = C1
      elif Ntp == 2:
          ans = C2
      elif Ntp == 4:
          ans = C4
      elif Ntp == 6:
          ans = C6
      elif Ntp == 8:
          ans = C8
      else:
          raise Exception('Only 1, 2, 4, 6, or 8 tube passes are supported')
      ans = int(ans)
      # In some cases, a negative number would be returned by these formulas
      if ans < 0:
          ans = 0 # pragma: no cover
      return ans