# Roger Palay copyright 2016-02-15
# Saline, MI 48176
#
pblumenkopf <- function( v, lower.tail=TRUE)
{ # compute the cumulative probability 
  # for the Blumenkopf distribution for vaules between 
  # -3.6 and 3.6 inclusive.
  # By default return P(X<v), but if lower.tail=FALSE
  # then return the upper tail
  
  # handle the extreme cases first
  if( v < -3.6 )
  {tail<-0}
  else if (v > 3.6 )
  {tail<-1}
  else
  { E<-121741/159168+16
    tail <- (v^7/7-26*v^5/5+169*v^3/3)/(E*72)+0.5
  }
  if( lower.tail)
  {return(tail)}
  else
  { return(1-tail)}
}

qblumenkopf <- function( q, lower.tail=TRUE)
{ # inverse function of pblumenkopf()
  if( q < 0 )
  { q <- 0 }
  else if ( q > 1 )
  {q <- 1 }
  if( !lower.tail )
  { q <- 1-q}
  # handle extreme cases
  if ( q == 0 ) {return(-3.6)}
  if ( q == 1 ) { return(3.6)}
  low <- -3.6
  high <- 3.6
  error <- 1
  while (error > 0.00001)
  { 
    next_v <- (low+high)/2
    next_p <- pblumenkopf( next_v )
    if ( next_p <= q )
    { 
      low <- next_v
      error <- q - next_p
    }
    if ( next_p >= q )
    { 
      high <- next_v
      error <-  next_p - q
    }     
  }
  return( next_v )
}