docs: def nozzle_expansibility(D, Do, P1, P2, k, beta=None):
      r'''Calculates the expansibility factor for a nozzle or venturi nozzle,
      based on the geometry of the plate, measured pressures of the orifice, and
      the isentropic exponent of the fluid.
      
      .. math::
          \epsilon = \left\{\left(\frac{\kappa \tau^{2/\kappa}}{\kappa-1}\right)
          \left(\frac{1 - \beta^4}{1 - \beta^4 \tau^{2/\kappa}}\right)
          \left[\frac{1 - \tau^{(\kappa-1)/\kappa}}{1 - \tau}
          \right] \right\}^{0.5}
          
      Parameters
      ----------
      D : float
          Upstream internal pipe diameter, [m]
      Do : float
          Diameter of orifice of the venturi or nozzle, [m]
      P1 : float
          Static pressure of fluid upstream of orifice at the cross-section of
          the pressure tap, [Pa]
      P2 : float
          Static pressure of fluid downstream of orifice at the cross-section of
          the pressure tap, [Pa]
      k : float
          Isentropic exponent of fluid, [-]
      beta : float, optional
          Optional `beta` ratio, which is useful to specify for wedge meters or
          flow meters which have a different beta ratio calculation, [-]
  
      Returns
      -------
      expansibility : float
          Expansibility factor (1 for incompressible fluids, less than 1 for
          real fluids), [-]
  
      Notes
      -----
      This formula was determined for the range of P2/P1 >= 0.75.
  
      Examples
      --------
      >>> nozzle_expansibility(D=0.0739, Do=0.0222, P1=1E5, P2=9.9E4, k=1.4)
      0.9945702344566746
  
      References
      ----------
      .. [1] American Society of Mechanical Engineers. Mfc-3M-2004 Measurement 
         Of Fluid Flow In Pipes Using Orifice, Nozzle, And Venturi. ASME, 2001.
      .. [2] ISO 5167-3:2003 - Measurement of Fluid Flow by Means of Pressure 
         Differential Devices Inserted in Circular Cross-Section Conduits Running
         Full -- Part 3: Nozzles and Venturi Nozzles.
      '''
      if beta is None:
          beta = Do/D
      beta2 = beta*beta
      beta4 = beta2*beta2
      tau = P2/P1
      term1 = k*tau**(2.0/k )/(k - 1.0)
      term2 = (1.0 - beta4)/(1.0 - beta4*tau**(2.0/k))
      try:
          term3 = (1.0 - tau**((k - 1.0)/k))/(1.0 - tau)
      except ZeroDivisionError:
          '''Obtained with:
              from sympy import *
              tau, k = symbols('tau, k')
              expr = (1 - tau**((k - 1)/k))/(1 - tau)
              limit(expr, tau, 1)
          '''
          term3 = (k - 1.0)/k
      return (term1*term2*term3)**0.5