docs: def diffuser_conical(Di1, Di2, l=None, angle=None, fd=None, Re=None,
                       roughness=0.0, method='Rennels'):
      r'''Returns the loss coefficient for any conical pipe diffuser.
      This calculation has four methods available.
      
      The 'Rennels' [1]_ formulas are as follows (three different formulas are
      used, depending on the angle and the ratio of diameters):
  
      For 0 to 20 degrees, all aspect ratios:
  
      .. math::
          K_1 = 8.30[\tan(\alpha/2)]^{1.75}(1-\beta^2)^2 + \frac{f(1-\beta^4)}{8\sin(\alpha/2)}
  
      For 20 to 60 degrees, beta < 0.5:
  
      .. math::
          K_1 = \left\{1.366\sin\left[\frac{2\pi(\alpha-15^\circ)}{180}\right]^{0.5}
          - 0.170 - 3.28(0.0625-\beta^4)\sqrt{\frac{\alpha-20^\circ}{40^\circ}}\right\}
          (1-\beta^2)^2 + \frac{f(1-\beta^4)}{8\sin(\alpha/2)}
  
      For 20 to 60 degrees, beta >= 0.5:
  
      .. math::
          K_1 = \left\{1.366\sin\left[\frac{2\pi(\alpha-15^\circ)}{180}\right]^{0.5}
          - 0.170 \right\}(1-\beta^2)^2 + \frac{f(1-\beta^4)}{8\sin(\alpha/2)}
  
      For 60 to 180 degrees, beta < 0.5:
  
      .. math::
          K_1 = \left[1.205 - 3.28(0.0625-\beta^4)-12.8\beta^6\sqrt{\frac
          {\alpha-60^\circ}{120^\circ}}\right](1-\beta^2)^2
  
      For 60 to 180 degrees, beta >= 0.5:
  
      .. math::
          K_1 = \left[1.205 - 0.20\sqrt{\frac{\alpha-60^\circ}{120^\circ}}
          \right](1-\beta^2)^2
          
      The Swamee [5]_ formula is:
          
      .. math::
          K = \left\{\frac{0.25}{\theta^3}\left[1 + \frac{0.6}{r^{1.67}}
          \left(\frac{\pi-\theta}{\theta} \right) \right]^{0.533r - 2.6} 
          \right\}^{-0.5}
  
      .. figure:: fittings/diffuser_conical.png
         :scale: 60 %
         :alt: diffuser conical; after [1]_
  
      Parameters
      ----------
      Di1 : float
          Inside diameter of original pipe (smaller), [m]
      Di2 : float
          Inside diameter of following pipe (larger), [m]
      l : float, optional
          Length of the contraction along the pipe axis, optional, [m]
      angle : float, optional
          Angle of contraction, [degrees]
      fd : float, optional
          Darcy friction factor [-]
      Re : float, optional
          Reynolds number of the pipe (used in Rennels method only if no friction
          factor given), [m]
      roughness : float, optional
          Roughness of bend wall (used in Rennel method if no friction factor 
          given), [m]   
      method : str
          The method to use for the calculation; one of 'Rennels', 'Crane',
          'Miller', 'Swamee', or 'Idelchik' [-]
  
      Returns
      -------
      K : float
          Loss coefficient with respect to smaller, upstream diameter [-]
  
      Notes
      -----
      The Miller method changes around quite a bit.
      
      There is quite a bit of variance in the predictions of the methods, as 
      demonstrated by the following figure.
  
      .. plot:: plots/diffuser_conical.py
  
      Examples
      --------
      >>> diffuser_conical(Di1=1/3., Di2=1.0, angle=50.0, Re=1E6)
      0.8027721093415322
  
      References
      ----------
      .. [1] Rennels, Donald C., and Hobart M. Hudson. Pipe Flow: A Practical
         and Comprehensive Guide. 1st edition. Hoboken, N.J: Wiley, 2012.
      .. [2] Idel’chik, I. E. Handbook of Hydraulic Resistance: Coefficients of 
         Local Resistance and of Friction (Spravochnik Po Gidravlicheskim 
         Soprotivleniyam, Koeffitsienty Mestnykh Soprotivlenii i Soprotivleniya
         Treniya). National technical information Service, 1966.
      .. [3] Crane Co. Flow of Fluids Through Valves, Fittings, and Pipe. Crane,
         2009.
      .. [4] Swamee, Prabhata K., and Ashok K. Sharma. Design of Water Supply 
         Pipe Networks. John Wiley & Sons, 2008.
      .. [5] Miller, Donald S. Internal Flow Systems: Design and Performance
         Prediction. Gulf Publishing Company, 1990.
      '''
      beta = Di1/Di2
      beta2 = beta*beta
      if angle is not None:
          angle_rad = radians(angle)
          l = (Di2 - Di1)/(2.0*tan(0.5*angle_rad))
      elif l is not None:
          angle_rad = 2.0*atan(0.5*(Di2-Di1)/l)
          angle = degrees(angle_rad)
      else:
          raise Exception('Either `l` or `angle` must be specified')
      if method is None:
          method == 'Rennels'
      if method == 'Rennels':
          if fd is None:
              if Re is None:
                  raise ValueError("The `Rennels` method requires either a "
                                   "specified friction factor or `Re`")
              fd = Colebrook(Re=Re, eD=roughness/Di2, tol=-1)
          
          if 0.0 < angle <= 20.0:
              K = 8.30*tan(0.5*angle_rad)**1.75*(1.0 - beta2)**2 + 0.125*fd*(1.0 - beta2*beta2)/sin(0.5*angle_rad)
          elif 20.0 < angle <= 60.0 and 0.0 <= beta < 0.5:
              K = (1.366*sin(2.0*pi*(angle - 15.0)/180.)**0.5 - 0.170
              - 3.28*(0.0625-beta**4)*(0.025*(angle-20.0))**0.5)*(1.0 - beta2)**2 + 0.125*fd*(1.0 - beta2*beta2)/sin(0.5*angle_rad)
          elif 20.0 < angle <= 60.0 and beta >= 0.5:
              K = (1.366*sin(2.0*pi*(angle - 15.0)/180.0)**0.5 - 0.170)*(1.0 - beta2)**2 + 0.125*fd*(1.0 - beta2*beta2)/sin(0.5*angle_rad)
          elif 60.0 < angle <= 180.0 and 0.0 <= beta < 0.5:
              beta4 = beta2*beta2
              K = (1.205 - 3.28*(0.0625 - beta4) - 12.8*beta4*beta2*((angle - 60.0)/120.)**0.5)*(1.0 - beta2)**2
          elif 60.0 < angle <= 180.0 and beta >= 0.5:
              K = (1.205 - 0.20*((angle - 60.0)/120.)**0.5)*(1.0 - beta**2)**2
          else:
              raise Exception('Conical diffuser inputs incorrect')
          return K
      elif method == 'Crane':
          return diffuser_conical_Crane(Di1=Di1, Di2=Di2, l=l, angle=angle)
      elif method == 'Miller':
          A_ratio = 1.0/beta2
          if A_ratio > 4.0:
              A_ratio = 4.0
          elif A_ratio < 1.1:
              A_ratio = 1.1
          
          l_R1_ratio = l/(0.5*Di1)
          if l_R1_ratio < 0.1:
              l_R1_ratio = 0.1
          elif l_R1_ratio > 20.0:
              l_R1_ratio = 20.0
          Kd = max(float(bisplev(log(l_R1_ratio), log(A_ratio), tck_diffuser_conical_Miller)), 0)
          return Kd
      elif method == 'Idelchik':
          A_ratio = beta2
          # Angles 0 to 20, ratios 0.05 to 0.06
          if angle > 20.0:
              angle_fric = 20.0
          elif angle < 2.0:
              angle_fric = 2.0
          else:
              angle_fric = angle
              
          A_ratio_fric = A_ratio
          if A_ratio_fric < 0.05:
              A_ratio_fric = 0.05
          elif A_ratio_fric > 0.6:
              A_ratio_fric = 0.6
          
          K_fr = float(contraction_conical_frction_Idelchik_obj(angle_fric, A_ratio_fric))
          K_exp = float(diffuser_conical_Idelchik_obj(min(0.6, A_ratio), max(3.0, angle)))
          return K_fr + K_exp
          
      elif method == 'Swamee':
          # Really starting to thing Swamee uses a different definition of loss coefficient!
          r = Di2/Di1
          K = (0.25*angle_rad**-3*(1.0 + 0.6*r**(-1.67)*(pi-angle_rad)/angle_rad)**(0.533*r - 2.6))**-0.5
          return K
      else:
          raise ValueError('Specified method not recognized; methods are %s'
                           %(diffuser_conical_methods))