z.hallyarboroughL <- function(pres.a, temp.f, gas.sg, pres.pc, temp.pc,
                              pres.pr, temp.pr,
                              n2.frac = 0, co2.frac = 0, h2s.frac = 0,
                              verbose = FALSE, ...)
{
    funcY <- function(y) {
        # implicit equation
        # in some literature A = A1, B = A2, C = A3, D = A4

        - A * pres.pr + (y + y^2 + y^3 - y^4) / (1 - y)^3  - B * y^2 + C * y^D
    }

    if (missing(gas.sg)) {
        if (verbose) cat("No gas.sg supplied. ")
        if (missing(pres.pc) || missing(temp.pc)) {
            if (missing(pres.pr) || missing(temp.pr)) stop()
            if (missing(pres.a) || missing(temp.f)) {
                #: user ONLY supplies pseudo-reduced P, T
                if (verbose) cat("Using ONLY pres.pr and temp.pr to calculate z \n")
                temp.r  <- 1 / temp.pr
                pres.pc <- NA
                temp.pc <- NA
            } else {
                #: providing pres.a, temp.f, pres.pr, temp.pr
                if (verbose) cat("Using instead pres.pr and temp.pr \n")
                temp.r <- 1 / temp.pr
                # calculate pseudo-criticals
                pres.pc <- pres.a / pres.pr
                temp.pc <- (temp.f + 460) / temp.pr
                # Guo's worksheet has bug in the Farenheit add
            }
        } else {
            if (missing(pres.pc) || missing(temp.pc)) stop()
            if (verbose)cat("Using instead Ppc and Tpc \n")
            crit <- calcGasPseudoReduced(pres.a, pres.pc, temp.f, temp.pc)
            pres.pr <- crit$pres.pr
            temp.pr <- crit$temp.pr
            temp.r  <- crit$temp.r
        }
    } else {
        #: the standard calculation when specific gravity of the gas is provided
        if (verbose) cat("gas.sg has been provided. Will calculate Ppc, Tpc, Ppr, Tpr \n")
        crit <- calcCriticals(pres.a, temp.f, gas.sg,
                              co2.frac = 0, h2s.frac = 0, n2.frac = 0)
        pres.pr <- crit$pres.pr
        temp.pr <- crit$temp.pr
        temp.r  <- crit$temp.r
        pres.pc <- crit$pres.pc
        temp.pc <- crit$temp.pc
    }
    #: calculate reduced temperature
    t <- temp.r   # make it easier to read in the equation below
    print(temp.pr)
    print(temp.r)
    #: calculate constants given pseudo-reduced temperature
    #: formulas have been checked against five sources: papers and books
    A <- 0.06125 * t * exp(-1.2 * (1 - t)^2)
    B <- t * (14.76 - 9.76 * t + 4.58 * t^2)
    C <- t * (90.7 - 242.2 * t + 42.4 * t^2)
    D <- 2.18 + 2.82 * t

    All <- rootSolve::uniroot.all(funcY, c(-100.01, 100.99)) # find the root
    Y <- min(All)                         # take the minimum value
    z <- A * pres.pr / Y
    #     #: prepare for table
    # zfactors <- named.list(z, Y, A, B, C, D,
    #                        pres.pr, temp.pr, pres.pc, temp.pc, temp.r)
    # return(zfactors)
    return(z)
}
z.hallyarboroughL(pres.pr = 0.5, temp.pr = 1.3)


f0nzie/old.zFactor documentation built on May 24, 2019, 4:03 a.m.