f_dst_norm_ppf_core Function

elemental function f_dst_norm_ppf_core(p, mu, sigma) result(x)

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

! ==== Declarations
  real(wp)   , intent(in) :: p                  !! probability between 0.0 - 1.0
  real(wp)   , intent(in) :: mu                 !! distribution location (mean)
  real(wp)   , intent(in) :: sigma              !! distribution dispersion/scale (standard deviation)
  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
  a = -20.0_wp
  b = 20.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_norm_cdf_core(x_mid, mu=0.0_wp, sigma=1.0_wp, tail="left") - p
     ! check if difference is acceptable, update section if not
     if (abs(p_mid) .lt. tol) then
        ! pass final x value and adjust for mu and sigma
        x = mu + sigma * 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_norm_ppf_core