R/windweibull.R

In pablorcw/pwind: Package for wind energy

#' windweibull
#'
#' Weibull distribution function calculates probablity for a given wind speed, a
#' shape factor and a scale factor
#'
#' @param c scale factor
#' @param k shape factor
#' @param ws wind speed
#' @param ro air density
#'
#' @return Function return a list of numeric values of a probablity for a given wind speed.
#' Additionaly for the given c and k the function returns mean wind speed, median wind speed,
#' maximum probability, wind speed for maximum probability, power density, scale factor,
#' shape factor
#' @export
#'
#' @examples
#' weibull(c = 7.0, k = 2.0, ws = 5.0, ro = 1.225)
#' weibull(7,2,5,1.225)
#'
windweibull <- function(x, ...){
  UseMethod("windweibull")
}

#' @rdname windweibull
#'
#' @export
#'
#' @method windweibull default
#'
windweibull.default <- function(c = 7.0, k = 2.0, ws = 5.0, ro = 1.225, ...){
  wfun <- function(cc = c, kk = k, wwss = ws) {
    round((kk/cc)*((wwss/cc)^(kk-1))*exp(-(wwss/cc)^kk),6)
  }
  # probability for the given parametersa and the wind speed
  probability <- round((k/c)*((ws/c)^(k-1))*exp(-(ws/c)^k),6)
  # mean wind speed
  wsmean <- round(c*gamma(1+1/k),2)
  # median wind speed
  wsmedian <- c*log(2)**(1/k)
  # maximum probability for the distribution with the given parameters
  df.w <- data.frame(V=seq(0,25,0.01),P=seq(0,25,0.01))
  df.w[2:2500,2] <- wfun(cc = c, kk = k, wwss=df.w[2:2500,1])
  mp <- max(df.w[2:2500,2])
  # wind speed at the maximum probability for the distribution with the given parameters
  mwsp <- df.w[which.max(df.w[2:2500,2]),1]
  # power density for the the distribution with the given parameters and the air density ro
  pd <- 0.5*ro*c**3*gamma(1+3/k)
  structure(list(probability = probability, mean.wind.speed = wsmean, median.wind.speed = wsmedian,
                 maximum.probability = mp, wind.speed.for.maximum.probability = mwsp,
                 power.density = pd, scale.factor = c, shape.factor = k), class = "windweibull")
}

#' @rdname summary
#'
#' @export
#'
#' @method summary windweibull
#'
summary.windweibull <- function(x,...){
  cat("Weibull distribution for the wind energy\n")
  cat("----------------------------------------\n")
  cat(paste("Scale factor:", format(x$scale.factor, digits = 3, zero.print = TRUE), "\n"))
  cat(paste("Shape factor:", format(x$shape.factor, digits = 3), "\n"))
  cat("Mean wind speed:", format(x$mean.wind.speed,2, digits = 3), "m/s\n")
  cat("Median wind speed:", format(x$median.wind.speed,2, digits = 3), "m/s\n")
  cat("Power density:", format(x$power.density,4, digits = 5), "W/m2\n")
  cat("Maximum probability:", format(x$maximum.probability, digits = 5), "\n")
  cat("Wind speed at the higest power density:", format(x$wind.speed.for.maximum.probability, digits = 3), " m/s\n")
}

#' @rdname mean
#'
#' @export
#'
#' @method mean windweibull
#'
mean.windweibull <- function(x,...){
  x$mean.wind.speed
}

#' @rdname median
#'
#' @export
#'
#' @method median windweibull
#'
median.windweibull <- function(x,...){
  x$median.wind.speed
}

#' @rdname plot
#'
#' @export
#'
#' @method plot windweibull
#'
plot.windweibull <- function(x,...){
  fstep <- 0.1
  rmin <- 0
  rmax <- 30
  df.dim = rmax / fstep
  df.wd <- data.frame(ws = seq(0, rmax, fstep), p = seq(0, rmax, fstep))
  df.wd[2:df.dim, 2] <- windweibull(k = x$shape.factor, c = x$scale.factor, ws = df.wd[2:df.dim, 1])$probability
  s.range <- rmin / fstep
  e.range <- rmax / fstep
  df.wdplot <- as.data.frame(df.wd[s.range:e.range, ])
  plot(df.wdplot$ws, df.wdplot$p,
       xlab = "Wind speed m/s",
       ylab = "Probability",
       ...)
}

#' @rdname hist
#'
#' @export
#'
#' @method hist windweibull
#'
hist.windweibull <- function(x, ...){
  fstep <- 0.1
  rmin <- 0
  rmax <- 30
  df.dim = rmax / fstep
  df.wd <- data.frame(ws = seq(0, rmax, fstep), p = seq(0, rmax, fstep))
  df.wd[2:df.dim, 2] <- windweibull(k = x$shape.factor, c = x$scale.factor, ws = df.wd[2:df.dim, 1])$probability
  s.range <- rmin / fstep
  e.range <- rmax / fstep
  df.wdplot <- as.data.frame(df.wd[s.range:e.range, ])
  df.wdplot$p <- round(df.wdplot$p * 1000,0)
  df.wdplot <- df.wdplot[df.wdplot$p != 0,]
  dfhist <- 0
  dfhist <- unlist(append(dfhist, apply(df.wdplot, 1, function(x) matrix(data = x[1], nrow = x[2])[,1])), use.names = FALSE)
  hist(dfhist,
       freq = FALSE,
       breaks = c(0:30),
       main = "Histogram of Weibull Distribution",
       xlab = "Wind speed m/s",
       ylab = "Probability",
       ...)
}