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)
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=dp), | intent(in) | :: | x |
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) end if end if ! return end function qe_erf