# Function to generate random observations from a
# hyperbolic distribution using the 
# transformed density rejection method
# Follows Hoermann, ACM ToMS (1995),21,182--193.
# Uses a squeeze to help with test

rhypTDRsq <- 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")
  
  # restrict to delta = 1, mu = 0
  lowerupper <- calculate.range(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)
  }
  
  h <- function(x,alpha,hyperbBeta){
    h <- -alpha*sqrt(1+x^2)+hyperbBeta*x
  }
  
  dh <- function(x,alpha,hyperbBeta){
    h <- -alpha*x/sqrt(1+x^2)+hyperbBeta
  }
  
  xl <- uniroot(function(x) -alpha*sqrt(1+x^2)+hyperbBeta*x+zeta+1,
                     lower=lower,upper=hyperbPi)$root
  xr <- uniroot(function(x) -alpha*sqrt(1+x^2)+hyperbBeta*x+zeta+1,
                     lower=hyperbPi,upper=upper)$root 
                     
  
  hxl <- h(xl,alpha,hyperbBeta)
  hxr <- h(xr,alpha,hyperbBeta)
  dhxl <- dh(xl,alpha,hyperbBeta)
  dhxr <- dh(xr,alpha,hyperbBeta)                 
                     
  bl <- xl+(-zeta-hxl)/dhxl
  br <- xr+(-zeta-hxr)/dhxr
  
  Fhm <- exp(-zeta)
  fm <- Fhm
    
  vl <- Fhm/dhxl
  vc <- fm*(br-bl)
  vr <- -Fhm/dhxr
  
  # define squeeze function
  sl <- (-zeta-hxl)/(hyperbPi-xl)
  sr <- (-zeta-hxr)/(hyperbPi-xr)
  
  output <- numeric(n)
  
  for(i in 1:n){
    need.value <- TRUE
    while(need.value==TRUE){
      U <- (vl+vc+vr)*runif(1) 
      if(U<=vl){
        X <- (log(-U*dhxl+Fhm)-hxl)/dhxl+xl
        lx <- exp(dhxl*(X-xl)+hxl)
      }else{
        if(U<=(vl+vc)){
          X <- ((U-vl)/vc)*(br-bl)+bl
          lx <- fm
        }else{
          # Beware error in Hoermann here
          X <- (log(((U-(vl+vc)))*dhxr+Fhm)-hxr)/dhxr+xr
          lx <- exp(dhxr*(X-xr)+hxr)
        } 
      }
      V <- lx*runif(1)
      fX <- quasiDens(X,alpha,hyperbBeta)
      if(X<hyperbPi){
        if((X>xl)&(V<=exp(-zeta-(hyperbPi-X)*sl))){
          need.value <- FALSE
        }else{
          if(V<=fX){
            need.value <- FALSE
          }
        }
      }else{      #X>=hyperbPi
        if((X<xr)&(V<=exp(-zeta-(hyperbPi-X)*sr))){
          need.value <- FALSE
        }else{
          if(V<=fX){
            need.value <- FALSE
          }
        }
      }
    }
    output[i] <- X
  }  
  delta*output+mu
}