def nozzle_expansibility(D, Do, P1, P2, k, beta=None):
    """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
