docs: def P_upstream_isothermal_critical_flow(P, fd, D, L):
      '''Not part of the public API. Reverses `P_isothermal_critical_flow`.
  
      Examples
      --------
      >>> P_upstream_isothermal_critical_flow(P=389699.7317645518, fd=0.00185,
      ... L=1000., D=0.5)
      1000000.0000000001
      '''
      lambertw_term = float(lambertw(-exp(-(fd*L+D)/D), -1).real)
      return exp(-0.5*(D*lambertw_term+fd*L+D)/D)*P