Manningpara.R

#' Parabolic cross-section for the Gauckler-Manning-Strickler equation
#'
#' This function solves for one missing variable in the Gauckler-Manning-
#' Strickler equation for a parabolic cross-section and uniform flow. The
#' \code{\link[stats]{uniroot}} function is used to obtain the missing parameters.
#'
#'
#'
#'
#' Gauckler-Manning-Strickler equation is expressed as
#'
#' \deqn{V = \frac{K_n}{n}R^\frac{2}{3}S^\frac{1}{2}}
#'
#' \describe{
#'   \item{\emph{V}}{the velocity (m/s or ft/s)}
#'   \item{\emph{n}}{Manning's roughness coefficient (dimensionless)}
#'   \item{\emph{R}}{the hydraulic radius (m or ft)}
#'   \item{\emph{S}}{the slope of the channel bed (m/m or ft/ft)}
#'   \item{\emph{K_n}}{the conversion constant -- 1.0 for SI and
#'        3.2808399 ^ (1 / 3) for English units -- m^(1/3)/s or ft^(1/3)/s}
#' }
#'
#'
#'
#'
#' This equation is also expressed as
#'
#' \deqn{Q = \frac{K_n}{n}\frac{A^\frac{5}{3}}{P^\frac{2}{3}}S^\frac{1}{2}}
#'
#' \describe{
#'   \item{\emph{Q}}{the discharge [m^3/s or ft^3/s (cfs)] is VA}
#'   \item{\emph{n}}{Manning's roughness coefficient (dimensionless)}
#'   \item{\emph{P}}{the wetted perimeters of the channel (m or ft)}
#'   \item{\emph{A}}{the cross-sectional area (m^2 or ft^2)}
#'   \item{\emph{S}}{the slope of the channel bed (m/m or ft/ft)}
#'   \item{\emph{K_n}}{the conversion constant -- 1.0 for SI and
#'        3.2808399 ^ (1 / 3) for English units -- m^(1/3)/s or ft^(1/3)/s}
#' }
#'
#'
#'
#'
#' Other important equations regarding the parabolic cross-section follow:
#' \deqn{R = \frac{A}{P}}
#'
#' \describe{
#'   \item{\emph{R}}{the hydraulic radius (m or ft)}
#'   \item{\emph{A}}{the cross-sectional area (m^2 or ft^2)}
#'   \item{\emph{P}}{the wetted perimeters of the channel (m or ft)}
#' }
#'
#'
#' \deqn{A = \frac{2}{3}By}
#'
#' \describe{
#'   \item{\emph{A}}{the cross-sectional area (m^2 or ft^2)}
#'   \item{\emph{y}}{the flow depth (normal depth in this function) [m or ft]}
#'   \item{\emph{B}}{the top width of the channel (m or ft)}
#' }
#'
#'
#'
#' \deqn{P = \left(\frac{B}{2}\right)\left[\sqrt{\left(1 + x^2\right)} + \left(\frac{1}{x}\right) \ln \left(x + \sqrt{\left(1 + x^2\right)}\right)\right]}
#'
#' \describe{
#'   \item{\emph{P}}{the wetted perimeters of the channel (m or ft)}
#'   \item{\emph{B}}{the top width of the channel (m or ft)}
#'   \item{\emph{x}}{4y/b (dimensionless)}
#' }
#'
#'
#'
#'
#' \deqn{x = \frac{4y}{B}}
#'
#' \describe{
#'   \item{\emph{x}}{4y/b (dimensionless)}
#'   \item{\emph{B}}{the top width of the channel (m or ft)}
#'   \item{\emph{y}}{the flow depth (normal depth in this function) [m or ft]}
#' }
#'
#'
#'
#'
#' \deqn{B = B_1 \left(\sqrt{\frac{y}{y_1}}\right)}
#'
#' \describe{
#'   \item{\emph{B}}{the top width of the channel (m or ft)}
#'   \item{\emph{y}}{the flow depth (normal depth in this function) [m or ft]}
#'   \item{\emph{\eqn{B_1}}}{the "bank-full width" (m or ft)}
#'   \item{\emph{\eqn{y_1}}}{the "bank-full depth" (m or ft)}
#' }
#'
#'
#'
#' \deqn{D = \frac{A}{B}}
#'
#' \describe{
#'   \item{\emph{D}}{the hydraulic depth (m or ft)}
#'   \item{\emph{A}}{the cross-sectional area (m^2 or ft^2)}
#'   \item{\emph{B}}{the top width of the channel (m or ft)}
#' }
#'
#'
#'
#' \deqn{Z = \frac{\sqrt{2}}{2}my^2.5}
#'
#' \describe{
#'   \item{\emph{Z}}{the Section factor (m or ft)}
#'   \item{\emph{y}}{the flow depth (normal depth in this function) [m or ft]}
#'   \item{\emph{m}}{the horizontal side slope}
#' }
#'
#'
#'
#' \deqn{E = y + \frac{Q^2}{2gA^2}}
#'
#' \describe{
#'   \item{\emph{E}}{the Specific Energy (m or ft)}
#'   \item{\emph{Q}}{the discharge [m^3/s or ft^3/s (cfs)] is VA}
#'   \item{\emph{g}}{gravitational acceleration (m/s^2 or ft/sec^2)}
#'   \item{\emph{A}}{the cross-sectional area (m^2 or ft^2)}
#'   \item{\emph{y}}{the flow depth (normal depth in this function) [m or ft]}
#' }
#'
#'
#'
#'
#' \deqn{VH = \frac{V^2}{2g}}
#'
#' \describe{
#'   \item{\emph{VH}}{the Velocity Head (m or ft)}
#'   \item{\emph{V}}{the velocity (m/s or ft/s)}
#'   \item{\emph{g}}{gravitational acceleration (m/s^2 or ft/sec^2)}
#' }
#'
#'
#'
#' A rough turbulent zone check is performed on the water flowing in the
#' channel using the Reynolds number (Re). The Re equation follows:
#'
#' \deqn{Re = \frac{\\rho RV}{\\mu}}
#'
#' \describe{
#'   \item{\emph{Re}}{Reynolds number (dimensionless)}
#'   \item{\emph{\\rho}}{density (kg/m^3 or slug/ft^3)}
#'   \item{\emph{R}}{the hydraulic radius (m or ft)}
#'   \item{\emph{V}}{the velocity (m/s or ft/s)}
#'   \item{\emph{\\mu}}{dynamic viscosity (* 10^-3 kg/m*s or * 10^-5 lb*s/ft^2)}
#' }
#'
#'
#'
#' A critical flow check is performed on the water flowing in the channel
#' using the Froude number (Fr). The Fr equation follows:
#'
#' \deqn{Fr = \frac{V}{\left(\sqrt{g * D}\right)}}
#'
#' \describe{
#'   \item{\emph{Fr}}{the Froude number (dimensionless)}
#'   \item{\emph{V}}{the velocity (m/s or ft/s)}
#'   \item{\emph{g}}{gravitational acceleration (m/s^2 or ft/sec^2)}
#'   \item{\emph{D}}{the hydraulic depth (m or ft)}
#' }
#'
#'
#'
#' @note
#' Assumptions: uniform flow, prismatic channel, and surface water temperature
#' of 20 degrees Celsius (68 degrees Fahrenheit) at atmospheric pressure
#'
#' Note: Units must be consistent
#'
#'
#'
#'
#' @param Q numeric vector that contains the discharge value (m^3/s or ft^3/s),
#'   if known.
#' @param n numeric vector that contains the Manning's roughness coefficient n,
#'   if known.
#' @param m numeric vector that contains the "cross-sectional side slope of m:1
#'   (horizontal:vertical)", if known.
#' @param Sf numeric vector that contains the bed slope (m/m or ft/ft),
#'   if known.
#' @param y numeric vector that contains the flow depth (m or ft), if known.
#' @param B1 numeric vector that contains the "bank-full width", if known.
#' @param y1 numeric vector that contains the "bank-full depth", if known.
#' @param Temp numeric vector that contains the temperature (degrees C or degrees
#'   Fahrenheit), if known.
#' @param units character vector that contains the system of units [options are
#'   \code{SI} for International System of Units or \code{Eng} for English units
#'   (United States Customary System in the United States and Imperial Units in
#'   the United Kingdom)]
#'   
#'   
#' @return the missing parameter (Q, n, m, Sf, B1, y1, or y) & area (A), wetted
#'   perimeter (P), velocity (V), top width (B), hydraulic radius (R),
#'   Reynolds number (Re), and Froude number (Fr) as a \code{\link[base]{list}}.
#'
#'
#'
#'
#' @references
#' \enumerate{
#'    \item Terry W. Sturm, \emph{Open Channel Hydraulics}, 2nd Edition, New York City, New York: The McGraw-Hill Companies, Inc., 2010, page 2, 8, 36, 102, 120, 153.
#'    \item Dan Moore, P.E., NRCS Water Quality and Quantity Technology Development Team, Portland Oregon, "Using Mannings Equation with Natural Streams", August 2011, \url{https://web.archive.org/web/20210416091858/https://www.wcc.nrcs.usda.gov/ftpref/wntsc/H&H/xsec/manningsNaturally.pdf}. Retrieved thanks to the Internet Archive: Wayback Machine
#'    \item Gilberto E. Urroz, Utah State University Civil and Environmental Engineering - OCW, CEE6510 - Numerical Methods in Civil Engineering, Spring 2006 (2006). Course 3. "Solving selected equations and systems of equations in hydraulics using Matlab", August/September 2004, \url{https://digitalcommons.usu.edu/ocw_cee/3/}.
#'    \item Tyler G. Hicks, P.E., \emph{Civil Engineering Formulas: Pocket Guide}, 2nd Edition, New York City, New York: The McGraw-Hill Companies, Inc., 2002, page 423, 425.
#'    \item Wikimedia Foundation, Inc. Wikipedia, 26 November 2015, "Manning formula", \url{https://en.wikipedia.org/wiki/Manning_formula}.
#'    \item John C. Crittenden, R. Rhodes Trussell, David W. Hand, Kerry J. Howe, George Tchobanoglous, \emph{MWH's Water Treatment: Principles and Design}, Third Edition, Hoboken, New Jersey: John Wiley & Sons, Inc., 2012, page 1861-1862.
#'    \item Andrew Chadwick, John Morfett and Martin Borthwick, \emph{Hydraulics in Civil and Environmental Engineering}, Fourth Edition, New York City, New York: Spon Press, Inc., 2004, page 133.
#'    \item Robert L. Mott and Joseph A. Untener, \emph{Applied Fluid Mechanics}, Seventh Edition, New York City, New York: Pearson, 2015, page 376.
#'    \item Ven Te Chow, Ph.D., \emph{Open-Channel Hydraulics}, McGraw-Hill Classic Textbook Reissue, New York City, New York: McGraw-Hill Book Company, 1988, pages 21, 40-41.
#'    \item Gary P. Merkley, "BIE6300 - Irrigation & Conveyance Control Systems, Spring 2004", 2004, Biological and Irrigation Engineering - OCW. Course 2, \url{https://digitalcommons.usu.edu/ocw_bie/2/}.
#'    \item The NIST Reference on Constants, Units, and Uncertainty, Fundamental Constants Data Center of the NIST Physical Measurement Laboratory, "standard acceleration of gravity g_n", \url{https://physics.nist.gov/cgi-bin/cuu/Value?gn}.
#'    \item Wikimedia Foundation, Inc. Wikipedia, 15 May 2019, "Conversion of units", \url{https://en.wikipedia.org/wiki/Conversion_of_units}.
#' }
#'
#'
#'
#'
#' @author Irucka Embry
#'
#'
#'
#' @encoding UTF-8
#'
#'
#'
#'
#'
#'
#'
#'
#'
#' @seealso \code{\link{Manningtrap}} for a trapezoidal cross-section, \code{\link{Manningrect}} for a
#'   rectangular cross-section, \code{\link{Manningtri}} for a triangular cross-section,
#'   and \code{\link{Manningcirc}} for a circular cross-section.
#'
#'
#'
#'
#' @examples
#'
#' library(iemisc)
#'
#' # Exercise 4.3 from Sturm (page 153)
#' y <- Manningpara(Q = 12.0, B1 = 10, y1 = 2.0, Sf = 0.005, n = 0.05, units = "SI")
#' # defines all list values within the object named y
#' # Q = 12.0 m^3/s, B1 = 10 m, y1 = 2.0 m, Sf = 0.005 m/m, n = 0.05, units = SI units
#' # This will solve for y since it is missing and y will be in m
#'
#' y$y # gives the value of y
#'
#'
#'
#' # Modified Exercise 4.3 from Sturm (page 153)
#' Manningpara(y = y$y, B1 = 10, y1 = 2.0, Sf = 0.005, n = 0.05, units = "SI")
#' # y = 1.254427 m, B1 = 10 m, y1 = 2.0 m, Sf = 0.005 m/m, n = 0.05, units = SI units
#' # This will solve for Q since it is missing and Q will be in m^3/s
#'
#'
#' 
#' @importFrom fpCompare %==%
#' @importFrom assertthat assert_that
#' @importFrom checkmate qtest
#' @importFrom units set_units make_units drop_units
#' @importFrom round round_r3
#' @importFrom stats uniroot
#'
#' @export
Manningpara <- function (Q = NULL, n = NULL, m = NULL, Sf = NULL, y = NULL, B1 = NULL, y1 = NULL, Temp = NULL, units = c("SI", "Eng")) {

checks <- c(Q, n, m, Sf, y, B1, y1)

units <- units



# Check
assert_that(!any(qtest(checks, "N+(0,)") == FALSE), msg = "Either Q, n, m, Sf, y, B1, or y1 is 0, NA, NaN, Inf, -Inf, empty, or a string. Please try again.")
# only process with finite values and provide an error message if the check fails

assert_that(qtest(units, "S==1"), msg = "There is not an unit type or more than 1 unit type. Please specify either 'SI' or 'Eng'.")
# only process with enough known variables and provide an error message if the check fails

assert_that(isTRUE(any(c("SI", "Eng") %in% units)), msg = "The unit system has not been identified correctly as either 'SI' or 'Eng'. Please try again.")
# only process with a specified unit and provide a stop warning if not







# units
if (units == "SI") {

# use the temperature to determine the density & absolute and kinematic viscosities
Temp <- ifelse(is.null(Temp), 20, Temp) # degrees C

assert_that(qtest(Temp, "N+(0,)"), msg = "Either Temp is equal to or less than 0 C, NA, NaN, Inf, -Inf, empty, or a string. The equation is valid only for temperatures greater than 0 C / 32 F. Please try again.")
# only process with specified, finite values and provide an error message if the check fails

# create a numeric vector with the units of degrees Celsius
T_C <- set_units(Temp, "degree_C")


# create a numeric vector to convert from degrees Celsius to Kelvin
T_K <- T_C


# create a numeric vector with the units of Kelvin
units(T_K) <- make_units(K)


# create viscosities based on temperatures
# saturated liquid density at given temperature in degrees Celsius (SI units)
rho_SI <- density_water(Temp, "SI")

# absolute or dynamic viscosity at given temperature in degrees Celsius and density of rho (SI units)
mu_SI <- dyn_visc_water(Temp, "SI")

# kinematic viscosity at given temperature in degrees Celsius and density of rho (SI units)
nu_SI <- kin_visc_water(rho_SI, mu_SI, rho_units = "kg/m^3", mu_units = "Pa*s or kg/m/s")


k <- 1

g <- 9.80665 # m / s^2

rho <- rho_SI

nu <- nu_SI

mu <- mu_SI

# unit weight of water at given temperature in degrees Celsius and density of rho (SI units)
gamma <- unit_wt(rho = rho, units = "SI")


density_water_units <- "kg/m^3"

dyn_visc_water_units <- "Pa * s or kg/m*s"

kin_visc_water_units <- "m^2/s"


result_units <- c("m", "m^2", "m", "m", "m", "m", "m", "m/s", "m^3/s", "dimensionless", "m/m", "degrees Celsius", "Kelvin", density_water_units, dyn_visc_water_units, kin_visc_water_units, "dimensionless", "dimensionless", "m/m", "m/m", "m/m", "m", "m", "m", "m", "m^3/s", "m", "m", "pascal (N/m^2)", "pascal (N/m^2)")


} else if (units == "Eng") {

# use the temperature to determine the density & absolute and kinematic viscosities
Temp <- ifelse(is.null(Temp), 68, Temp) # degrees F

assert_that(qtest(Temp, "N+(32,)"), msg = "Either Temp is equal to or less than 32 F, NA, NaN, Inf, -Inf, empty, or a string. The equation is valid only for temperatures greater than 32 F / 0 C. Please try again.")
# only process with specified, finite values and provide an error message if the check fails

T_F <- Temp

# create a numeric vector with the units of degrees Fahrenheit
T_F <- set_units(T_F, "degree_F")


# create a numeric vector to convert from degrees Fahrenheit to Kelvin
T_K <- T_F


# create a numeric vector with the units of Kelvin
units(T_K) <- make_units(K)


# create viscosities based on temperatures
# saturated liquid density at given temperature in degrees Fahrenheit (US Customary units)
rho_Eng <- density_water(Temp, "Eng", Eng_units = "slug/ft^3")


# absolute or dynamic viscosity at given temperature in degrees Fahrenheit and density of rho (US Customary units)
mu_Eng <- dyn_visc_water(Temp, "Eng", Eng_units = "slug/ft/s")


# kinematic viscosity at given temperature in degrees Fahrenheit and density of rho (US Customary units)
nu_Eng <- kin_visc_water(rho_Eng, mu_Eng, rho_units = "slug/ft^3", mu_units = "slug/ft/s")



k <- 3.2808399 ^ (1 / 3)

g <- 9.80665 * (3937 / 1200) # ft / sec^2

gc <- 9.80665 * (3937 / 1200) # lbm-ft/lbf-sec^2

rho <- rho_Eng

nu <- nu_Eng

mu <- mu_Eng

# unit weight of water at given temperature in degrees Fahrenheit and density of rho (US Customary units)
gamma <- unit_wt(rho = rho, units = "Eng", Eng_units = "slug/ft^3")


density_water_units <- "slug/ft^3"

dyn_visc_water_units <- "slug/ft*s"

kin_visc_water_units <- "ft^2/s"


result_units <- c("ft", "ft^2", "ft", "ft", "ft", "ft", "ft", "ft/sec (fps)", "ft^3/sec (cfs)", "dimensionless", "ft/ft", "degrees Fahrenheit", "Kelvin", density_water_units, dyn_visc_water_units, kin_visc_water_units, "dimensionless", "dimensionless", "ft/ft", "ft/ft", "ft/ft", "ft", "ft", "ft", "ft", "ft^3/sec (cfs)", "ft", "ft", "lb/ft^2", "lb/ft^2")

}






if (missing(Q)) {

Qfun <- function(Q) {Q - (((((2 / 3) * B1 * sqrt(y / y1) * y)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((B1 * sqrt(y / y1) / 2) * (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2) + (1 / (4 * y) / B1 * sqrt(y / y1)) * log((4 * y) / B1 * sqrt(y / y1) + (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2))))) ^ (2 / 3))}

Quse <- uniroot(Qfun, interval = c(0.0000001, 200), extendInt = "yes")

Q <- Quse$root

B <- B1 * sqrt(y / y1)
x <- (4 * y) / B
A <- (2 / 3) * B * y
P <- (B / 2) * (sqrt(1 + x ^ 2) + (1 / x) * log(x + (sqrt(1 + x ^ 2))))
R <- A / P
D <- A / B

V <- Q / A

w <- y * sqrt(m ^ 2 + 1)

w1 <- NA_real_

w2 <- NA_real_

m1 <- NA_real_

m2 <- NA_real_

Z <- A * R ^ (2 / 3)

E <- y + ((Q ^ 2) / (2 * g * A ^ 2))

Vel_Head <- (V ^ 2) / (2 * g)

Re <- (rho * R * V) / mu

# average shear stress [mean boundary shear stress]
tau0 <- gamma * R * Sf

# maximum shear stress [shear stress in channel at maximum depth]
taud <- gamma * y * Sf

# conveyance
K <- (k / n) * (A * R ^ (2 / 3))


if (Re > 2000) {

cat("\nFlow IS in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is acceptable to use.\n\n")

} else {

cat("\nFlow is NOT in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is not acceptable to use.\n\n")

}

Fr <- V / (sqrt(g * D))

if (Fr %==% 1) {

cat("\nThis is critical flow.\n\n")

} else if (Fr < 1) {

cat("\nThis is subcritical flow.\n\n")

} else if (Fr > 1) {

cat("\nThis is supercritical flow.\n\n")

}

return(list(Q = Q, V = V, A = A, P = P, R = R, B = B, D = D, Re = Re, Fr = Fr))


} else if (missing(n)) {

nfun <- function(n) {Q - (((((2 / 3) * B1 * sqrt(y / y1) * y)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((B1 * sqrt(y / y1) / 2) * (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2) + (1 / (4 * y) / B1 * sqrt(y / y1)) * log((4 * y) / B1 * sqrt(y / y1) + (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2))))) ^ (2 / 3))}

nuse <- uniroot(nfun, interval = c(0.0000001, 200), extendInt = "yes")

n <- nuse$root

B <- B1 * sqrt(y / y1)
x <- (4 * y) / B
A <- (2 / 3) * B * y
P <- (B / 2) * (sqrt(1 + x ^ 2) + (1 / x) * log(x + (sqrt(1 + x ^ 2))))
R <- A / P
D <- A / B

V <- Q / A

Re <- (rho * R * V) / mu

if (Re > 2000) {

cat("\nFlow IS in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is acceptable to use.\n\n")

} else {

cat("\nFlow is NOT in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is not acceptable to use.\n\n")

}

Fr <- V / (sqrt(g * D))

if (Fr %==% 1) {

cat("\nThis is critical flow.\n\n")

} else if (Fr < 1) {

cat("\nThis is subcritical flow.\n\n")

} else if (Fr > 1) {

cat("\nThis is supercritical flow.\n\n")

}

return(list(n = n, V = V, A = A, P = P, R = R, B = B, D = D, Re = Re, Fr = Fr))


} else if (missing(y1)) {

y1fun <- function(y1) {Q - (((((2 / 3) * B1 * sqrt(y / y1) * y)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((B1 * sqrt(y / y1) / 2) * (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2) + (1 / (4 * y) / B1 * sqrt(y / y1)) * log((4 * y) / B1 * sqrt(y / y1) + (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2))))) ^ (2 / 3))}

y1use <- uniroot(y1fun, interval = c(0.0000001, 200), extendInt = "yes")

y1 <- y1use$root

B <- B1 * sqrt(y / y1)
x <- (4 * y) / B
A <- (2 / 3) * B * y
P <- (B / 2) * (sqrt(1 + x ^ 2) + (1 / x) * log(x + (sqrt(1 + x ^ 2))))
R <- A / P
D <- A / B

V <- Q / A

Re <- (rho * R * V) / mu

if (Re > 2000) {

cat("\nFlow IS in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is acceptable to use.\n\n")

} else {

cat("\nFlow is NOT in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is not acceptable to use.\n\n")

}

Fr <- V / (sqrt(g * D))

if (Fr %==% 1) {

cat("\nThis is critical flow.\n\n")

} else if (Fr < 1) {

cat("\nThis is subcritical flow.\n\n")

} else if (Fr > 1) {

cat("\nThis is supercritical flow.\n\n")

}

return(list(y1 = y1, V = V, A = A, P = P, R = R, B = B, D = D, Re = Re, Fr = Fr))


} else if (missing(B1)) {

B1fun <- function(B1) {Q - (((((2 / 3) * B1 * sqrt(y / y1) * y)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((B1 * sqrt(y / y1) / 2) * (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2) + (1 / (4 * y) / B1 * sqrt(y / y1)) * log((4 * y) / B1 * sqrt(y / y1) + (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2))))) ^ (2 / 3))}

B1use <- uniroot(B1fun, interval = c(0.0000001, 200), extendInt = "yes")

B1 <- B1use$root

B <- B1 * sqrt(y / y1)
x <- (4 * y) / B
A <- (2 / 3) * B * y
P <- (B / 2) * (sqrt(1 + x ^ 2) + (1 / x) * log(x + (sqrt(1 + x ^ 2))))
R <- A / P
D <- A / B

V <- Q / A

Re <- (rho * R * V) / mu

if (Re > 2000) {

cat("\nFlow IS in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is acceptable to use.\n\n")

} else {

cat("\nFlow is NOT in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is not acceptable to use.\n\n")

}

Fr <- V / (sqrt(g * D))

if (Fr %==% 1) {

cat("\nThis is critical flow.\n\n")

} else if (Fr < 1) {

cat("\nThis is subcritical flow.\n\n")

} else if (Fr > 1) {

cat("\nThis is supercritical flow.\n\n")

}

return(list(B1 = B1, V = V, A = A, P = P, R = R, B = B, D = D, Re = Re, Fr = Fr))


} else if (missing(y)) {

yfun <- function(y) {Q - (((((2 / 3) * B1 * sqrt(y / y1) * y)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((B1 * sqrt(y / y1) / 2) * (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2) + (1 / (4 * y) / B1 * sqrt(y / y1)) * log((4 * y) / B1 * sqrt(y / y1) + (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2))))) ^ (2 / 3))}

yuse <- uniroot(yfun, interval = c(0.0000001, 200), extendInt = "yes")

y <- yuse$root

B <- B1 * sqrt(y / y1)
x <- (4 * y) / B
A <- (2 / 3) * B * y
P <- (B / 2) * (sqrt(1 + x ^ 2) + (1 / x) * log(x + (sqrt(1 + x ^ 2))))
R <- A / P
D <- A / B

V <- Q / A

Re <- (rho * R * V) / mu

if (Re > 2000) {

cat("\nFlow IS in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is acceptable to use.\n\n")

} else {

cat("\nFlow is NOT in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is not acceptable to use.\n\n")

}

Fr <- V / (sqrt(g * D))

if (Fr %==% 1) {

cat("\nThis is critical flow.\n\n")

} else if (Fr < 1) {

cat("\nThis is subcritical flow.\n\n")

} else if (Fr > 1) {

cat("\nThis is supercritical flow.\n\n")

}

return(list(y = y, V = V, A = A, P = P, R = R, B = B, D = D, Re = Re, Fr = Fr))


} else if (missing(Sf)) {

Sffun <- function(Sf) {Q - (((((2 / 3) * B1 * sqrt(y / y1) * y)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((B1 * sqrt(y / y1) / 2) * (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2) + (1 / (4 * y) / B1 * sqrt(y / y1)) * log((4 * y) / B1 * sqrt(y / y1) + (sqrt(1 + (4 * y) / B1 * sqrt(y / y1) ^ 2))))) ^ (2 / 3))}

Sfuse <- uniroot(Sffun, interval = c(0.0000001, 200), extendInt = "yes")

Sf <- Sfuse$root

B <- B1 * sqrt(y / y1)
x <- (4 * y) / B
A <- (2 / 3) * B * y
P <- (B / 2) * (sqrt(1 + x ^ 2) + (1 / x) * log(x + (sqrt(1 + x ^ 2))))
R <- A / P
D <- A / B

V <- Q / A

Re <- (rho * R * V) / mu

if (Re > 2000) {

cat("\nFlow IS in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is acceptable to use.\n\n")

} else {

cat("\nFlow is NOT in the rough turbulent zone so the Gauckler-Manning-Strickler equation\n is not acceptable to use.\n\n")

}

Fr <- V / (sqrt(g * D))

if (Fr %==% 1) {

cat("\nThis is critical flow.\n\n")

} else if (Fr < 1) {

cat("\nThis is subcritical flow.\n\n")

} else if (Fr > 1) {

cat("\nThis is supercritical flow.\n\n")

}

return(list(Sf = Sf, V = V, A = A, P = P, R = R, B = B, D = D, Re = Re, Fr = Fr))
}
}