## Probability density function for the Generalised Normal Laplace distribution
dgnl <- function(x, mu = 0, sigma = 1, alpha = 1, beta = 1, rho = 1,
                 param = c(mu, sigma, alpha, beta, rho)) {

  ## check parameters
  parResult <- gnlCheckPars(param)
  case <- parResult$case
  errMessage <- parResult$errMessage

  if (case == "error")
    stop(errMessage)

  mu <- param[1]
  sigma <- param[2]
  alpha <- param[3]
  beta <- param[4]
  rho <- param[5]

  ## Shifting by mu
  x <- x - mu

  ## Initialising result vector
  pdfValues <- rep(0, length(x))

  ## Because 'integrate' doesn't take vectors as input, we need to iterate over
  ## x to evaluate densities
  for (i in 1:length(x)) {
    ## Modified characteristic function. Includes minor calculation regarding
    ## complex numbers to ensure the function returns a real number
    chfn <- function(s) {
      result <- (alpha * beta * exp(-((sigma^2 * s^2) / 2))) /
                (complex(real = alpha, imaginary = -s) *
                 complex(real = beta, imaginary = s))
      result <- result^rho  ## Scaling result by rho
      r <- Mod(result)
      theta <- Arg(result)
      r * cos(theta - (s * x[i]))
    }

    ## Integrating modified characteristic function
    pdfValues[i] <- (1 / pi) * integrate(chfn, 0, Inf)$value
  }

  ## Returning vector of densities
  pdfValues
}