# Function to generate random observations from a
# hyperbolic distribution using the 
# ratio of uniforms method

rhypRoU <- function(n,Theta){
  hyperbPi <- Theta[1]
  zeta <- Theta[2]
  delta <- Theta[3]
  mu <- Theta[4]

  if(zeta<=0) stop("zeta must be positive")
  if(delta<=0) stop("delta must be positive")
  
  lowerupper <- hyperbCalcRange(c(hyperbPi,zeta,1,0))
  lower <- lowerupper[1]
  upper <- lowerupper[2]
  

  alpha <- zeta * sqrt(1 + hyperbPi^2)
  hyperbBeta <- zeta * hyperbPi

  
  quasiDens <- function(x,alpha,hyperbBeta){
    quasiDens <- exp(-alpha*sqrt(1+x^2)+hyperbBeta*x)
  }
  
  vfun <- function(x,alpha,hyperbBeta){
    vfun <- -alpha*x/sqrt(1+x^2)+hyperbBeta
  }
  
  vplusloc <- uniroot(function(x)vfun(x,alpha,hyperbBeta),
                     lower=lower,upper=upper)$root
  
  vplus <- sqrt(quasiDens(vplusloc,alpha,hyperbBeta))

  ufun <- function(x,alpha,hyperbBeta){
    ufun <- x*exp(-(1/2)*alpha*sqrt(1+x^2)+(1/2)*hyperbBeta*x)
  }

  dufun <- function(x,alpha,hyperbBeta){
    dufun <- (-(1/2)*exp(-alpha*sqrt(1+x^2)/2+hyperbBeta*x/2)*
             (-2*sqrt(1+x^2)+alpha*x*x-x*hyperbBeta*sqrt(1+x^2))/
             sqrt(1+x^2))
  }

  uminusloc <- uniroot(function(x)dufun(x,alpha,hyperbBeta),
                      lower=lower,upper=0)$root
  uplusloc <- uniroot(function(x)dufun(x,alpha,hyperbBeta),
                      lower=0,upper=upper)$root
  
  uminus <- uminusloc*sqrt(quasiDens(uminusloc,alpha,hyperbBeta))
  uplus <- uplusloc*sqrt(quasiDens(uplusloc,alpha,hyperbBeta))
  
  # Make sure the maximum and minimum are the right way around
  if (uminus > uplus){
    temp <- uminus
    uminus <- uplus
    uplus <- temp
  }

  output <- numeric(n)
  
  for(i in 1:n){
    need.value <- TRUE
    while(need.value==TRUE){
      U <- uminus+(uplus-uminus)*runif(1)
      V <- vplus*runif(1)
      X <- U/V
      fX <- quasiDens(X,alpha,hyperbBeta)
      if(V^2<fX){
          need.value <- FALSE
      }
    }
    output[i] <- X
  }
  delta*output+mu
}