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