docs: def bend_rounded_Miller(Di, angle, Re, rc=None, bend_diameters=None, 
                          roughness=0.0, L_unimpeded=None):
      r'''Calculates the loss coefficient for a rounded pipe bend according to 
      Miller [1]_. This is a sophisticated model which uses corrections for
      pipe roughness, the length of the pipe downstream before another 
      interruption, and a correction for Reynolds number. It interpolates several
      times using several corrections graphs in [1]_.
      
      Parameters
      ----------
      Di : float
          Inside diameter of pipe, [m]
      angle : float
          Angle of bend, [degrees]
      Re : float
          Reynolds number of the pipe (no specification if inlet or outlet
          properties should be used), [m]
      rc : float, optional
          Radius of curvature of the entrance, [m]
      bend_diameters : float, optional
          Number of diameters of pipe making up the bend radius  (used if rc not 
          provided; defaults to 5), [-]
      roughness : float, optional
          Roughness of bend wall, [m]   
      L_unimpeded : float, optional
          The length of unimpeded pipe without any fittings, instrumentation,
          or flow disturbances downstream (assumed 20 diameters if not 
          specified), [m]
  
      Returns
      -------
      K : float
          Loss coefficient [-]
  
      Notes
      -----
      When inputting bend diameters, note that manufacturers often specify
      this as a multiplier of nominal diameter, which is different than actual
      diameter. Those require that rc be specified.
      
      `rc` is limited to 0.5 or above; which represents a sharp, square, inner 
      edge - and an outer bend radius of 1.0. Losses are at a minimum when this
      value is large.
      
      This was developed for bend angles between 10 and 180 degrees; and r/D
      ratios between 0.5 and 10. Both smooth and rough data was used in its 
      development from several sources.
      
      Note the loss coefficient includes the surface friction of the pipe as if
      it was straight.
     
      Examples
      --------
      >>> bend_rounded_Miller(Di=.6, bend_diameters=2, angle=90,  Re=2e6,
      ... roughness=2E-5, L_unimpeded=30*.6)
      0.152618207051459
      
      References
      ----------
      .. [1] Miller, Donald S. Internal Flow Systems: Design and Performance
         Prediction. Gulf Publishing Company, 1990.
      '''
      if not rc:
          if bend_diameters is None:
              bend_diameters = 5
          rc = Di*bend_diameters
      
      radius_ratio = rc/Di
      
      if L_unimpeded is None:
          # Assumption - smooth outlet
          L_unimpeded = 20*Di
      
      # Graph is defined for angles 10 to 180 degrees, ratios 0.5 to 10
      if radius_ratio < 0.5:
          radius_ratio = 0.5
      if radius_ratio > 10.0:
          radius_ratio = 10.0
      if angle < 10.0:
          angle = 10.0
      
      # Curve fit in terms of degrees
      # Caching could work here - angle, radius ratio does not change often
      Kb = bend_rounded_Miller_Kb(radius_ratio, angle) 
      
      C_roughness = Miller_bend_roughness_correction(Re=Re, Di=Di, 
                                                     roughness=roughness)
  
      '''Section 9.2.2 - Reynolds Number Correction
      Allow some extrapolation up to 1E8 (1E7 max in graph but the trend looks good)
      '''
      Re_C_Re = min(max(Re, 1E4), 1E8)
      if radius_ratio >= 2.0:
          if Re_C_Re == 1E8:
              C_Re = 0.4196741237602154 # bend_rounded_Miller_C_Re(1e8, 2.0)
          elif Re_C_Re == 1E4:
              C_Re = 2.1775876405173977 # bend_rounded_Miller_C_Re(1e4, 2.0)
          else:
              C_Re = bend_rounded_Miller_C_Re(Re_C_Re, 2.0)
      elif radius_ratio <= 1.0:
          # newton(lambda x: bend_rounded_Miller_C_Re(x, 1.0)-1, 2e5) to get the boundary value
          C_Re_1 = bend_rounded_Miller_C_Re(Re_C_Re, 1.0) if Re_C_Re < 207956.58904584477 else 1.0
          if radius_ratio > 0.7 or Kb < 0.4:
              C_Re = C_Re_1
          else:
              C_Re = Kb/(Kb - 0.2*C_Re_1 + 0.2)
              if C_Re > 2.2 or C_Re < 0:
                  C_Re = 2.2
      else:
          # regardless of ratio - 1
          if Re_C_Re > 1048884.4656835075:
              C_Re = 1.0
          elif Re_C_Re > horner(bend_rounded_Miller_C_Re_limit_1, radius_ratio):
              C_Re = 1.0
  #            ps = np.linspace(1, 2)
  #            qs = [newton(lambda x: bend_rounded_Miller_C_Re(x, i)-1, 2e5) for i in ps]
  #            np.polyfit(ps, qs, 4).tolist()
              # Line of C_Re=1 as a function of r_d between 0 and 1
          else:
              C_Re = bend_rounded_Miller_C_Re(Re_C_Re, radius_ratio)
      C_o =  Miller_bend_unimpeded_correction(Kb=Kb, Di=Di, L_unimpeded=L_unimpeded)
      
  #    print('Kb=%g, C Re=%g, C rough =%g, Co=%g' %(Kb, C_Re, C_roughness, C_o))
      return Kb*C_Re*C_roughness*C_o