Python scipy.log10() Examples

def _visco0(self, rho, T, fase=None):
        a = [17.67484, -2.78751, 311498.7, -48826500, 3938774000, -1.654629e11,
             2.86561e12]
        Tr = T/0.29944
        y = 0.68321*(a[0] + a[1]*log10(Tr) + a[2]/Tr**2 + a[3]/Tr**3 +
                     a[4]/Tr**4 + a[5]/Tr**5 + a[6]/Tr**6)
        nt = 266.93*(T*self.M)**0.5/y
        om = rho/1673.0
        c = [1.03010, -0.99175, 2.47127, -3.11864, 1.57066]
        b = [0.48148, -1.18732, 2.80277, -5.41058, 7.04779, -3.76608]
        sum1 = sum([ci*om**i for i, ci in enumerate(c)])
        sum2 = sum([bi*om**i for i, bi in enumerate(b)])
        sigma = 3.05e-10*(sum1-sum2*log10(T/122.1))
        br = 2.0/3.0*pi*Avogadro*sigma**3
        brho = rho/self.M*1000*br
        d = [1, 0.27676, 0.014355, 2.6480, -1.9643, 0.89161]
        nd = sum([di*brho**i for i, di in enumerate(d)])
        return unidades.Viscosity(nd*nt/100, "muPas") 

def addqqplotinfo(qnull,M,xl='-log10(P) observed',yl='-log10(P) expected',xlim=None,ylim=None,alphalevel=0.05,legendlist=None,fixaxes=False):    
    distr='log10'
    pl.plot([0,qnull.max()], [0,qnull.max()],'k')
    pl.ylabel(xl)
    pl.xlabel(yl)
    if xlim is not None:
        pl.xlim(xlim)
    if ylim is not None:
        pl.ylim(ylim)        
    if alphalevel is not None:
        if distr == 'log10':
            betaUp, betaDown, theoreticalPvals = _qqplot_bar(M=M,alphalevel=alphalevel,distr=distr)
            lower = -sp.log10(theoreticalPvals-betaDown)
            upper = -sp.log10(theoreticalPvals+betaUp)
            pl.fill_between(-sp.log10(theoreticalPvals),lower,upper,color="grey",alpha=0.5)
            #pl.plot(-sp.log10(theoreticalPvals),lower,'g-.')
            #pl.plot(-sp.log10(theoreticalPvals),upper,'g-.')
    if legendlist is not None:
        leg = pl.legend(legendlist, loc=4, numpoints=1)
        # set the markersize for the legend
        for lo in leg.legendHandles:
            lo.set_markersize(10)

    if fixaxes:
        fix_axes() 

def Viscosidad_ASTM(self, T, T1=unidades.Temperature(100, "F"), T2=unidades.Temperature(210, "F"), mu1=0, mu2=0):
        """Cálculo de la viscosidad cinemática a cualquier temperatura, conociendo otros dos valores de viscosidad a otras temperaturas, API procedure 11A4.4, pag 1063
        Parámetros:
        T:Temperatura a la que se quiere calcular la viscosidad
        T1,T2:opcional, temperatura a la que se conoce la viscosidad
        mu1,mu2:opcionales, valores de viscosidad conocidos
        Si no se suministran los parámetros opcionales se consideran los valores a 100 y 210ºF
        """
        if mu1==0:
            mu1=self.v100.cSt
        if mu2==0:
            mu2=self.v210.cSt
        t=unidades.Temperature(T)
        Z1=mu1+0.7+exp(-1.47-1.84*mu1-0.51*mu1**2)
        Z2=mu2+0.7+exp(-1.47-1.84*mu2-0.51*mu2**2)
        B=(log10(log10(Z1))-log10(log10(Z2)))/(log10(T1.R)-log10(T2.R))
        Z=10**(10**(log10(log10(Z1))+B*(log10(t.R)-log10(T1.R))))
        return unidades.Diffusivity(Z-0.7-exp(-0.7487-3.295*(Z-0.7)+0.6119*(Z-0.7)**2-0.3191*(Z-0.7)**3), "cSt") 

def w_Korsten(Tb, Tc, Pc):
    """Calculate petroleum fractions acentric factor with the Korsten (2000)
    correlation

    Parameters
    ------------
    Tb : float
        Normal boiling temperature, [K]
    Tc : float
        Critical temperature, [K]
    Pc : float
        Critical pressure, [Pa]

    Returns
    -------
    w: float
        Acentric factor, [-]
    """
    tbr = Tb/Tc
    # Eq 29
    w = 0.5899*tbr**1.3/(1-tbr**1.3)*log10(Pc/101325)-1.
    return {"w": unidades.Dimensionless(w)} 

def tdft(audio, srate, windowsize, windowshift,fftsize):

    """Calculate the real valued fast Fourier transform of a segment of audio multiplied by a 
    a Hamming window.  Then, convert to decibels by multiplying by 20*log10.  Repeat for all
    segments of the audio."""
    
    windowsamp = int(windowsize*srate)
    shift = int(windowshift*srate)
    window = scipy.hamming(windowsamp)
    spectrogram = scipy.array([20*scipy.log10(abs(np.fft.rfft(window*audio[i:i+windowsamp],fftsize))) 
                     for i in range(0, len(audio)-windowsamp, shift)])
    return spectrogram 

def get_base(unit='bit'):
    if unit == 'bit':
        log = sp.log2
    elif unit == 'nat':
        log = sp.log
    elif unit in ('digit', 'dit'):
        log = sp.log10
    else:
        raise ValueError('The "unit" "%s" not understood' % unit)
    return log 

def _qqplot_bar(M=1000000, alphalevel = 0.05,distr = 'log10'):
    '''
    calculate error bars for a QQ-plot
    --------------------------------------------------------------------
    Input:
    -------------   ----------------------------------------------------
    M               number of points to compute error bars
    alphalevel      significance level for the error bars (default 0.05)
    distr           space in which the error bars are implemented
                    Note only log10 is implemented (default 'log10')
    --------------------------------------------------------------------
    Returns:
    -------------   ----------------------------------------------------
    betaUp          upper error bars
    betaDown        lower error bars
    theoreticalPvals    theoretical P-values under uniform
    --------------------------------------------------------------------
    '''


    #assumes 'log10'

    mRange=10**(sp.arange(sp.log10(0.5),sp.log10(M-0.5)+0.1,0.1));#should be exp or 10**?
    numPts=len(mRange);
    betaalphaLevel=sp.zeros(numPts);#down in the plot
    betaOneMinusalphaLevel=sp.zeros(numPts);#up in the plot
    betaInvHalf=sp.zeros(numPts);
    for n in xrange(numPts):
        m=mRange[n]; #numplessThanThresh=m;
        betaInvHalf[n]=st.beta.ppf(0.5,m,M-m);
        betaalphaLevel[n]=st.beta.ppf(alphalevel,m,M-m);
        betaOneMinusalphaLevel[n]=st.beta.ppf(1-alphalevel,m,M-m);
        pass
    betaDown=betaInvHalf-betaalphaLevel;
    betaUp=betaOneMinusalphaLevel-betaInvHalf;

    theoreticalPvals=mRange/M;
    return betaUp, betaDown, theoreticalPvals 

def calculate(config):
    fanning = config.getboolean("Moody", "fanning")
    method = config.getint("Moody", "method")
    eD = eval(config.get("Moody", "eD"))
    F = f_list[method]

    dat = {}
    dat["fanning"] = fanning
    dat["method"] = method
    # laminar
    if fanning:
        dat["laminar"] = [16./R for R in Re_laminar]
        x = 4
    else:
        dat["laminar"] = [64./R for R in Re_laminar]
        x = 1

    # turbulent
    turb = {}
    for e in eD:
        turb[e] = [F(Rei, e)/x for Rei in Re_turbulent]
        dat["turbulent"] = turb

    # Line to define the fully desarrolled turbulent flux
    dat["fully"] = [(1/(1.14-2*log10(3500/R)))**2/x for R in Re_fully]

    # Save to file
    with open(conf_dir+"moody.dat", "w") as file:
        json.dump(dat, file, indent=4) 

def Viscosity_Index(self):
        """Método de calculo del indice de viscosidad, API procedure 11A6.1, pag 1083"""
        if self.v210.cSt>70:
            L=0.8353*self.v210.cSt**2+14.67*self.v210.cSt-216
            H=0.1684*self.v210.cSt**2+11.85*self.v210.cSt-97
        else: #Ajuste de los datos de la tabla, producción propia
            """Polynomial Fit of dataset: Table1_L, using function: a0+a1*x+a2*x^2+a3*x^3+a4*x^4+a5*x^5
            --------------------------------------------------------------------------------------
            Chi^2/doF = 1,4854228443647e+00
            R^2 = 0,9999990418201
            Adjusted R^2 = 0,9999990229086
            RMSE (Root Mean Squared Error) = 1,218779243491
            RSS (Residual Sum of Squares) = 453,0539675312
            ---------------------------------------------------------------------------------------

            Polynomial Fit of dataset: Table1_H, using function: a0+a1*x+a2*x^2+a3*x^3+a4*x^4+a5*x^5
            --------------------------------------------------------------------------------------
            Chi^2/doF = 1,7949865589785e-01
            R^2 = 0,999998895674
            Adjusted R^2 = 0,9999988738781
            RMSE (Root Mean Squared Error) = 0,4236728170391
            RSS (Residual Sum of Squares) = 54,74709004884
            ---------------------------------------------------------------------------------------
            """
            L=-1.2756691045812e+01+6.1466654146190*self.v210.cSt+9.9774520581931e-01*self.v210.cSt**2-2.5045430263656e-03*self.v210.cSt**3+3.1748553181177e-05*self.v210.cSt**4-1.8264076604682e-07*self.v210.cSt**5
            H=-7.6607640162508+5.4845434135144*self.v210.cSt+0.38222987985934*self.v210.cSt**2-4.6556329076069e-03*self.v210.cSt**3+5.1653200038471e-05*self.v210.cSt**4-2.3274903246922e-07*self.v210.cSt**5
        if self.v100.cSt>H:
            VI=(L-self.v100.cSt)/(L-H)*100
        else:
            N=(log10(H)-log10(self.v100.cSt))/log10(self.v210.cSt)
            VI=(10**N-1)/0.00715+100
        return VI 

def Mu_Singh(self, T, mu):
        """Calculo de la viscosidad cinemática de fracciones petrolíferas a baja presión, conocida la viscosidad cinemática a 100ºF, API procedure 11A4.1 pag 1054
        mu: viscosidad cinematica experimental a 100ºF
        valor obtenido en centistokes"""
        t=unidades.Temperature(T)
        S=0.28008*log10(mu)+1.8616
        B=log10(mu)+0.86960
        return 10**(B*(559.67/t.R)**S-0.86960) 

def Pv_Maxwell_Bonnell(self, T, mod=False):
        """Maxwell, J. B. and Bonnell, L. S., Vapor Pressure Charts for Petroleum Engineers, Exxon Research and Engineering Company, Florham Park, NJ, 1955. Reprinted in 1974."Deviation and Precision of a New Vapor Pressure Correlation for Petroleum Hydrocarbons," Industrial and Engineering Chemistry, Vol. 49, 1957, pp. 1187-1196.
        modificación de Tsonopoulos para coal liquids:  Tsonopoulos, C., Heidman, J. L., and Hwang, S.-C.,Thermodynamic and Transport Properties of Coal Liquids, An Exxon Monograph, Wiley, New York, 1986."""
        if mod:
            if self.Tb<=366.5:
                F1=0
            else:
                F1=-1+0.009*(self.Tb-255.37)
            F2=(self.watson-12)-0.01304*(self.watson-12)**2
        else:
            if  self.Tb<367:
                F=0
            elif self.Tb<478:
                F=-3.2985+0.0009*self.Tb
            else:
                F=-3.2985+0.009*self.Tb

        pvap=0.
        pvapcalc=1.
        while abs(pvap-pvapcalc)<1e-6:
            pvap=pvapcalc
            if mod:
                if pvap<=760:
                    F3=1.47422*log10(pvap/760.)
                else:
                    F3=1.47422*log10(pvap/760.)+1.190833*(log10(pvap/760.))**2
                DTb=F1*F2*F3
            else:
                DTb=1.3889*F*(self.watson-12)*log10(pvap/760.)
            Tb_=self.Tb-DTb
            Q=(Tb_/T-0.00051606*Tb_)/(748.1-0.3861*Tb_)
            if Q>0.0022:
                pvapcalc=10**((3000.538*Q-6.76156)/(43*Q-0.987672))
            elif Q>=0.0013:
                pvapcalc=10**((2663.129*Q-5.994296)/(95.76*Q-0.972546))
            else:
                pvapcalc=10**((2770.085*Q-6.412631)/(36*Q-0.989679))

        return unidades.Pressure(pvapcalc, "mmHg") 

def MuL_Singh(T, v100):
    """Calculate kinematic viscosity of liquid petroleum fractions at low
    pressure by Singh correlation, also referenced in API Procedure 11A4.1,
    pag 1031

    Parameters
    ------------
    T : float
        Temperature, [K]
    v100 : float
        Kinematic viscosity at 100ºF, [cSt]

    Returns
    -------
    v : float
        Kinematic viscosity at T, [cSt]

    Examples
    --------
    Example from [20]_, Sumatran crude at 210ºF

    >>> T = unidades.Temperature(210, "F")
    >>> "%0.3f" % MuL_Singh(T, 1.38).cSt
    '0.705'
    """
    # Convert input temperature to Rankine
    t_R = unidades.K2R(T)
    S = 0.28008*log10(v100) + 1.8616
    B = log10(v100) + 0.86960
    mu = 10**(B*(559.67/t_R)**S - 0.86960)
    return unidades.Diffusivity(mu, "cSt")


# Component predition 

def prop_Tsonopoulos(SG, Tb):
    """Calculate petroleum fractions properties with the Tsonopoulos (1990)
    correlations using molecular weight and specific gravity as input parameter

    Parameters
    ------------
    SG: float
        Specific gravity, [-]
    Tb: float
        Boiling temperature, [K]

    Returns
    -------
    prop : dict
        A dict with the calculated properties:

            * Tc: Critic temperature, [K]
            * Pc: Critic pressure, [MPa]
    """
    # TODO: Search reference
    Tc = 10**(1.20016 - 0.61954*log10(Tb) + 0.48262*log10(SG) +
              0.67365*log10(SG)**2)
    Pc = 10**(7.37498 - 2.15833*log10(Tb) + 3.35417*log10(SG) +
              5.64019*log10(SG)**2)

    prop = {}
    prop["Tb"] = unidades.Temperature(Tb)
    prop["SG"] = unidades.Dimensionless(SG)
    prop["Tc"] = unidades.Temperature(Tc)
    prop["Pc"] = unidades.Pressure(Pc, "bar")
    return prop 

def prop_Standing(M, SG):
    """Calculate petroleum fractions properties with the Standing (1977)
    correlations based in the graphical plot of Mathews et al. (1942) using
    molecular weight and specific gravity as input parameter

    Parameters
    ------------
    M: float
        Molecular weight, [-]
    SG: float
        Specific gravity, [-]

    Returns
    -------
    prop : dict
        A dict with the calculated properties:

            * Tc: Critic temperature, [K]
            * Pc: Critic pressure, [atm]

    Examples
    --------
    Example 2.3 from [9]_: Petroleum fraction with M=216 and SG=0.8605

    >>> p = prop_Standing(216, 0.8605)
    >>> "%.1f %.0f" % (p["Tc"].R, p["Pc"].psi)
    '1269.3 270'
    """
    Tc = 608 + 364*log10(M-71.2) + (2450*log10(M)-3800)*log10(SG)       # Eq 25
    Pc = 1188 - 431*log10(M-61.1) + (2319-852*log10(M-53.7))*(SG-0.8)   # Eq 26

    prop = {}
    prop["M"] = unidades.Dimensionless(M)
    prop["SG"] = unidades.Dimensionless(SG)
    prop["Tc"] = unidades.Temperature(Tc, "R")
    prop["Pc"] = unidades.Pressure(Pc, "psi")
    return prop 

def _Pv(T, Tcm, Pcm, wm):
    """Pseudo vapor presure from pseudocritical properties to use in
    compressed liquid density correlations

    Parameters
    ----------
    T : float
        Temperature, [K]
    Tcm : float
        Pseudo critical temperature of mixture, [K]
    Pcm : float
        Pseudo critical pressure of mixture, [Pa]
    wm : float
        Acentric factor of mixture, [-]

    Returns
    -------
    Pbp : float
        Boiling point pressure, [Pa]
    """
    Trm = T/Tcm
    alpha = 35 - 36/Trm - 96.736*log10(Trm) + Trm**6                   # Eq 22
    beta = log10(Trm) + 0.03721754*alpha                               # Eq 23
    Pro = 5.8031817*log10(Trm)+0.07608141*alpha                        # Eq 20

    # Using Aalto modificated equation Eq 8
    Pr1 = 4.86601*beta+0.03721754*alpha                                # Eq 21

    Pbp = Pcm*10**(Pro+wm*Pr1)
    return unidades.Pressure(Pbp)


# Liquid viscosity correlations 

def _nio(self, T, P):
        """Liquid fugacity coefficient"""

        nio = []
        for cmp in self.componente:
            Tr = T/cmp.Tc
            Pr = P/cmp.Pc

            # Modified Parameters from [2]_
            if cmp.id == 1:  # Hydrogen
                A = [1.50709, 2.74283, -.0211, .00011, 0, .008585, 0, 0, 0, 0]
            elif cmp.id == 2:  # Methane
                A = [1.36822, -1.54831, 0, 0.02889, -0.01076, 0.10486,
                     -0.02529, 0, 0, 0]
            else:
                A = [2.05135, -2.10899, 0, -0.19396, 0.02282, 0.08852, 0,
                     -0.00872, -0.00353, 0.00203]

            # Eq 3
            logn0 = A[0] + A[1]/Tr + A[2]*Tr + A[3]*Tr**2 + A[4]*Tr**3 + \
                (A[5]+A[6]*Tr+A[7]*Tr**2)*Pr + (A[8]+A[9]*Tr)*Pr**2 - log10(Pr)

            # Eq 4
            logn1 = -4.23893 + 8.65808*Tr - 1.2206/Tr - 3.15224*Tr**3 - \
                0.025*(Pr-0.6)

            # Eq 2
            nio.append(10**(logn0 + cmp.f_acent_mod*logn1))
            # The correlation use the modified acentric factor table

        return nio 

def h_anulli_Turbulent(Re, Pr, a, dhL=0, boundary=0):
    """VDI Heat Atlas G2 Pag.703"""
    if boundary == 0:         #Inner surface heated
        Fann = 0.75*a**-0.17
    elif boundary == 1:       #Outer surface heated
        Fann = (0.9-0.15*a**0.6)
    elif boundary == 2:       #Both surfaces heated
        Fann = (0.75*a**-0.17+(0.9-0.15*a**0.6))/(1+a)

    Re_ = Re*((1+a**2)*log(a)+(1-a**2))/((1-a)**2*log(a))
    Xann = (1.8*log10(Re_)-1.5)**-2
    k1 = 1.07+900/Re-0.63/(1+10*Pr)
    Nu = Xann/8*Re*Pr/(k1+12.7*(Xann/8)**0.5*(Pr**(2./3.)-1))*(1+dhL**(2./3.))*Fann
    return Nu 

def Viscosity_Carr_Kobayashi_Burrows(self, T):
        """Carr, N., R. Kobayashi, and D. Burrows. “Viscosity of Hydrocarbon Gases under Pressure.” Transactions of the AIME 201 (1954): 270–275."""
        muo=8.118e-3-6.15e-3*log10(self.SG)+(1.709e-5-2.062e-6*self.SG)*unidades.Temperature(T, "R").F
        muN2=self.N2*(8.49e-3*log10(self.SG)+9.59e-3)
        muCO2=self.CO2*(9.08e-3*log10(self.SG)+6.24e-3)
        muH2S=self.H2S*(8.49e-3*log10(self.SG)+3.73e-3)
        return unidades.Viscosity(muo+muN2+muCO2+muH2S, "cP") 

def _Sublimation_Pressure(cls, T):
        """Special sublimation pressure correlation"""
        # Use decimal logarithm
        P = 10**(-43.39/T+2.5*log10(T)+2.047)
        return unidades.Pressure(P, "mmHg") 

def Mu_Kartoatmodjo_Schmidt(self, T):
        """Kartoatmodjo, T. and Schmidt, Z.: "Large Data Bank Improve Crude Physical Property Correlations," Oil and Gas J. (July 4, 1994) 51-55"""
        t=unidades.Temperature(T)
        mu=16e8*t.F**-2.8177*log10(self.API)**(5.7526*log10(t.F)-26.9718)
        return unidades.Viscosity(mu, "cP") 

def Mu_Glaso(self, T):
        """Glaso, O.: "Generalized Pressure-Volume-Temperature Correlations," J. Pet. Tech. (May 1980), 785-795"""
        t=unidades.Temperature(T)
        mu=3.141e10*t.F**-3.444*log10(self.API)**(10.313*log10(t.F)-36.447)
        return unidades.Viscosity(mu, "cP") 

def Z_Brill_Beggs(Tr, Pr):
    """Calculate gas compressibility factor using the correlation of Brill-
    Beggs (1974)

    Parameters
    ------------
    Tr : float
        Reduced temperature [-]
    Pr : float
        Reduced pressure [-]

    Returns
    -------
    Z : float
        Gas compressibility factor [-]

    Notes
    -----
    Raise :class:`NotImplementedError` if input pair isn't in limit:

        * 1.15 ≤ Tr ≤ 2.4
        * 0.2 ≤ Pr ≤ 15
    """
    # Check input in range of validity
    if Tr < 1.15 or Tr > 2.4 or Pr < 0.2 or Pr > 15:
        raise NotImplementedError("Incoming out of bound")

    A = 1.39*(Tr-0.92)**0.5 - 0.36*Tr - 0.101
    B = (0.62-0.23*Tr)*Pr + (0.066/(Tr-0.86)-0.037)*Pr**2 + \
        0.32/10**(9*(Tr-1))*Pr**6
    C = 0.132 - 0.32*log10(Tr)
    D = 10**(0.3016-0.49*Tr+0.1824*Tr**2)
    Z = A + (1-A)/exp(B) + C*Pr**D
    return unidades.Dimensionless(Z) 

def pb_Vazquez_Beggs(self, T, ts=350, ps=10):
        """Vázquez, M.E. and Beggs, H.D.: "Correlations for Fluid Physical Property Prediction," J.Pet. Tech. (June 1980), 968-970"""
        t=unidades.Temperature(T)
        ts=unidades.Temperature(ts)
        ps=unidades.Pressure(ps, "atm")
        if self.API<=30:
            C1, C2, C3=0.0362, 1.0937, 25.724
        else:
            C1, C2, C3=0.0178, 1.187, 23.931

        gravity_corr=self.gas.SG*(1.+5.912e-5*self.API*ts.F*log10(ps.psi/114.7))
        return unidades.Pressure((self.Rgo.ft3bbl/C1/gravity_corr/exp(C3*self.API/T.R))**(1./C2), "psi") 

def B_Kartoatmodjo_Schmidt(self, T, ts=350, ps=10):
        """Kartoatmodjo, T. and Schmidt, Z.: "Large Data Bank Improve Crude Physical Property Correlations," Oil and Gas J. (July 4, 1994) 51-55"""
        t=unidades.Temperature(T)
        ts=unidades.Temperature(ts)
        ps=unidades.Pressure(ps, "atm")

        gravity_corr=self.gas.SG*(1.+0.1595*self.API**0.4078*ts.F**-0.2506*log10(ps.psi/114.7))
        F=self.Rgo.ft3bbl**0.755*gravity_corr**0.25*self.SG**-1.5+0.45*t.F
        return 0.98496+1e-4*F**1.5 

def B_Glaso(self, T):
        """Glaso, O.: "Generalized Pressure-Volume-Temperature Correlations," J. Pet. Tech. (May 1980), 785-795"""
        t=unidades.Temperature(T)
        F=self.Rgo.ft3bbl(self.gas.SG/self.SG)**0.526*+0.968*t.F
        return 1.+10**(-6.58511+2.91329*log10(F)-0.27683*log10(F)**2) 

def B_Vazquez_Beggs(self, T, ts=350, ps=10):
        """Vázquez, M.E. and Beggs, H.D.: "Correlations for Fluid Physical Property Prediction," J.Pet. Tech. (June 1980), 968-970"""
        t=unidades.Temperature(T)
        ts=unidades.Temperature(ts)
        ps=unidades.Pressure(ps, "atm")
        if self.API<=30:
            C1, C2, C3=4.667e-4, 1.751e-5, -1.8106e-6
        else:
            C1, C2, C3=4.67e-4, 1.1e-5, 1.337e-9

        gravity_corr=self.gas.SG*(1.+5.912e-5*self.API*ts.F*log10(ps.psi/114.7))
        return 1.+C1*self.Rgo.ft3bbl+C2*(t.F-60)*self.API/gravity_corr+C3*self.Rgo.ft3bbl*(t.F-60)*self.API/gravity_corr 

def pb_Kartoatmodjo_Schmidt(self, T, ts=350, ps=10):
        """Kartoatmodjo, T. and Schmidt, Z.: "Large Data Bank Improve Crude Physical Property Correlations," Oil and Gas J. (July 4, 1994) 51-55"""
        t=unidades.Temperature(T)
        ts=unidades.Temperature(ts)
        ps=unidades.Pressure(ps, "atm")
        if self.API<=30:
            C1, C2, C3, C4=0.05958, 0.7972, 13.1405, 0.9986
        else:
            C1, C2, C3, C4=0.0315, 0.7587, 11.2895, 0.9143

        gravity_corr=self.gas.SG*(1.+0.1595*self.API**0.4078*ts.F**-0.2506*log10(ps.psi/114.7))
        return unidades.Pressure((self.Rgo.ft3bbl/C1/gravity_corr**C2/10**(C3*self.API/t.R))**C4, "psi") 

def nichols_plot(sys_list, omega=None, grid=None):
    """Nichols plot for a system

    Plots a Nichols plot for the system over a (optional) frequency range.

    Parameters
    ----------
    sys_list : list of LTI, or LTI
        List of linear input/output systems (single system is OK)
    omega : array_like
        Range of frequencies (list or bounds) in rad/sec
    grid : boolean, optional
        True if the plot should include a Nichols-chart grid. Default is True.

    Returns
    -------
    None
    """
    # Get parameter values
    grid = config._get_param('nichols', 'grid', grid, True)


    # If argument was a singleton, turn it into a list
    if not getattr(sys_list, '__iter__', False):
        sys_list = (sys_list,)

    # Select a default range if none is provided
    if omega is None:
        omega = default_frequency_range(sys_list)

    for sys in sys_list:
        # Get the magnitude and phase of the system
        mag_tmp, phase_tmp, omega = sys.freqresp(omega)
        mag = np.squeeze(mag_tmp)
        phase = np.squeeze(phase_tmp)

        # Convert to Nichols-plot format (phase in degrees,
        # and magnitude in dB)
        x = unwrap(sp.degrees(phase), 360)
        y = 20*sp.log10(mag)

        # Generate the plot
        plt.plot(x, y)

    plt.xlabel('Phase (deg)')
    plt.ylabel('Magnitude (dB)')
    plt.title('Nichols Plot')

    # Mark the -180 point
    plt.plot([-180], [0], 'r+')

    # Add grid
    if grid:
        nichols_grid() 

def CetaneIndex(API, Tb):
    """Calculate cetane index of petroleum fractions using the procedure 2B13
    from API technical databook, pag. 245
    Cetane index is the number equal to the porcentage of cetane in a blend of
    cetane and α-methyl naphthalene having the same ignition quality as a
    sample of the petroleum fraction.

    Parameters
    ------------
    API : float
        API gravity, [-]
    Tb : float
        Mean average boiling point, [ºF]

    Returns
    -------
    CI : float
        Cetane Index, [-]

    Notes
    ------
    Raise :class:`NotImplementedError` if input isn't in limit:

        * 27 ≤ API ≤ 47
        * 360ºF ≤ Tb ≤ 700ºF

    Examples
    --------
    >>> T = unidades.Temperature(617, "F")
    >>> "%0.1f" % CetaneIndex(32.3, T)
    '57.1'
    """
    # Convert input Tb in Kelvin to Fahrenheit to use in the correlation
    Tb_F = unidades.K2F(Tb)

    # Check input in range of validity
    if API < 27 or API > 47 or Tb_F < 360 or Tb_F > 700:
        raise NotImplementedError("Incoming out of bound")

    # Eq 2B13.1-1
    CI = 415.26 - 7.673*API + 0.186*Tb_F + 3.503*API*log10(Tb_F) - \
        193.816*log10(Tb.F)

    return unidades.Dimensionless(CI)


# PNA composition procedures 

def CloudPoint(SG, Tb):
    """Calculate cloud point of petroleum fractions using the procedure 2B12
    from API technical databook, pag. 243
    The cloud point of is the temperature at which its solid paraffin content
    begins to solidify and separate in tiny crystals, causing the oil to appear
    cloudy.

    Parameters
    ------------
    SG : float
        Specific gravity, [-]
    Tb : float
        Mean average boiling point, [ºR]

    Returns
    -------
    CP : float
        Cloud point, [ºR]

    Notes
    ------
    Raise :class:`NotImplementedError` if input isn't in limit:

        * 0.77 ≤ SG ≤ 0.93
        * 800ºR ≤ Tb ≤ 1225ºR

    Examples
    --------
    >>> T = unidades.Temperature(811.5, "R")
    >>> "%0.1f" % CloudPoint(0.787, T).R
    '383.4'
    """
    # Convert input Tb in Kelvin to Rankine to use in the correlation
    Tb_R = unidades.K2R(Tb)

    # Check input in range of validity
    if SG < 0.77 or SG > 0.93 or Tb_R < 800 or Tb_R > 1225:
        raise NotImplementedError("Incoming out of bound")

    # Eq 2B12.1-1
    CP = 10**(-7.41+5.49*log10(Tb_R)-0.712*Tb_R**0.315-0.133*SG)

    return unidades.Temperature(CP, "R")