/
operpoint.R
80 lines (73 loc) · 2.89 KB
/
operpoint.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#' Uses input pump and system curves to find the operating point
#' for a pump and create a plot.
#'
#' @param pcurve A pumpcurve object
#' @param scurve A systemcurve object
#'
#' @return Returns a list including:
#' \itemize{
#' \item Qop - flow at the operating point [\eqn{m^3 s^{-1}}{m^3/s} or \eqn{ft^3 s^{-1}}{ft^3/s}]
#' \item hop - head at the operating point [\eqn{m}{m} or \eqn{ft}{ft}]
#' \item p - a plot object of the curves
#' }
#'
#' @author Ed Maurer
#'
#' @seealso \code{\link{pumpcurve}} and \code{\link{systemcurve}} for input object preparation
#'
#' @importFrom stats uniroot optimise
#'
#' @name operpoint
NULL
#' @export
#' @rdname operpoint
operpoint <- function (pcurve = NULL, scurve = NULL) {
if(!inherits(pcurve, "pumpcurve")) {
stop("Pump curve not a valid pumpcurve object")
}
if(!inherits(scurve, "systemcurve")) {
stop("System curve not a valid systemcurve object")
}
if(scurve$curve(0) > pcurve$curve(0)) {
stop("System static head exceeds pump shutoff head")
}
if(scurve$units != pcurve$units) {
stop("Unit systems for curves do not match")
}
units <- scurve$units
if (length(units) != 1) stop("Incorrect unit system. Specify either SI or Eng.")
if (units == "SI") {
txtx <- expression("Flow, "~m^{3}/s)
txty <- sprintf("Head, m")
} else if (units == "Eng") {
txtx <- expression("Flow, "~ft^{3}/s)
txty <- sprintf("Head, ft")
}
#find intersection in interval of flows from 0 to somewhere beyond
#point where pump curve reaches zero
maxQ <- uniroot(pcurve$curve, interval = c(0, 100000))$root
maxh <- max(scurve$curve(maxQ), pcurve$curve(maxQ))
Qop = optimise(f=function(x) abs(scurve$curve(x)-pcurve$curve(x)), c(0,maxQ))$minimum
hop = as.numeric(scurve$curve(Qop))
#oppt <- paste0( 'Op pt:(',Qop,',',hop,')')
oppt1 <- paste0( "Q:",signif(Qop,3) )
oppt2 <- paste0( "h:",signif(hop,3) )
#create a plot to export
p <- ggplot2::ggplot() +
ggplot2::scale_x_continuous(txtx, limits = c(0, maxQ), expand = c(0,0)) +
ggplot2::scale_y_continuous(txty, limits = c(0, maxh), expand = c(0,0)) +
ggplot2::geom_point(ggplot2::aes(x=Qop, y=hop)) +
ggplot2::geom_function(fun = pcurve$curve, colour = "blue") +
ggplot2::annotate("text", x=maxQ*0.05,y=maxh*0.05,
label = pcurve$eqn, parse=T, hjust = 0, colour = "blue") +
ggplot2::geom_function(fun = scurve$curve, colour = "black") +
ggplot2::annotate("text", x=maxQ*0.05,y=maxh*0.10,
label = scurve$eqn, parse=T, hjust = 0, colour = "black") +
ggplot2::annotate("text", x=maxQ*0.95,y=hop+maxh*0.03,
label = oppt1, parse=T, hjust = 1, colour = "black") +
ggplot2::annotate("text", x=maxQ*0.95,y=hop-maxh*0.03,
label = oppt2, parse=T, hjust = 1, colour = "black") +
ggplot2::theme_bw()
output <- list(Qop = Qop, hop = hop, p = p)
return(output)
}