code: def _saturation(self, T):
          
          rhoc = self._constants.get("rhoref", self.rhoc)
          Tc = self._constants.get("Tref", self.Tc)
  
          if T > Tc:
              T = Tc
          tau = Tc/T
  
          rhoLo = self._Liquid_Density(T)
          rhoGo = self._Vapor_Density(T)
  
          def f(parr):
              rhol, rhog = parr
              deltaL = rhol/rhoc
              deltaG = rhog/rhoc
              phirL = _phir(tau, deltaL, self._constants)
              phirG = _phir(tau, deltaG, self._constants)
              phirdL = _phird(tau, deltaL, self._constants)
              phirdG = _phird(tau, deltaG, self._constants)
              Jl = deltaL*(1+deltaL*phirdL)
              Jv = deltaG*(1+deltaG*phirdG)
              Kl = deltaL*phirdL+phirL+log(deltaL)
              Kv = deltaG*phirdG+phirG+log(deltaG)
              return Kv-Kl, Jv-Jl
  
          rhoL, rhoG = fsolve(f, [rhoLo, rhoGo])
          if rhoL == rhoG:
              Ps = self.Pc
          else:
              deltaL = rhoL/self.rhoc
              deltaG = rhoG/self.rhoc
              firL = _phir(tau, deltaL, self._constants)
              firG = _phir(tau, deltaG, self._constants)
  
              Ps = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(firL-firG+log(deltaL/deltaG))
          return rhoL, rhoG, Ps
docs: Saturation calculation for two phase search