#
# Airy function
#
Airy <- function(x){
  result <- x * 0
  for( i in 1:length(x) ){
    if( x[i] >= 0 ){
      result[i] <- 1/pi * sqrt(x[i]/3) * besselK( 2/3 * x[i]^(3/2), 1/3)
    }else{
      result[i] <- sqrt(-x[i])/2 * ( besselJ( 2/3 * (-x[i])^(3/2), 1/3) - 1/sqrt(3) * besselY(2/3 * (-x[i])^(3/2), 1/3) )
    }
  }

  return(result)
}

#
# Airy function
#
Airy <- function(x){
  result <- x * 0
  result[x>0]  <- 1/pi * sqrt(x[x>0]/3) * besselK( 2/3 * x[x>0]^(3/2), 1/3)
  result[x<0]  <- sqrt(-x[x<0])/2 * ( besselJ( 2/3 * (-x[x<0])^(3/2), 1/3) - 1/sqrt(3) * besselY(2/3 * (-x[x<0])^(3/2), 1/3) )
  result[x==0] <- 1.115353525912245791218/pi # please note, that this value is not the analytical one.
  return(result)
}