docs: def P_isothermal_critical_flow(P, fd, D, L):
      r'''Calculates critical flow pressure `Pcf` for a fluid flowing
      isothermally and suffering pressure drop caused by a pipe's friction factor.
  
      .. math::
          P_2 = P_{1} e^{\frac{1}{2 D} \left(D \left(\operatorname{LambertW}
          {\left (- e^{\frac{1}{D} \left(- D - L f_d\right)} \right )} + 1\right)
          + L f_d\right)}
  
      Parameters
      ----------
      P : float
          Inlet pressure [Pa]
      fd : float
          Darcy friction factor for flow in pipe [-]
      D : float
          Diameter of pipe, [m]
      L : float
          Length of pipe, [m]
  
      Returns
      -------
      Pcf : float
          Critical flow pressure of a compressible gas flowing from `P1` to `Pcf`
          in a tube of length L and friction factor `fd` [Pa]
  
      Notes
      -----
      Assumes isothermal flow. Developed based on the `isothermal_gas` model,
      using SymPy.
  
      The isothermal gas model is solved for maximum mass flow rate; any pressure
      drop under it is impossible due to the formation of a shock wave.
  
      Examples
      --------
      >>> P_isothermal_critical_flow(P=1E6, fd=0.00185, L=1000., D=0.5)
      389699.7317645518
  
      References
      ----------
      .. [1] Wilkes, James O. Fluid Mechanics for Chemical Engineers with
         Microfluidics and CFD. 2 edition. Upper Saddle River, NJ: Prentice Hall,
         2005.
      '''
      # Correct branch of lambertw found by trial and error
      lambert_term = float(lambertw(-exp((-D - L*fd)/D), -1).real)
      return P*exp((D*(lambert_term + 1.0) + L*fd)/(2.0*D))