Calculate U239 fission product decay heat power from ANS-5.1-1979 Eq.14
 ts: Time after remove (sec)
 t0: Cumulative reactor operating time (sec)

 calc_thermal_fission_functions_with_U239(t0, ts, read_data)

 return fU239
176042
def calc_thermal_fission_functions_with_U239(t0, ts, read_data)
      read_data.theEU239 * read_data.theR * (1 - Math.exp(-read_data.lamda1 * t0)) * Math.exp(-read_data.lamda1 * ts)
    end
832236
def thePi
      hash = HashWithThermalFission.new

      hash.thermal_fission[:U235]  = 0.98
      hash.thermal_fission[:Pu239] = 0.01
      hash.thermal_fission[:U238]  = 0.01

      hash.thermal_fission
    end
920164
def eta_from_tau0_tau3(tau0, tau3, f_lower):
    """
    mtotal = mtotal_from_tau0_tau3(tau0, tau3, f_lower,
                                   in_seconds=True)
    eta = mtotal**(-2./3.) * (_a3(f_lower) / tau3)
    return eta
1236783
def calculate_perf_counter_bulk_count(previous, current, property_name):
    
    n0 = previous[property_name]
    n1 = current[property_name]
    d0 = previous["Timestamp_Sys100NS"]
    d1 = current["Timestamp_Sys100NS"]
    f = current["Frequency_Sys100NS"]

    if n0 is None or n1 is None:
        return

    return (n1 - n0) / ((d1 - d0) / f)
1245322
def tau0_from_mtotal_eta(mtotal, eta, f_lower):
    
    """
    # convert to seconds
    mtotal = mtotal * lal.MTSUN_SI
    # formulae from arxiv.org:0706.4437
    return _a0(f_lower) / (mtotal**(5./3.) * eta)
1005033
def Pgas(rho,T,mu):
    ''' 
    

    '''
    R = old_div(boltzmann_constant, atomic_mass_unit)
    return (old_div(R,mu)) * rho * T
225853
def ultimate_oscillator(close_data, low_data):
    
    a7 = 4 * average_7(close_data, low_data)
    a14 = 2 * average_14(close_data, low_data)
    a28 = average_28(close_data, low_data)
    uo = 100 * ((a7 + a14 + a28) / 7)
    return uo
1223539
def T(t, step=30):
    
    if hasattr(t, 'timestamp'):
        timestamp = t.timestamp()
    else:
        # python 2
        if t.tzinfo is None:
            timestamp = time.mktime(t.timetuple())
        else:
            utc_naive  = t.replace(tzinfo=None) - t.utcoffset()
            timestamp = (utc_naive - datetime.datetime(1970, 1, 1)).total_seconds()

    return int(timestamp) // step
286434
def debye_temperature(self, structure):
        
        v0 = (structure.volume * 1e-30 / structure.num_sites)
        vl, vt = self.long_v(structure), self.trans_v(structure)
        vm = 3**(1./3.) * (1 / vl**3 + 2 / vt**3)**(-1./3.)
        td = 1.05457e-34 / 1.38065e-23 * vm * (6 * np.pi**2 / v0) ** (1./3.)
        return td
855504
def dec2hms(dec):
    
    DEGREE = 360.
    HOUR = 24.
    MINUTE = 60.
    SECOND = 3600.
    
    dec = float(dec)
    fhour = dec*(HOUR/DEGREE)
    hour = int(fhour)

    fminute = (fhour - hour)*MINUTE
    minute = int(fminute)
    
    second = (fminute - minute)*MINUTE
    return (hour, minute, second)
1031282
def get_abundance_iso_decay(self,cycle):
        

        import nuutils as u

        masses_for_this_cycle = self.se.get(cycle,'mass')
        self._read_iso_abund_marco([min(masses_for_this_cycle),max(masses_for_this_cycle)],cycle)

        u.stable_specie()
        self.decay(self.mass_frac)

        self.index_for_all_species = u.cl
        self.index_for_stable_species = u.back_ind

        self.decayed_stable_isotopes_per_cycle = decayed_multi_d

        # from here read solar abundances
        solar_factor = 2.
        u.solar('iniab1.0E-02.ppn_GN93',solar_factor)

        self.stable_isotope_identifier=u.jjdum
        self.stable_isotope_list=u.stable

        self.isotopic_production_factors=[]
        for i in range(len(masses_for_this_cycle)):
            pf_dum=[]
            jj=0
            for j in range(len(self.stable_isotope_identifier)):
                if self.stable_isotope_identifier[j] == 1:
                    pf_dum.append(float(old_div(self.mass_frac[i][self.index_for_all_species[self.stable_isotope_list
[jj].capitalize()]],u.solar_abundance[self.stable_isotope_list[jj].lower()])))
                    jj=jj+1
                #elif self.stable_isotope_identifier[j] == 0:
                #       pf_dum.append(float(0.))
            self.isotopic_production_factors.append(pf_dum)

        self.isotopic_production_factors_decayed=[]
        for i in range(len(masses_for_this_cycle)):
            pf_dum_d=[]
            jj=0
            for j in range(len(self.stable_isotope_identifier)):
                if self.stable_isotope_identifier[j] == 1:
                    pf_dum_d.append(float(old_div(self.decayed_stable_isotopes_per_cycle[i][self.index_for_stable_species[self.stable_isotope_list
[jj].upper()]],u.solar_abundance[self.stable_isotope_list[jj].lower()])))
                    jj=jj+1
            self.isotopic_production_factors_decayed.append(pf_dum_d)
317964
def calc_uncorrelated_signal_power(params):
    
    ''''''


    for i, data_type in enumerate(['LFP','CSD']):
        if i % SIZE == RANK:
    
            # Determine size of PSD matrix
    
            f = h5py.File(os.path.join(params.savefolder, data_type + 'sum.h5'),'r')
            data = f['data'].value
            srate = f['srate'].value
            if ana_params.mlab:
                Psum, freqs = plt.mlab.psd(data[0], NFFT=ana_params.NFFT, Fs=srate,
                                    noverlap=ana_params.noverlap, window=ana_params.window)
            else:
                [freqs, Psum] = hlp.powerspec([data[0]], tbin= 1./srate*1000.,
                                       Df=ana_params.Df, pointProcess=False)
            f.close()
            P = np.zeros((data.shape[0],Psum.shape[0]))
               
            for y in params.y:
    
                print 'processing ', y
                # Load data
                f = h5py.File(os.path.join(params.populations_path, '%s_%ss' %
                                           (y,data_type) + '.h5'),'r')
                data_y = f['data'].value[:,:, ana_params.transient:]
                # subtract mean
                for j in range(len(data_y)):
                    data_yT = data_y[j].T - data_y[j].mean(axis=1)
                    data_y[j] = data_yT.T
                srate = f['srate'].value
                tvec = np.arange(data_y.shape[2]) * 1000. / srate
                f.close()
    
            
                for j in range(len(data_y)): # loop over cells
                    if ana_params.mlab:
                        for ch in range(len(params.electrodeParams['z'])): # loop over channels
                            P_j_ch, freqs = plt.mlab.psd(data_y[j,ch],
                                                         NFFT=ana_params.NFFT, Fs=srate,
                                                         noverlap=ana_params.noverlap,
                                                         window=ana_params.window)
                                       
                            P[ch] += P_j_ch
                    else:
                        [freqs, P_j] = hlp.powerspec(data_y[j], tbin= 1./srate*1000.,
                                                     Df=ana_params.Df, pointProcess=False)
                        mask = np.where(freqs >= 0.)
                        freqs = freqs[mask]
                        P_j = P_j[:,mask][:,0,:]
                        P_j = P_j/tvec[tvec >= ana_params.transient].size**2
                        
                        P += P_j
    
            #rescale PSD as they may be computed from a fraction of single cell LFPs
            P /= params.recordSingleContribFrac
            
            # Save data
            f = h5py.File(os.path.join(params.savefolder, ana_params.analysis_folder,
                                       data_type +  ana_params.fname_psd_uncorr),'w')
            f['freqs']=freqs
            f['psd']=P
            f['transient']=ana_params.transient
            f['mlab']=ana_params.mlab
            f['NFFT']=ana_params.NFFT
            f['noverlap']=ana_params.noverlap
            f['window']=str(ana_params.window)
            f['Df']=str(ana_params.Df)        
            f.close()

    return
328003
def mw(mol, ndigits=2):
    
    mol.require("Valence")
    return round(sum(a.mw() for _, a in mol.atoms_iter()), ndigits)
33440
def _decay(self):
        """"""
        costs = []
        for var in tf.trainable_variables():
            if var.op.name.find(r"DW") > 0:
                costs.append(tf.nn.l2_loss(var))

        return tf.multiply(self.hps.weight_decay_rate, tf.add_n(costs))
191668
def _compute_heating_rates(self):
        ''''''
        (ncol, nlay, icld, iaer, permuteseed, irng,
         play, plev, tlay, tlev, tsfc,
         h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
         aldif, aldir, asdif, asdir, coszen, adjes, dyofyr, scon, isolvar,
         indsolvar, bndsolvar, solcycfrac,
         inflgsw, iceflgsw, liqflgsw,
         cldfrac, ciwp, clwp, reic, relq, tauc, ssac, asmc, fsfc,
         tauaer, ssaaer, asmaer, ecaer,) = self._prepare_sw_arguments()
        if icld == 0:  # clear-sky only
            cldfmcl = np.zeros((ngptsw,ncol,nlay))
            ciwpmcl = np.zeros((ngptsw,ncol,nlay))
            clwpmcl = np.zeros((ngptsw,ncol,nlay))
            reicmcl = np.zeros((ncol,nlay))
            relqmcl = np.zeros((ncol,nlay))
            taucmcl = np.zeros((ngptsw,ncol,nlay))
            ssacmcl = np.zeros((ngptsw,ncol,nlay))
            asmcmcl = np.zeros((ngptsw,ncol,nlay))
            fsfcmcl = np.zeros((ngptsw,ncol,nlay))
        else:
            #  Call the Monte Carlo Independent Column Approximation (McICA, Pincus et al., JC, 2003)
            (cldfmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl, taucmcl,
            ssacmcl, asmcmcl, fsfcmcl) = _rrtmg_sw.climlab_mcica_subcol_sw(
                            ncol, nlay, icld, permuteseed, irng, play,
                            cldfrac, ciwp, clwp, reic, relq, tauc, ssac, asmc, fsfc)
        #  Call the RRTMG_SW driver to compute radiative fluxes
        (swuflx, swdflx, swhr, swuflxc, swdflxc, swhrc) = \
            _rrtmg_sw.climlab_rrtmg_sw(ncol, nlay, icld, iaer,
                play, plev, tlay, tlev, tsfc,
                h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
                asdir, asdif, aldir, aldif,
                coszen, adjes, dyofyr, scon, isolvar,
                inflgsw, iceflgsw, liqflgsw, cldfmcl,
                taucmcl, ssacmcl, asmcmcl, fsfcmcl,
                ciwpmcl, clwpmcl, reicmcl, relqmcl,
                tauaer, ssaaer, asmaer, ecaer,
                bndsolvar, indsolvar, solcycfrac)
        #  Output is all (ncol,nlay+1) or (ncol,nlay)
        self.SW_flux_up = _rrtm_to_climlab(swuflx) + 0.*self.SW_flux_up
        self.SW_flux_down = _rrtm_to_climlab(swdflx) + 0.*self.SW_flux_down
        self.SW_flux_up_clr = _rrtm_to_climlab(swuflxc) + 0.*self.SW_flux_up_clr
        self.SW_flux_down_clr = _rrtm_to_climlab(swdflxc) + 0.*self.SW_flux_down_clr
        #  Compute quantities derived from fluxes, including ASR
        self._compute_SW_flux_diagnostics()
        #  calculate heating rates from flux divergence
        SWheating_Wm2 = np.array(-np.diff(self.SW_flux_net, axis=-1)) + 0.*self.Tatm
        SWheating_clr_Wm2 = np.array(-np.diff(self.SW_flux_net_clr, axis=-1)) + 0.*self.Tatm
        self.heating_rate['Ts'] = np.array(self.SW_flux_net[..., -1, np.newaxis]) + 0.*self.Ts
        self.heating_rate['Tatm'] = SWheating_Wm2
        #  Convert to K / day
        Catm = self.Tatm.domain.heat_capacity
        self.TdotSW = SWheating_Wm2 / Catm * const.seconds_per_day
        self.TdotSW_clr = SWheating_clr_Wm2 / Catm * const.seconds_per_day
284869
def tau3_from_mtotal_eta(mtotal, eta, f_lower):
    
    """
    # convert to seconds
    mtotal = mtotal * lal.MTSUN_SI
    # formulae from arxiv.org:0706.4437
    return _a3(f_lower) / (mtotal**(2./3.) * eta)
1187092
def calculo(self):
        """"""
        T = self.kwargs["T"]
        rho = self.kwargs["rho"]
        P = self.kwargs["P"]
        s = self.kwargs["s"]
        h = self.kwargs["h"]
        u = self.kwargs["u"]
        x = self.kwargs["x"]

        # Initial values
        T0 = self.kwargs["T0"]
        rho0 = self.kwargs["rho0"]

        if T0 or rho0:
            To = T0
            rhoo = rho0
        elif self.name == "air":
            To = 300
            rhoo = 1e-3
        else:
            try:
                st0 = IAPWS97(**self.kwargs)
            except NotImplementedError:
                To = 300
                rhoo = 900
            else:
                if st0.status:
                    To = st0.T
                    rhoo = st0.rho
                else:
                    To = 300
                    rhoo = 900

        self.R = self._constants["R"]/self._constants.get("M", self.M)
        rhoc = self._constants.get("rhoref", self.rhoc)
        Tc = self._constants.get("Tref", self.Tc)

        propiedades = None

        if self._mode not in ("Tx", "Px"):
            # Method with iteration necessary to get x
            if self._mode == "TP":
                try:
                    if self.name == "air":
                        raise ValueError
                    st0 = IAPWS97(**self.kwargs)
                    rhoo = st0.rho
                except NotImplementedError:
                    if rho0:
                        rhoo = rho0
                    elif T < self.Tc and P < self.Pc and \
                            self._Vapor_Pressure(T) < P:
                        rhoo = self._Liquid_Density(T)
                    elif T < self.Tc and P < self.Pc:
                        rhoo = self._Vapor_Density(T)
                    else:
                        rhoo = self.rhoc*3
                except ValueError:
                    rhoo = 1e-3

                def f(rho):
                    delta = rho/rhoc
                    tau = Tc/T

                    fird = _phird(tau, delta, self._constants)
                    Po = (1+delta*fird)*self.R*T*rho
                    return Po-P*1000

                rho = fsolve(f, rhoo)[0]

                # Calculate quality
                if T > self.Tc:
                    x = 1
                else:
                    Ps = self._Vapor_Pressure(T)
                    if Ps*0.95 < P < Ps*1.05:
                        rhol, rhov, Ps = self._saturation(T)
                        Ps *= 1e-3

                    if Ps > P:
                        x = 1
                    else:
                        x = 0

            elif self._mode == "Th":
                tau = Tc/T
                ideal = self._phi0(tau, 1)
                fiot = ideal["fiot"]

                def f(rho):
                    delta = rho/rhoc
                    fird = _phird(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
                    return ho-h

                if T >= self.Tc:
                    rhoo = self.rhoc
                    rho = fsolve(f, rhoo)[0]
                else:
                    x0 = self.kwargs["x0"]
                    rhov = self._Vapor_Density(T)
                    rhol = self._Liquid_Density(T)
                    deltaL = rhol/rhoc
                    deltaG = rhov/rhoc

                    firdL = _phird(tau, deltaL, self._constants)
                    firtL = _phirt(tau, deltaL, self._constants)
                    firdG = _phird(tau, deltaG, self._constants)
                    firtG = _phirt(tau, deltaG, self._constants)
                    hl = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
                    hv = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
                    if x0 not in (0, 1) and hl <= h <= hv:
                        rhol, rhov, Ps = self._saturation(T)
                        vapor = self._Helmholtz(rhov, T)
                        liquido = self._Helmholtz(rhol, T)
                        hv = vapor["h"]
                        hl = liquido["h"]
                        x = (h-hl)/(hv-hl)
                        rho = 1/(x/rhov+(1-x)/rhol)
                        P = Ps/1000
                    else:
                        if h > hv:
                            rhoo = rhov
                        else:
                            rhoo = rhol
                        rho = fsolve(f, rhoo)[0]

            elif self._mode == "Ts":
                tau = Tc/T

                def f(rho):
                    if rho < 0:
                        rho = 1e-20
                    delta = rho/rhoc

                    ideal = self._phi0(tau, delta)
                    fio = ideal["fio"]
                    fiot = ideal["fiot"]
                    fir = _phir(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    so = self.R*(tau*(fiot+firt)-fio-fir)
                    return so-s

                if T >= self.Tc:
                    rhoo = self.rhoc
                    rho = fsolve(f, rhoo)[0]
                else:
                    rhov = self._Vapor_Density(T)
                    rhol = self._Liquid_Density(T)
                    deltaL = rhol/rhoc
                    deltaG = rhov/rhoc

                    idealL = self._phi0(tau, deltaL)
                    idealG = self._phi0(tau, deltaG)
                    fioL = idealL["fio"]
                    fioG = idealG["fio"]
                    fiot = idealL["fiot"]
                    firL = _phir(tau, deltaL, self._constants)
                    firtL = _phirt(tau, deltaL, self._constants)
                    sl = self.R*(tau*(fiot+firtL)-fioL-firL)
                    firG = _phir(tau, deltaG, self._constants)
                    firtG = _phirt(tau, deltaG, self._constants)
                    sv = self.R*(tau*(fiot+firtG)-fioG-firG)

                    if sl <= s <= sv:
                        rhol, rhov, Ps = self._saturation(T)
                        vapor = self._Helmholtz(rhov, T)
                        liquido = self._Helmholtz(rhol, T)
                        sv = vapor["s"]
                        sl = liquido["s"]
                        x = (s-sl)/(sv-sl)
                        rho = 1/(x/rhov+(1-x)/rhol)
                        P = Ps/1000
                    else:
                        if s > sv:
                            rhoo = rhov
                        else:
                            rhoo = rhol
                        rho = fsolve(f, rhoo)[0]

            elif self._mode == "Tu":
                tau = Tc/T
                ideal = self._phi0(tau, 1)
                fiot = ideal["fiot"]

                def f(rho):
                    delta = rho/rhoc

                    fird = _phird(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    Po = (1+delta*fird)*self.R*T*rho
                    ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)

                    return ho-Po/rho-u

                if T >= self.Tc:
                    rhoo = self.rhoc
                    rho = fsolve(f, rhoo)[0]
                else:
                    rhov = self._Vapor_Density(T)
                    rhol = self._Liquid_Density(T)
                    deltaL = rhol/rhoc
                    deltaG = rhov/rhoc

                    firdL = _phird(tau, deltaL, self._constants)
                    firtL = _phirt(tau, deltaL, self._constants)
                    firdG = _phird(tau, deltaG, self._constants)
                    firtG = _phirt(tau, deltaG, self._constants)
                    PoL = (1+deltaL*firdL)*self.R*T*rhol
                    PoG = (1+deltaG*firdG)*self.R*T*rhov
                    hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
                    hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)

                    uv = hoG-PoG/rhov
                    ul = hoL-PoL/rhol
                    if ul <= u <= uv:
                        rhol, rhov, Ps = self._saturation(T)
                        vapor = self._Helmholtz(rhov, T)
                        liquido = self._Helmholtz(rhol, T)
                        uv = vapor["h"]-vapor["P"]/rhov
                        ul = liquido["h"]-liquido["P"]/rhol
                        x = (u-ul)/(uv-ul)
                        rho = 1/(x/rhov-(1-x)/rhol)
                        P = Ps/1000
                    else:
                        if u > uv:
                            rhoo = rhov
                        else:
                            rhoo = rhol
                        rho = fsolve(f, rhoo)[0]

            elif self._mode == "Prho":
                delta = rho/rhoc

                def f(T):
                    tau = Tc/T

                    fird = _phird(tau, delta, self._constants)
                    Po = (1+delta*fird)*self.R*T*rho
                    return Po-P*1000

                T = fsolve(f, To)[0]
                rhol = self._Liquid_Density(T)
                rhov = self._Vapor_Density(T)
                if T == To or rhov <= rho <= rhol:

                    def f(parr):
                        T, rhol, rhog = parr
                        tau = Tc/T
                        deltaL = rhol/self.rhoc
                        deltaG = rhog/self.rhoc

                        firL = _phir(tau, deltaL, self._constants)
                        firdL = _phird(tau, deltaL, self._constants)
                        firG = _phir(tau, deltaG, self._constants)
                        firdG = _phird(tau, deltaG, self._constants)

                        Jl = rhol*(1+deltaL*firdL)
                        Jv = rhog*(1+deltaG*firdG)
                        K = firL-firG
                        Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
                        return (Jl-Jv,
                                Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
                                Ps - P*1000)

                    for to in [To, 300, 400, 500, 600]:
                        rhoLo = self._Liquid_Density(to)
                        rhoGo = self._Vapor_Density(to)
                        sol = fsolve(f, [to, rhoLo, rhoGo], full_output=True)
                        T, rhoL, rhoG = sol[0]
                        x = (1./rho-1/rhoL)/(1/rhoG-1/rhoL)
                        if sol[2] == 1 and 0 <= x <= 1 and \
                                sum(abs(sol[1]["fvec"])) < 1e-5:
                            break

                    if sum(abs(sol[1]["fvec"])) > 1e-5:
                        raise(RuntimeError(sol[3]))

                    liquido = self._Helmholtz(rhoL, T)
                    vapor = self._Helmholtz(rhoG, T)
                    P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
                        liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000

            elif self._mode == "Ph":
                def funcion(parr):
                    rho, T = parr
                    delta = rho/rhoc
                    tau = Tc/T

                    ideal = self._phi0(tau, delta)
                    fiot = ideal["fiot"]
                    fird = _phird(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    Po = (1+delta*fird)*self.R*T*rho
                    ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
                    return Po-P*1000, ho-h

                rho, T = fsolve(funcion, [rhoo, To])
                rhol = self._Liquid_Density(T)
                rhov = self._Vapor_Density(T)
                if rho == rhoo or rhov <= rho <= rhol:

                    def f(parr):
                        T, rhol, rhog, x = parr
                        tau = Tc/T
                        deltaL = rhol/self.rhoc
                        deltaG = rhog/self.rhoc

                        ideal = self._phi0(tau, deltaL)
                        fiot = ideal["fiot"]

                        firL = _phir(tau, deltaL, self._constants)
                        firdL = _phird(tau, deltaL, self._constants)
                        firtL = _phirt(tau, deltaL, self._constants)
                        hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
                        firG = _phir(tau, deltaG, self._constants)
                        firdG = _phird(tau, deltaG, self._constants)
                        firtG = _phirt(tau, deltaG, self._constants)
                        hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)

                        Jl = rhol*(1+deltaL*firdL)
                        Jv = rhog*(1+deltaG*firdG)
                        K = firL-firG
                        Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))

                        return (Jl-Jv,
                                Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
                                hoL*(1-x)+hoG*x - h,
                                Ps - P*1000)

                    for to in [To, 300, 400, 500, 600]:
                        rLo = self._Liquid_Density(to)
                        rGo = self._Vapor_Density(to)
                        sol = fsolve(f, [to, rLo, rGo, 0.5], full_output=True)
                        T, rhoL, rhoG, x = sol[0]
                        if sol[2] == 1 and 0 <= x <= 1 and \
                                sum(abs(sol[1]["fvec"])) < 1e-5:
                            break

                    if sum(abs(sol[1]["fvec"])) > 1e-5:
                        raise(RuntimeError(sol[3]))

                    liquido = self._Helmholtz(rhoL, T)
                    vapor = self._Helmholtz(rhoG, T)
                    P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
                        liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000

            elif self._mode == "Ps":
                try:
                    x0 = st0.x
                except NameError:
                    x0 = None

                if x0 is None or x0 == 0 or x0 == 1:
                    def f(parr):
                        rho, T = parr
                        delta = rho/rhoc
                        tau = Tc/T

                        ideal = self._phi0(tau, delta)
                        fio = ideal["fio"]
                        fiot = ideal["fiot"]
                        fird = _phird(tau, delta, self._constants)
                        fir = _phir(tau, delta, self._constants)
                        firt = _phirt(tau, delta, self._constants)
                        Po = (1+delta*fird)*self.R*T*rho
                        so = self.R*(tau*(fiot+firt)-fio-fir)
                        return Po-P*1000, so-s

                    rho, T = fsolve(f, [rhoo, To])

                else:
                    def funcion(parr):
                        rho, T = parr
                        rhol, rhov, Ps = self._saturation(T)
                        vapor = self._Helmholtz(rhov, T)
                        liquido = self._Helmholtz(rhol, T)
                        x = (1./rho-1/rhol)/(1/rhov-1/rhol)
                        return Ps-P*1000, vapor["s"]*x+liquido["s"]*(1-x)-s
                    rho, T = fsolve(funcion, [2., 500.])
                    rhol, rhov, Ps = self._saturation(T)
                    vapor = self._Helmholtz(rhov, T)
                    liquido = self._Helmholtz(rhol, T)
                    sv = vapor["s"]
                    sl = liquido["s"]
                    x = (s-sl)/(sv-sl)

            elif self._mode == "Pu":
                def f(parr):
                    rho, T = parr
                    delta = rho/rhoc
                    tau = Tc/T

                    ideal = self._phi0(tau, delta)
                    fiot = ideal["fiot"]
                    fird = _phird(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    Po = (1+delta*fird)*self.R*T*rho
                    ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
                    return ho-Po/rho-u, Po-P*1000

                sol = fsolve(f, [rhoo, To], full_output=True)
                rho, T = sol[0]
                rhol = self._Liquid_Density(T)
                rhov = self._Vapor_Density(T)
                if rho == rhoo or sol[2] != 1:

                    def f(parr):
                        T, rhol, rhog, x = parr
                        tau = Tc/T
                        deltaL = rhol/self.rhoc
                        deltaG = rhog/self.rhoc

                        ideal = self._phi0(tau, deltaL)
                        fiot = ideal["fiot"]

                        firL = _phir(tau, deltaL, self._constants)
                        firdL = _phird(tau, deltaL, self._constants)
                        firtL = _phirt(tau, deltaL, self._constants)
                        hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
                        firG = _phir(tau, deltaG, self._constants)
                        firdG = _phird(tau, deltaG, self._constants)
                        firtG = _phirt(tau, deltaG, self._constants)
                        hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)

                        Jl = rhol*(1+deltaL*firdL)
                        Jv = rhog*(1+deltaG*firdG)
                        K = firL-firG
                        Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
                        vu = hoG-Ps/rhog
                        lu = hoL-Ps/rhol
                        return (Jl-Jv,
                                Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
                                lu*(1-x)+vu*x - u,
                                Ps - P*1000)

                    for to in [To, 300, 400, 500, 600]:
                        rLo = self._Liquid_Density(to)
                        rGo = self._Vapor_Density(to)
                        sol = fsolve(f, [to, rLo, rGo, 0.5], full_output=True)
                        T, rhoL, rhoG, x = sol[0]
                        if sol[2] == 1 and 0 <= x <= 1 and \
                                sum(abs(sol[1]["fvec"])) < 1e-5:
                            break

                    if sum(abs(sol[1]["fvec"])) > 1e-5:
                        raise(RuntimeError(sol[3]))

                    liquido = self._Helmholtz(rhoL, T)
                    vapor = self._Helmholtz(rhoG, T)
                    P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
                        liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000

            elif self._mode == "rhoh":
                delta = rho/rhoc

                def f(T):
                    tau = Tc/T

                    ideal = self._phi0(tau, delta)
                    fiot = ideal["fiot"]
                    fird = _phird(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
                    return ho-h

                T = fsolve(f, To)[0]
                rhol = self._Liquid_Density(T)
                rhov = self._Vapor_Density(T)
                if T == To or rhov <= rho <= rhol:
                    def f(parr):
                        T, rhol, rhog = parr
                        tau = Tc/T
                        deltaL = rhol/self.rhoc
                        deltaG = rhog/self.rhoc

                        ideal = self._phi0(tau, deltaL)
                        fiot = ideal["fiot"]
                        firL = _phir(tau, deltaL, self._constants)
                        firdL = _phird(tau, deltaL, self._constants)
                        firtL = _phirt(tau, deltaL, self._constants)
                        hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
                        firG = _phir(tau, deltaG, self._constants)
                        firdG = _phird(tau, deltaG, self._constants)
                        firtG = _phirt(tau, deltaG, self._constants)
                        hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)

                        Jl = rhol*(1+deltaL*firdL)
                        Jv = rhog*(1+deltaG*firdG)
                        K = firL-firG
                        x = (1./rho-1/rhol)/(1/rhog-1/rhol)
                        return (Jl-Jv,
                                Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
                                hoL*(1-x)+hoG*x - h)

                    for to in [To, 300, 400, 500, 600]:
                        rhoLo = self._Liquid_Density(to)
                        rhoGo = self._Vapor_Density(to)
                        sol = fsolve(f, [to, rhoLo, rhoGo], full_output=True)
                        T, rhoL, rhoG = sol[0]
                        x = (1./rho-1/rhoL)/(1/rhoG-1/rhoL)
                        if sol[2] == 1 and 0 <= x <= 1 and \
                                sum(abs(sol[1]["fvec"])) < 1e-5:
                            break

                    if sum(abs(sol[1]["fvec"])) > 1e-5:
                        raise(RuntimeError(sol[3]))

                    liquido = self._Helmholtz(rhoL, T)
                    vapor = self._Helmholtz(rhoG, T)
                    P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
                        liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000

            elif self._mode == "rhos":
                delta = rho/rhoc

                def f(T):
                    tau = Tc/T
                    ideal = self._phi0(tau, delta)
                    fio = ideal["fio"]
                    fiot = ideal["fiot"]
                    fir = _phir(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    so = self.R*(tau*(fiot+firt)-fio-fir)
                    return so-s

                T = fsolve(f, To)[0]
                rhol = self._Liquid_Density(T)
                rhov = self._Vapor_Density(T)
                if T == To or rhov <= rho <= rhol:

                    def f(parr):
                        T, rhol, rhog = parr
                        tau = Tc/T
                        deltaL = rhol/self.rhoc
                        deltaG = rhog/self.rhoc

                        idealL = self._phi0(tau, deltaL)
                        fioL = idealL["fio"]
                        fiot = idealL["fiot"]
                        idealG = self._phi0(tau, deltaG)
                        fioG = idealG["fio"]

                        firL = _phir(tau, deltaL, self._constants)
                        firdL = _phird(tau, deltaL, self._constants)
                        firtL = _phirt(tau, deltaL, self._constants)
                        soL = self.R*(tau*(fiot+firtL)-fioL-firL)
                        firG = _phir(tau, deltaG, self._constants)
                        firdG = _phird(tau, deltaG, self._constants)
                        firtG = _phirt(tau, deltaG, self._constants)
                        soG = self.R*(tau*(fiot+firtG)-fioG-firG)

                        Jl = rhol*(1+deltaL*firdL)
                        Jv = rhog*(1+deltaG*firdG)
                        K = firL-firG
                        x = (1./rho-1/rhol)/(1/rhog-1/rhol)
                        return (Jl-Jv,
                                Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
                                soL*(1-x)+soG*x - s)

                    for to in [To, 300, 400, 500, 600]:
                        rhoLo = self._Liquid_Density(to)
                        rhoGo = self._Vapor_Density(to)
                        sol = fsolve(f, [to, rhoLo, rhoGo], full_output=True)
                        T, rhoL, rhoG = sol[0]
                        x = (1./rho-1/rhoL)/(1/rhoG-1/rhoL)
                        if sol[2] == 1 and 0 <= x <= 1 and \
                                sum(abs(sol[1]["fvec"])) < 1e-5:
                            break

                    if sum(abs(sol[1]["fvec"])) > 1e-5:
                        raise(RuntimeError(sol[3]))

                    liquido = self._Helmholtz(rhoL, T)
                    vapor = self._Helmholtz(rhoG, T)
                    P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
                        liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000

            elif self._mode == "rhou":
                delta = rho/rhoc

                def f(T):
                    tau = Tc/T

                    ideal = self._phi0(tau, delta)
                    fiot = ideal["fiot"]
                    fird = _phird(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    Po = (1+delta*fird)*self.R*T*rho
                    ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
                    return ho-Po/rho-u

                T = fsolve(f, To)[0]
                rhol = self._Liquid_Density(T)
                rhov = self._Vapor_Density(T)
                if T == To or rhov <= rho <= rhol:
                    def f(parr):
                        T, rhol, rhog = parr
                        tau = Tc/T
                        deltaL = rhol/self.rhoc
                        deltaG = rhog/self.rhoc

                        ideal = self._phi0(tau, deltaL)
                        fiot = ideal["fiot"]
                        firL = _phir(tau, deltaL, self._constants)
                        firdL = _phird(tau, deltaL, self._constants)
                        firtL = _phirt(tau, deltaL, self._constants)
                        hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
                        firG = _phir(tau, deltaG, self._constants)
                        firdG = _phird(tau, deltaG, self._constants)
                        firtG = _phirt(tau, deltaG, self._constants)
                        hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)

                        Jl = rhol*(1+deltaL*firdL)
                        Jv = rhog*(1+deltaG*firdG)
                        K = firL-firG
                        x = (1./rho-1/rhol)/(1/rhog-1/rhol)
                        Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
                        vu = hoG-Ps/rhog
                        lu = hoL-Ps/rhol
                        return (Jl-Jv,
                                Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
                                lu*(1-x)+vu*x - u)

                    for to in [To, 300, 400, 500, 600]:
                        rhoLo = self._Liquid_Density(to)
                        rhoGo = self._Vapor_Density(to)
                        sol = fsolve(f, [to, rhoLo, rhoGo], full_output=True)
                        T, rhoL, rhoG = sol[0]
                        x = (1./rho-1/rhoL)/(1/rhoG-1/rhoL)
                        if sol[2] == 1 and 0 <= x <= 1 and \
                                sum(abs(sol[1]["fvec"])) < 1e-5:
                            break

                    if sum(abs(sol[1]["fvec"])) > 1e-5:
                        raise(RuntimeError(sol[3]))

                    liquido = self._Helmholtz(rhoL, T)
                    vapor = self._Helmholtz(rhoG, T)
                    P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
                        liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000

            elif self._mode == "hs":
                def f(parr):
                    rho, T = parr
                    delta = rho/rhoc
                    tau = Tc/T

                    ideal = self._phi0(tau, delta)
                    fio = ideal["fio"]
                    fiot = ideal["fiot"]
                    fird = _phird(tau, delta, self._constants)
                    fir = _phir(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
                    so = self.R*(tau*(fiot+firt)-fio-fir)
                    return ho-h, so-s

                rho, T = fsolve(f, [rhoo, To])
                rhol = self._Liquid_Density(T)
                rhov = self._Vapor_Density(T)
                if rhov <= rho <= rhol:

                    def f(parr):
                        T, rhol, rhog, x = parr
                        tau = Tc/T
                        deltaL = rhol/self.rhoc
                        deltaG = rhog/self.rhoc

                        idealL = self._phi0(tau, deltaL)
                        fiot = idealL["fiot"]
                        fioL = idealL["fio"]
                        idealG = self._phi0(tau, deltaG)
                        fioG = idealG["fio"]

                        firL = _phir(tau, deltaL, self._constants)
                        firdL = _phird(tau, deltaL, self._constants)
                        firtL = _phirt(tau, deltaL, self._constants)
                        hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
                        soL = self.R*(tau*(fiot+firtL)-fioL-firL)
                        firG = _phir(tau, deltaG, self._constants)
                        firdG = _phird(tau, deltaG, self._constants)
                        firtG = _phirt(tau, deltaG, self._constants)
                        hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
                        soG = self.R*(tau*(fiot+firtG)-fioG-firG)

                        Jl = rhol*(1+deltaL*firdL)
                        Jv = rhog*(1+deltaG*firdG)
                        K = firL-firG
                        return (Jl-Jv,
                                Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
                                hoL*(1-x)+hoG*x - h,
                                soL*(1-x)+soG*x - s)

                    for to in [To, 300, 400, 500, 600]:
                        rLo = self._Liquid_Density(to)
                        rGo = self._Vapor_Density(to)
                        sol = fsolve(f, [to, rLo, rGo, 0.5], full_output=True)
                        T, rhoL, rhoG, x = sol[0]
                        if sol[2] == 1 and 0 <= x <= 1 and \
                                sum(abs(sol[1]["fvec"])) < 1e-5:
                            break

                    if sum(abs(sol[1]["fvec"])) > 1e-5:
                        raise(RuntimeError(sol[3]))

                    liquido = self._Helmholtz(rhoL, T)
                    vapor = self._Helmholtz(rhoG, T)
                    P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
                        liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000

            elif self._mode == "hu":
                def f(parr):
                    rho, T = parr
                    delta = rho/rhoc
                    tau = Tc/T

                    ideal = self._phi0(tau, delta)
                    fiot = ideal["fiot"]
                    fird = _phird(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    Po = (1+delta*fird)*self.R*T*rho
                    ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)

                    return ho-Po/rho-u, ho-h

                sol = fsolve(f, [rhoo, To], full_output=True)
                rho, T = sol[0]
                rhol = self._Liquid_Density(T)
                rhov = self._Vapor_Density(T)
                if sol[2] != 1 or rhov <= rho <= rhol:

                    def f(parr):
                        T, rhol, rhog, x = parr
                        tau = Tc/T
                        deltaL = rhol/self.rhoc
                        deltaG = rhog/self.rhoc

                        ideal = self._phi0(tau, deltaL)
                        fiot = ideal["fiot"]

                        firL = _phir(tau, deltaL, self._constants)
                        firdL = _phird(tau, deltaL, self._constants)
                        firtL = _phirt(tau, deltaL, self._constants)
                        hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
                        firG = _phir(tau, deltaG, self._constants)
                        firdG = _phird(tau, deltaG, self._constants)
                        firtG = _phirt(tau, deltaG, self._constants)
                        hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)

                        Jl = rhol*(1+deltaL*firdL)
                        Jv = rhog*(1+deltaG*firdG)
                        K = firL-firG

                        Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
                        vu = hoG-Ps/rhog
                        lu = hoL-Ps/rhol

                        return (Jl-Jv,
                                Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
                                hoL*(1-x)+hoG*x - h,
                                lu*(1-x)+vu*x - u)

                    for to in [To, 300, 400, 500, 600]:
                        rLo = self._Liquid_Density(to)
                        rGo = self._Vapor_Density(to)
                        sol = fsolve(f, [to, rLo, rGo, 0.5], full_output=True)
                        T, rhoL, rhoG, x = sol[0]
                        if sol[2] == 1 and 0 <= x <= 1 and \
                                sum(abs(sol[1]["fvec"])) < 1e-5:
                            break

                    if sum(abs(sol[1]["fvec"])) > 1e-5:
                        raise(RuntimeError(sol[3]))

                    liquido = self._Helmholtz(rhoL, T)
                    vapor = self._Helmholtz(rhoG, T)
                    P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
                        liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000

            elif self._mode == "su":
                def f(parr):
                    rho, T = parr
                    delta = rho/rhoc
                    tau = Tc/T

                    ideal = self._phi0(tau, delta)
                    fio = ideal["fio"]
                    fiot = ideal["fiot"]
                    fird = _phird(tau, delta, self._constants)
                    fir = _phir(tau, delta, self._constants)
                    firt = _phirt(tau, delta, self._constants)
                    ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
                    so = self.R*(tau*(fiot+firt)-fio-fir)
                    Po = (1+delta*fird)*self.R*T*rho
                    return ho-Po/rho-u, so-s

                sol = fsolve(f, [rhoo, To], full_output=True)
                rho, T = sol[0]
                rhol = self._Liquid_Density(T)
                rhov = self._Vapor_Density(T)
                if sol[2] != 1 or rhov <= rho <= rhol:

                    def f(parr):
                        T, rhol, rhog, x = parr
                        tau = Tc/T
                        deltaL = rhol/self.rhoc
                        deltaG = rhog/self.rhoc

                        idealL = self._phi0(tau, deltaL)
                        fiot = idealL["fiot"]
                        fioL = idealL["fio"]
                        idealG = self._phi0(tau, deltaG)
                        fioG = idealG["fio"]

                        firL = _phir(tau, deltaL, self._constants)
                        firdL = _phird(tau, deltaL, self._constants)
                        firtL = _phirt(tau, deltaL, self._constants)
                        hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
                        soL = self.R*(tau*(fiot+firtL)-fioL-firL)
                        firG = _phir(tau, deltaG, self._constants)
                        firdG = _phird(tau, deltaG, self._constants)
                        firtG = _phirt(tau, deltaG, self._constants)
                        hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
                        soG = self.R*(tau*(fiot+firtG)-fioG-firG)

                        Jl = rhol*(1+deltaL*firdL)
                        Jv = rhog*(1+deltaG*firdG)
                        K = firL-firG

                        Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
                        vu = hoG-Ps/rhog
                        lu = hoL-Ps/rhol

                        return (Jl-Jv,
                                Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
                                soL*(1-x)+soG*x - s,
                                lu*(1-x)+vu*x - u)

                    for to in [To, 300, 400, 500, 600]:
                        rLo = self._Liquid_Density(to)
                        rGo = self._Vapor_Density(to)
                        sol = fsolve(f, [to, rLo, rGo, 0.5], full_output=True)
                        T, rhoL, rhoG, x = sol[0]
                        if sol[2] == 1 and 0 <= x <= 1 and \
                                sum(abs(sol[1]["fvec"])) < 1e-5:
                            break

                    if sum(abs(sol[1]["fvec"])) > 1e-5:
                        raise(RuntimeError(sol[3]))

                    liquido = self._Helmholtz(rhoL, T)
                    vapor = self._Helmholtz(rhoG, T)
                    P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
                        liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000

            elif self._mode == "Trho":
                if T < self.Tc:
                    rhov = self._Vapor_Density(T)
                    rhol = self._Liquid_Density(T)
                    if rhol > rho > rhov:
                        rhol, rhov, Ps = self._saturation(T)
                        if rhol > rho > rhov:
                            vapor = self._Helmholtz(rhov, T)
                            liquido = self._Helmholtz(rhol, T)
                            x = (1/rho-1/rhol)/(1/rhov-1/rhol)
                            P = Ps/1000

            rho = float(rho)
            T = float(T)
            propiedades = self._Helmholtz(rho, T)

            if T > self.Tc:
                x = 1
            elif x is None:
                x = 0

            if not P:
                P = propiedades["P"]/1000.

        elif self._mode == "Tx":
            # Check input T in saturation range
            if self.Tt > T or self.Tc < T or x > 1 or x < 0:
                raise NotImplementedError("Incoming out of bound")

            rhol, rhov, Ps = self._saturation(T)
            vapor = self._Helmholtz(rhov, T)
            liquido = self._Helmholtz(rhol, T)
            if x == 0:
                propiedades = liquido
            elif x == 1:
                propiedades = vapor
            P = Ps/1000.

        elif self._mode == "Px":
            # Check input P in saturation range
            if self.Pc < P or x > 1 or x < 0:
                raise NotImplementedError("Incoming out of bound")

            # Iterate over saturation routine to get T
            def f(T):
                rhol = self._Liquid_Density(T)
                rhog = self._Vapor_Density(T)

                deltaL = rhol/self.rhoc
                deltaG = rhog/self.rhoc

                tau = Tc/T

                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 Ps/1000-P

            if T0:
                To = T0
            elif self.name == "water":
                To = _TSat_P(P)
            else:
                To = (self.Tc+self.Tt)/2
            T = fsolve(f, To)[0]
            rhol, rhov, Ps = self._saturation(T)
            vapor = self._Helmholtz(rhov, T)
            liquido = self._Helmholtz(rhol, T)
            if x == 0:
                propiedades = liquido
            elif x == 1:
                propiedades = vapor

        self.T = T
        self.Tr = T/self.Tc
        self.P = P
        self.Pr = self.P/self.Pc
        self.x = x
        if self._mode in ["Tx", "Px"] or 0 < x < 1:
            region = 4
        else:
            region = 0
        self.phase = getphase(self.Tc, self.Pc, self.T, self.P, self.x, region)

        self.Liquid = _fase()
        self.Gas = _fase()
        if x == 0:
            # liquid phase
            self.fill(self.Liquid, propiedades)
            self.fill(self, propiedades)
        elif x == 1:
            # vapor phase
            self.fill(self.Gas, propiedades)
            self.fill(self, propiedades)
        else:
            self.fill(self.Liquid, liquido)
            self.fill(self.Gas, vapor)

            self.v = x*self.Gas.v+(1-x)*self.Liquid.v
            self.rho = 1./self.v

            self.h = x*self.Gas.h+(1-x)*self.Liquid.h
            self.s = x*self.Gas.s+(1-x)*self.Liquid.s
            self.u = x*self.Gas.u+(1-x)*self.Liquid.u
            self.a = x*self.Gas.a+(1-x)*self.Liquid.a
            self.g = x*self.Gas.g+(1-x)*self.Liquid.g

            self.Z = x*self.Gas.Z+(1-x)*self.Liquid.Z
            self.f = x*self.Gas.f+(1-x)*self.Liquid.f

            self.Z_rho = x*self.Gas.Z_rho+(1-x)*self.Liquid.Z_rho
            self.IntP = x*self.Gas.IntP+(1-x)*self.Liquid.IntP

        # Calculate special properties useful only for one phase
        if self._mode in ("Px", "Tx") or (x < 1 and self.Tt <= T <= self.Tc):
            self.sigma = self._surface(T)
        else:
            self.sigma = None

        vir = self._virial(T)
        self.virialB = vir["B"]/self.rhoc
        self.virialC = vir["C"]/self.rhoc**2

        if 0 < x < 1:
            self.Hvap = vapor["h"]-liquido["h"]
            self.Svap = vapor["s"]-liquido["s"]
        else:
            self.Hvap = None
            self.Svap = None

        self.invT = -1/self.T

        # Ideal properties
        cp0 = self._prop0(self.rho, self.T)
        self.v0 = self.R*self.T/self.P/1000
        self.rho0 = 1./self.v0
        self.h0 = cp0.h
        self.u0 = self.h0-self.P*self.v0
        self.s0 = cp0.s
        self.a0 = self.u0-self.T*self.s0
        self.g0 = self.h0-self.T*self.s0
        self.cp0 = cp0.cp
        self.cv0 = cp0.cv
        self.cp0_cv = self.cp0/self.cv0
        cp0.v = self.v0
        self.gamma0 = -self.v0/self.P/1000*self.derivative("P", "v", "s", cp0)