qe_erf Function

  1. utility.F90
  2. w90_utility
  3. qe_erf

private function qe_erf(x)

Uses

  • proc~~qe_erf~~UsesGraph proc~qe_erf qe_erf module~w90_constants w90_constants proc~qe_erf->module~w90_constants

Error function - computed from the rational approximations of W. J. Cody, Math. Comp. 22 (1969), pages 631-637.

 for abs(x) le 0.47 erf is calculated directly
 for abs(x) gt 0.47 erf is calculated via erf(x)=1-erfc(x)

Arguments

Type IntentOptional AttributesName
real(kind=dp), intent(in) :: x

Return Value real(kind=dp)


Calls

proc~~qe_erf~~CallsGraph proc~qe_erf qe_erf proc~qe_erfc qe_erfc proc~qe_erf->proc~qe_erfc proc~qe_erfc->proc~qe_erf

Called by

proc~~qe_erf~~CalledByGraph proc~qe_erf qe_erf proc~qe_erfc qe_erfc proc~qe_erf->proc~qe_erfc proc~utility_wgauss utility_wgauss proc~utility_wgauss->proc~qe_erf proc~gauss_freq gauss_freq proc~utility_wgauss->proc~gauss_freq proc~qe_erfc->proc~qe_erf proc~gauss_freq->proc~qe_erfc

Contents

Source Code

qe_erf

Source Code

  function qe_erf(x)
    !---------------------------------------------------------------------
    !
    !! Error function - computed from the rational approximations of
    !! W. J. Cody, Math. Comp. 22 (1969), pages 631-637.
    !!
    !!     for abs(x) le 0.47 erf is calculated directly
    !!     for abs(x) gt 0.47 erf is calculated via erf(x)=1-erfc(x)
    !
    use w90_constants, only: dp

    implicit none
    real(kind=dp), intent(in) :: x
    real(kind=dp) :: x2, p1(4), q1(4)
    real(kind=dp) :: qe_erf
    data p1/2.426679552305318E2_dp, 2.197926161829415E1_dp, &
      6.996383488619136_dp, -3.560984370181538E-2_dp/
    data q1/2.150588758698612E2_dp, 9.116490540451490E1_dp, &
      1.508279763040779E1_dp, 1.000000000000000_dp/
    !
    if (abs(x) > 6.0_dp) then
      !
      !  erf(6)=1-10^(-17) cannot be distinguished from 1
      !
      qe_erf = sign(1.0_dp, x)
    else
      if (abs(x) <= 0.47_dp) then
        x2 = x**2
        qe_erf = x*(p1(1) + x2*(p1(2) + x2*(p1(3) + x2*p1(4)))) &
                 /(q1(1) + x2*(q1(2) + x2*(q1(3) + x2*q1(4))))
      else
        qe_erf = 1.0_dp - qe_erfc(x)
      endif
    endif
    !
    return
  end function qe_erf