f_dst_exp_ppf_core Function

elemental function f_dst_exp_ppf_core(p, lambda, loc) result(x)

! ==== Description
!! Percent point function/quantile function for exponential distribution.

! ==== Declarations
  real(wp)   , intent(in) :: p                  !! probability between 0.0 - 1.0
  real(wp)   , intent(in) :: loc                !! location parameter
  real(wp)   , intent(in) :: lambda             !! lambda parameter, beta(scale) = 1/lambda = mu/mean
  integer(i4), parameter  :: i_max = c_bisect_i !! max. number of iterations
  real(wp)   , parameter  :: tol = c_bisect_tol !! p deviation tolerance
  real(wp)                :: a, b               !! section bounds for bisection algorithm
  real(wp)                :: x_mid, p_mid       !! x and p mid points in bisection algorithm
  integer(i4)             :: i                  !! for iteration
  real(wp)                :: x                  !! sample position

! ==== Instructions

! ---- compute PPF

  ! set initial section (min possible approaches 0)
  a = loc
  b = loc - log(1.0_wp - p) / lambda * 10.0_wp

  ! iteratively refine with bisection method
  do i = 1, i_max
     x_mid = 0.5_wp * (a + b)
     ! difference between passed p and new mid point p
     p_mid = f_dst_exp_cdf_core(x_mid, lambda=lambda, loc=loc, tail="left") - p
     ! check if difference is acceptable, update section if not
     if (abs(p_mid) .lt. tol) then
        ! pass final x value
        x = x_mid
        return
     elseif (p_mid .lt. 0.0_wp) then
        a = x_mid
     else
        b = x_mid
     endif
  enddo

  ! if x not found in iterations, pass sentinel
  if (i .ge. i_max) x = c_sentinel_r

end function f_dst_exp_ppf_core